Splash wave and crown breakup after disc impact on a liquid surface
aa r X i v : . [ phy s i c s . f l u - dyn ] F e b Under consideration for publication in J. Fluid Mech. Splash wave and crown breakup after discimpact on a liquid surface.
I V O R. P E T E R S , D E V A R A J v a n d e r M E E R A N D
J. M. G O R D I L L O † Department of Applied Physics and J.M. Burgers Centre for Fluid Dynamics, University ofTwente, P.O. Box 217, 7500 AE Enschede, The Netherlands ´Area de Mec´anica de Fluidos, Departamento de Ingener´ıa Aeroespacial y Mec´anica deFluidos, Universidad de Sevilla, Avenida de los Descubrimientos s/n 41092, Sevilla, Spain.(Received ?? and in revised form ??) In this paper we analyze the impact of a circular disc on a free surface using experiments,potential flow numerical simulations and theory. We focus our attention both on the studyof the generation and possible breakup of the splash wave created after the impact andon the calculation of the force on the disc. We have experimentally found that dropsare only ejected from the rim located at the top part of the splash –giving rise to whatis known as the crown splash– if the impact Weber number exceeds a threshold value We crit ≃ Bo tip based onthe rim deceleration and its radius of curvature, with which we show using both numericalsimulations and experiments that a crown splash only occurs when Bo tip &
1, revealingthat the rim disrupts due to a Rayleigh-Taylor instability. Neglecting the effect of air,we show that the flow in the region close to the disc edge possesses a Weber-number-dependent self-similar structure for every Weber number. From this we demonstrate that Bo tip ∝ We , explaining both why the transition to crown splash can be characterized interms of the impact Weber number and why this transition occurs for W e crit ≃ F D . The existenceof the air layer also limits the range of times in which the self-similar solution is validand, accordingly, the maximum deceleration experienced by the liquid rim, what sets thelength scale of the splash drops ejected when W e > We crit .
1. Introduction
The seemingly straightforward experiment of impacting a solid or liquid object on aliquid surface exhibits many challenges for physical understanding, as was already no-ticed more than a century ago (Worthington & Cole 1896, 1900; Worthington 1908). Thecreation and collapse of an underwater cavity has been studied extensively for a solid im-pacting a liquid surface (Richardson 1948; Gaudet 1998; Lee et al. et al. et al. et al. et al. et al. et al. et al. † Corresponding author: [email protected]
Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo is formed by the liquid that is moving upwards close to the downwards moving object,which initially can be considered as a Wagner problem, under the condition that the ob-ject can locally be approximated as flat (Wagner 1932; Scolan & Korobkin 2001). Whenthe object is not locally flat, a splash will only be formed provided that certain conditionsare fulfilled. These conditions can sometimes be as subtle as the wetting properties of asmooth sphere, as was shown by May (1951); Duez et al. (2007).The impact of a flat plate on a free surface allows for analytical investigations andcan also be studied using experiments and numerical simulations. Self-similar solutionsfor the case without surface tension ( We → ∞ ) have been found by Iafrati & Korobkin(2004) for the initial stage after impact, whose scaling exponents were already noticed byYakimov (1973). Later, this analysis was expanded in order to calculate the hydrodynamicload on a flat surface close to impact by Iafrati & Korobkin (2008, 2011). One of theexisting differences between our study and the above mentioned ones is that we retainsurface tension effects in the analysis. We show that self-similar solutions exists for any value of the Weber number, a fact corroborated with the aid of boundary integral (BI)simulations. We extend our analysis by accounting for air effects, and show that theseenormously influence both the splash wave velocity during the initial instants after impactand the force on the disc. Surprisingly, the gas-to-liquid density ratio also determine thesizes of the drops generated when the impact speed is above a certain threshold. Indeed,the liquid that is thrown upwards due to the impact of an object, develops in a thin sheetwith a cylindrical rim on top of it (Yarin 2006) which is susceptible to instabilities thatmay result in the ejection of droplets. Finding the nature of instabilities that result inthe formation of droplets has been the motivation for many studies, examples of whichare given below.In Krechetnikov & Homsy (2009), it is argued that the crown formation is the re-sult of a Richtmyer-Meshkov instability, and the effect of interfacial curvature on liquidrims is discussed in Krechetnikov (2009), again in the light of Richtmyer-Meshkov andRayleigh-Taylor instabilities. Deegan et al. (2008) showed that in a narrow range of pa-rameters, this instability can develop in a very regular pattern, which was later usedby Zhang et al. (2010) to experimentally show that the wavelength of this regular pat-tern agrees with the predicted value corresponding to a Rayleigh-Plateau instability.In an analytical study, Krechetnikov (2010) elaborates more on the interplay betweenthe Rayleigh-Plateau and the Rayleigh-Taylor instability. Recently, Lister et al. (2011)specifically isolated the Rayleigh-Taylor instability on a cylinder, tuning the body forcenormal to its surface by tilting a liquid cylinder inside a liquid of higher density. Finally,Lhuissier & Villermaux (2012) found that a Rayleigh-Taylor instability is responsible forthe formation of ligaments from the liquid rim observed in bursting bubbles.In this paper we argue that the splash wave deceleration on its top part plays a keyrole in the ejection of droplets. Indeed, we find that drops will be detached when theeffective weight per unit length of the liquid in the rim, ρ ( A tip − g ) R c , overweighs thesurface tension forces per unit length σ , namely, when the local Bond number satisfiesthe condition Bo tip = ρ ( A tip − g ) R C /σ > , (1.1)where A tip is the rim deceleration, ρ and σ are, respectively, the density and surfacetension of water and R C is the radius of curvature of the rim. We combine experiments,numerical simulations and theory to show that the local Bond number at the rim is therelevant parameter to predict the transition to crown splash.The paper is organized as follows. In § plash wave and crown breakup after disc impact on a liquid surface. ( a )( b ) ( c ) Figure 1.
The splash generated at T = 0 .
01 s after the normal impact against a water interfaceof a circular disc with R D = 20 mm. ( a ) A snapshot of the disk just before the impact takesplace. ( b ) The splash created by an impact below the critical Weber number impact velocity V D = 0 . / s; We = 99). No drops are ejected from the tip of the rim since We < We crit . ( c )The splash created by an impact above the critical Weber number ( V D = 1 . / s; We = 274).Since We > We crit , the rim breaks into drops with sizes of the order of the rim width. experimental results can be reproduced by using potential flow numerical simulations.Making use of potential flow theory, we find in § We for dimensionlesstimes satisfying t ≪
1. In section § We crit above which the crown breaks up and droplets are ejected. In § We > We crit . In addition, in this section we study the force exerted on the disc duringthe impact process and show how it is regularized by the air layer in between the diskand the liquid. Conclusions are drawn in §
2. Experiments and comparison with boundary integral simulations.
Our experimental setup consists of a disc with radius R D which we pull down througha water surface at a constant speed V D using a linear motor. The linear motor ensuresthat our disc always moves with a constant prescribed velocity. We record the eventsusing a Photron SA1 high speed camera. A more detailed description of the experimentalsetup can be found in Bergmann et al. (2009). From now on, distances, times, velocitiesand pressures are made dimensionless using R D , R D /V D , V D and ρV D as characteristiclength, time, velocity and pressure respectively. Also, variables in capital letters willdenote dimensional quantities, whereas their corresponding dimensionless counterpartswill be expressed using lower case letters.Figure 1 shows the surface deformation that occurs immediately after a circular discimpacts perpendicularly with a constant velocity onto a free surface bounding a deepwater layer. From these images, it can be appreciated that a circular liquid sheet –thesplash– is ejected out of the liquid bulk all along the perimeter of the disc. The sheet thenpropagates radially outwards, ‘informing’ the rest of the fluid of the solid body impact.The liquid speeds inside the splash are much larger than the disc impact velocity V D :indeed, for a given instant in time, the distance traveled out of the liquid bulk by thetop part of the liquid sheet is much larger than the distance traveled within the liquid bythe disc. We observe that, depending on the impact velocity, the liquid sheet can eitherbreakup into drops or just retract into the liquid bulk without breaking. Indeed, if V D or, equivalently, if the impact Weber number We = ρR D V D /σ is sufficiently large, the Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo rim at the highest part of the sheet breaks into drops, provoking what we refer to as thecrown breakup or crown splash [figure 1( c )]. If, in contrast, V D (or We ) is sufficientlysmall, surface tension forces and gravity pull back the edge of the rim into the liquid andno breakup occurs [figure 1( b )].Note in figures 1–2 that initially the splash wave develops in a confined spatial region,which grows in time. More specifically, this region has a characteristic length X P ( T ) ≪ R D located very close to the disc edge, with X P the distance from the disc edge to thepoint P indicated in figure 2 where the slope in the splash is -1. Figure 3( a ) shows thedistance traveled by the point P in the splash wave as a function of time for several valuesof the impact velocity. Clearly, the propagation velocity of the splash wave increasesas the impact velocity of the disc increases. The different curves in figure 3( a ), whenrepresented using dimensionless units, collapse onto a single one as it is depicted infigure 3( b ), meaning that the dimensionless velocity of the splash wave is independent ofthe impact Weber number. The results in figure 3( b ) also indicate that the speed of thesplash wave is initially constant since ˙ x P ≃ . → ˙ X P ≃ . V D for t ≃
0. However, aftera short time interval, the splash velocity is no longer constant but decreases in time.Due to the fact that the impact Reynolds number Re = V D R D /ν (where ν indicatesthe kinematic viscosity of the liquid) is such that Re & O (10 ), we expect that viscouseffects are confined to a thin boundary layer developing at the disc bottom. Thus, thevelocity field in most of the liquid volume is expected to be described using a velocitypotential. As a consequence of this, the time evolution of the free surface will be correctlypredicted using a potential flow boundary integral method of the type used to describethe collapse of cavities (Longuet-Higgins & Oguz 1995; Gekle et al. et al. V D against a free surface which very well reproduces the experimentonce the origin of times is shifted in time by a quantity t exp . It will be shown belowthat t exp plays the role of a virtual origin in time which accounts for the presence ofair between the bottom of the plate and the free surface, as will be discussed in § t & t exp using a single fluid potential flow description.
3. Theoretical description of the splash
In this section we will provide a theoretical description for the formation of the splashneglecting the effect of air density, which is in part analogous to Iafrati & Korobkin(2004). Let us first notice from figures 1–2 that there exist two well-differentiated spa-tial regions after a solid impacts a free surface. Indeed, using the (dimensional) polarcoordinates (
R, θ ) centered at the disc edge shown in figure 2, we observe that thereexists an inner region ∼ X P ( T ) ≪ R D , where the interface deforms appreciably and anouter region X P ( T ) ≪ L ( T ) ≪ R D where the interface hardly changes from its initialposition. We start with deriving the flow field in the outer region, i.e., at distances r fromthe edge of the disc such that x P ≪ r ≪
1, for times close to the moment of impact( t ≪ r ∼ x P ,where the deformations of the free interface are no longer small. We will show theoreti-cally that this specific type of far field boundary condition leads to the flow field withinthe splash region to be self-similar. We confirm this result by rescaling the splash wave plash wave and crown breakup after disc impact on a liquid surface. R D X P ( T )V D L ( T ) xz q r P X P ( T ) W P P Figure 2.
Schematic drawing of the impact of the disc and generation of the splash wave. Wedefine the origin of our coordinate system at the edge of the disc for both the cartesian andpolar coordinate system. The point P, which we use as a reference point to measure velocitiesand lengths in the splash, is defined where the slope of the splash wave is −
1. The length ofthe splash region, X P ( T ), is the distance between the disc edge and the point P. This variableis used to indicate at which distance from the edge of the disc the free surface has deformedsignificantly. We define the intermediate length L ( T ), such that L ( T ) ≫ X P ( T ), but smallenough when compared to R D to approximate the flow near the disc edge as 2-dimensional, asexplained in § T (s) X P ( m ) We = 134 We = 1097 ( a ) x P = . t ( b ) t x P Figure 3. ( a ) Position of the point P defined in figure 2 as a function of time for several valuesof the disc impact velocity. ( b ) The same curves as in ( a ) but expressed in dimensionless units.Observe that, initially, the speed of the splash wave is constant and proportional to V D since,initially, ˙ x P ≃ .
6. After a short time interval, the wave speed is no longer constant but decreasesin time. profiles computed using the numerical boundary integral method for several values ofthe impact Weber number (Section 3.2). Finally, we show that the expected self-similarscalings related to the vertical positions and vertical velocities are fully recovered afterwe compensate for the downwards motion of the disc, which introduces a non self-similarcorrection (Section 3.3). 3.1.
Flow field
Since x P ( t ) ≪
1, it is possible to define an intermediate length ℓ ( t ) such that x P ( t ) ≪ ℓ ( t ) ≪ ℓ ( t ) ≪
1, the velocity potential φ at the intermediate region Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo T = 0.6 ms T = 4.0 ms T = 2.0 ms T = 6.0 ms T = 8.0 ms T = 10.0 ms1 cm ( ) b ( ) c ( ) d ( ) e ( ) f ( ) a Figure 4.
Six snapshots of the experiment at different times T after disc impact, compared withthe corresponding results from the boundary integral simulations (yellow lines). Disc radius is20 mm, the impact speed is 1 m / s. The agreement between numerical simulations and exper-iments is excellent once two precisions are made: First, in the numerical simulations, the tipof the splash is unstable and therefore breaks earlier than in experiments. For this reason, thesplash in the experiment appears higher than in the simulations. Second, the times correspond-ing to the experimental profiles are those of the numerical ones plus a constant T exp ≃ . can be described using a two dimensional approach, which satisfies the following equation: ∇ φ = ∂ φ∂x + ∂ φ∂z = 0 , (3.1)subject to the following boundary conditions: ∂φ∂z = − z = − t ≃ x < , (3.2)which is the kinematic boundary condition imposed by the downward moving disc, φ ≃ z ≃ x > , (3.3)denoting the dynamic boundary condition at the free surface † , and φ → p x + z → ∞ , (3.4)enforcing the fluid far away from the impact to be at rest. † Equation (3.3) has been obtained by approximating the time-integrated unsteady Bernoulliequation for very small values of t , which for z = 0 and x > φ ≃ φ ( t = 0) − |∇ φ ( t = 0) | t ≃ φ ( t = 0), and taking into account that φ ( t = 0) = 0 for z = 0 and x > plash wave and crown breakup after disc impact on a liquid surface. x, z ) defined in figure 2 and have taken intoaccount the observations in figure 4: that for t ≪ T ≪
20 ms in figure 4), the splashdevelops close to the edge of the disc and that the interface is not appreciably distortedin the intermediate region r ∼ ℓ ( t ). The solution to the system (3.1)-(3.4) can be foundusing standard conformal mapping techniques, yielding the complex flow field dωdζ ≡ ∂φ∂x − i ∂φ∂z = i + i A r − / e − iθ/ (cid:0) re iθ (cid:1) (cid:0) re iθ (cid:1) − / , (3.5)where i = √− A is determined by matching the velocity field given inequation (3.5) with the numerical solution of the tangential velocity field at the disc edge,which takes into account the true, i.e., three-dimensional, geometry of the impactor andthe corresponding flow. Figure 5 shows a comparison for two instants of time, such that t ≪
1, between the numerical calculated velocity field and the one given by equation (3.5),which is valid in the intermediate region r ∼ ℓ ( t ). The agreement is almost perfect whenthe deformation of the free surface is minimal, as in figure 5( a ). Figure 5( b ) shows theinfluence of the splash region r ∼ x P , where there is a clear discrepancy between theapproach of equation (3.5) and the numerical solution. However, looking at the region r ∼ ℓ ( t ), where deformation can be neglected, the agreement is fully recovered.Now, since the solution (3.5) constitutes the outer boundary condition for the velocityfield in the inner region r ∼ x P ( t ) ≪ ℓ ( t ) ≪
1, we are interested in the approximation ofthis equation in the limit r ≪
1, which yields, ∂ φ∂ x ≃ Ar − / sin( θ/ , ∂ φ∂ z ≃ − − Ar − / cos( θ/
2) (3.6)and the potential φ ≃ − z − Ar / sin( θ/ , (3.7)which will serve as a boundary condition to the problem of determining the free surfaceshape and potential in the splash region r ∼ x P ( t ). Notice from the first equation in(3.6) that the modulus of the velocity component tangent to the disc bottom surface( θ = π ) diverge as r − / . We will show that the huge tangential velocities near the discedge, which have to accommodate to the mass of fluid at rest, are responsible for thegeneration of the splash wave. 3.2. Self similarity
In the frame of reference moving at the disc velocity, the splash region r ∼ x P ( t ) ≪ ℓ ( t ) can be described by solving the Laplace equation (3.1) subjected to the followingboundary conditions at a given instant in time: ∂φ∂z = 0 at z = 0 x < , (3.8)which is the kinematic boundary condition at the bottom of the disc, and ∂φ∂t + |∇ φ | κ We + z Fr = 0 at z = f ( x, t ) x > , (3.9)which is the dynamic boundary condition at the free surface, with We = ρV D R D /σ theWeber number and Fr = V D / ( gR D ) the Froude number. These equations need to becomplemented with the far field velocity potential in the region where the interface isvirtually undisturbed, i.e., φ → − Ar / sin( θ/
2) for r → ∞ . (3.10) Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo −0.05 0 0.05−0.05−0.04−0.03−0.02−0.0100.010.02 t = 4.86 10 −8 z ( a ) numerical solutionanalytical solution−0.05 0 0.05−0.05−0.04−0.03−0.02−0.0100.010.02 t = 2.00 10 −4 z x ( b ) numerical solutionanalytical solution Figure 5.
The analytical flow field given in equation (3.5) compares very favorably with thenumerical one once the constant A is set to 0 .
44. ( a ) At an extremely short time after impact,when the free surface has not yet deformed, the analytical solution agrees with the boundaryintegral result in the full inner domain where r ≪
1. ( b ) At a later point in time we observe thatclose to the splash region, where the deformation is appreciable, the analytical flow field deviatesfrom the numerical solution. Away from the splash region, the agreement improves again. The function f ( x, t ) defines the position of the free interface, which satisfies the kinematicboundary condition ∂ f∂ t = ∂φ∂z − ∂φ∂x ∂f∂x at z = f ( x, t ) x > , (3.11)where f ( x, t = 0) = 0 and f ( x → ∞ , t ) → t. (3.12)Note that in the Bernoulli equation (3.9) we have used the interfacial curvature κ = ∂ f∂x (cid:18) ∂f∂x (cid:19) ! − / (3.13) plash wave and crown breakup after disc impact on a liquid surface. p = κ We − . Since there is nocharacteristic length scale in the system of equations (3.1) and (3.8)-(3.11), we expectthe existence of self-similar solutions of the type φ = ( t − t ) β ¯ φ (cid:18) x ( t − t ) α , z ( t − t ) α (cid:19) , (3.14)which we now write as φ = τ β ¯ φ ( χ, η ) , (3.15)with χ = x/τ α , η = z/τ α , and τ = t − t . The dimensionless quantity t can be anarbitrary constant, and the shape of the free surface can be expressed as f ( x, t ) = τ α F ( χ ) . (3.16)From the Bernoulli equation (3.9) we find, by comparing the first two terms, that self-similar solutions can only exist if β = 2 α −
1. Moreover, the matching with the boundarycondition at infinity (3.10) imposes the additional restriction 2 β = α . Combining thesetwo conditions then results in β = 1 / α = 2 /
3. Thus, lengths are expected to scalewith τ / and velocities with τ − / . Indeed, the system of equations that we solve forboth ¯ φ and F reads, with relative errors ∼ O ( τ / ) ≪ ∂ ¯ φ∂χ + ∂ ¯ φ∂η = 0 , (3.17) ∂ ¯ φ∂η = 0 at η = 0 , χ < φ → − A ¯ r / sin( θ/
2) for ¯ r → ∞ , with ¯ r = p η + χ (3.19)13 ¯ φ − (cid:18) χ ∂ ¯ φ∂χ + η ∂ ¯ φ∂η (cid:19) + "(cid:18) ∂ ¯ φ∂χ (cid:19) + (cid:18) ∂ ¯ φ∂η (cid:19) + ¯ κ We + ητ / Fr = 0at η = F ( χ ) , (3.20)with ¯ κ = d Fdχ (cid:18) dFdχ (cid:19) ! − / , and23 F − χ dFdχ = ∂ ¯ φ∂η − dFdχ ∂ ¯ φ∂χ at η = F ( χ ) and F ( χ ) → τ + t τ / for χ → ∞ . (3.21)Note that there are two terms breaking the self-similarity of equations (3.17)-(3.21)namely, the last term of equation (3.20) and the far field boundary condition in equation(3.21) due to the presence of τ in them. There will therefore exist a self-similar solutionwhenever τ & t and both τ / ≪ τ / ≪ Fr / ( ∼ O (1)) are satisfied. In otherwords, assuming that t is sufficiently small, the self similar solution will be valid duringa short time interval after impact. Under the conditions in which the flow field is selfsimilar, the disc has moved down only slightly, the normal velocity boundary conditionis then approximately applied at the undisturbed free surface and gravitational effectscan be neglected.It is however crucial to notice that the Laplace pressure term κ We − in equation(3.20) is self-similar. This is because the result of balancing the Laplace pressure term0 Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo −3.5 −3 −2.5 −2 −1.5 −1 −0.5−2.5−2−1.5−1−0.5 We = 274 ( a )log t slope = 2/3log z p log x p −3 −2 −1 0 1−2−1.5−1−0.500.51 log t − t l og x p We = 274 ( b ) simulation, t = 0experiment, t = 0experiment, t = 0.03 Figure 6. ( a ) Time evolution of the horizontal and vertical coordinates of point P. Numerical re-sults reveal that the scaling for both x p and z p follows the prediction of equation equation (3.14)for α = 2 /
3. As explained in the main text, the results for z p slightly deviate from a pure powerlaw due to the real boundary condition at infinity in (3.10). ( b ) Comparison between the exper-iment and the simulation for the radial position of point P. The green diamonds indicate theunmodified experimental data. Once a time shift i.e, a virtual origin for times t = t exp = 0 . t ∼
1, because of geometrical and gravitationaleffects that break the self-similar scaling. in the Bernoulli equation (3.20) with the inertial terms on the left hand side using theself-similar ansatz (3.15)-(3.16) will give us 2 β = α and β = 1 − α , which is solved byexactly the same exponents α = 2 / β = 1 / We . To check our theoretical predictions, figure 6 represents the function x P ( t ) calculated using the single fluid potential flow numerical code finding that, asexpected, x P ( t ) ∝ t / for t ≪
1. For times t ∼ O (1) the t / self-similar behavior breaksbecause the non self-similar terms in equations (3.17)-(3.21) are no longer negligible. Theexperimental data included in figure 6( b ) shows that the radial position of point P can befitted with a function proportional to ( t − t exp ) / , with t exp ∼ .
03 the virtual originof times, which accounts for physical effects not taken into account in the theoreticalderivation and the numerical simulations.Figure 6( a ) shows that, contrary to x P ( t ), the time evolution of the vertical coordinateof point P calculated numerically, z P ( t ), is not purely proportional to t / . As will beshown in § x P ∝ τ / scaling is obtained by making useof the fact that the boundary condition at infinity (3.19) implies that the disc acts as aconstant ‘source’ of momentum in time. Indeed, making use of the velocity field calcu-lated from the potential (3.19), the variation per unit time of the horizontal momentum plash wave and crown breakup after disc impact on a liquid surface. e x e z n n L S R X P S D Figure 7.
The surfaces Σ R (semicircle of radius L ≪ X P ), Σ D (line of length L located belowthe disc) and the free interface are the boundaries limiting the volume Ω( t ) to which the integralbalance of horizontal momentum is applied. enclosed in the volume Ω defined in figure 7 satisfies the following integral balance: e x · ddt Z Ω ρ v dω = e x · Z ππ ρ vv · n R dθ + ρ V D Z R v dx − Z ππ p n R dθ ! = − ρ A V D R D Z ππ sin( θ/
2) (sin( θ/
2) cos θ − cos( θ/
2) sin θ ) dθ + ρ V D Z R ∂ φ∂x dx ++ 12 ρ A V D R D Z ππ cos θ dθ = ρ π A V D R D + ρ V D Z R ∂φ∂x dx = ρ π A V D R D ++ O ( ρ V D R D ( R /R D ) / ) , (3.22)where n is the outward normal vector to the surface Σ R defined in figure 7, ( e x , e z ) areused to indicate the unit vectors along the coordinates ( x, z ) and use of the steady Euler-Bernoulli equation namely, p + 1 / ρ v · v = C , has been made to express the pressurefield p over the surface Σ R of figure 7 as a function of the velocity field. Notice that thecontribution of the dimensionless flux of momentum through the surface Σ D defined infigure 7 is of the order of ( R /R D ) / ≪ t ) is constant. We could have used this argumentto deduce that there exists a self-similar solution for short times after impact. Indeed,the characteristic length scale of the spatial region where the flow is time-dependent is ∼ X P ( T ) and since there is not a characteristic time scale either, dimensional analysisindicates that velocities within this part of the splash region should scale as v ∝ V D x P /τ .Therefore, e x · ddt Z Ω ρ v dω ∝ ρV D R D x P τ ≃ ρ V D R D π A → x P ∝ τ / , (3.23)where it has been taken into account that the area in the two-dimensional region wherethe flow is unsteady in the sketch of figure 7 is R Ω dω ∼ R D x P , concluding our alternativederivation of the scaling law for x P .Clearly, the solution of the system (3.17)-(3.21) needs to be found numerically. Itis, however, easier to solve the Laplace equation subjected to the unsteady boundaryconditions given by (3.9)-(3.11) and then express the solution in terms of the variables χ , η and F defined in equations (3.15)-(3.16), i.e., by just using the boundary integralmethod.Figure 8 shows the solution from the boundary integral method for one specific Webernumber at three instances in time, which all are in the regime τ ≪ Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo x z We = 274 τ = 0.005 τ = 0.04 τ = 0.1 Figure 8.
Time evolution of the splash region for a given value of the Weber number. Theorigin of the coordinate z is located at the disc position. Notice that the height of the splashwave is much larger than the depth at which the disc is located. Considering a disc with a radius R D = 20 mm and impacting at a speed V D = 1 m / s, the dimensionless times shown in thisfigure, τ = 5 × − , τ = 4 × − and τ = 10 − , correspond to T = 0 . T = 0 . T = 2 ms respectively. to find self similar solutions. In figure 9( c ) we have rescaled the solution of figure 8,according to (3.15). Inspecting all shapes in figure 9( a-d ), we see that indeed we find selfsimilar solutions for a large range of Weber numbers, that only depend on We . Figure 10shows that for We → ∞ the solution becomes independent of We , confirming the resultsof Iafrati & Korobkin (2004).3.3. Correction for non self-similar terms
From figure 6 it is clear that the vertical position z p of point P does not follow the scaling τ / strictly, except in the limit τ →
0. This difference is caused by the downward motionof the disc (or, equivalently, the upward motion of the undisturbed free surface, sincethe problem (3.17-3.21) was formulated in the frame of reference comoving with thedisc), which introduces a non self-similar term in the far field boundary condition ofequation (3.21). Nevertheless, we can compensate for the effect of the downward motionof the disc by introducing a constant B z such that( z p + B z τ ) ∝ τ / . (3.24)A similar correction is needed for the vertical velocity u p at point P:( u p + B u ) ∝ τ − / . (3.25)The inset in figure 11 shows the position, velocity, and width of the splash at point P.All values become independent of time after the proper rescaling and figure 11 showsthe scaling of lengths ( a ) and velocities ( b ) for a wide range of Weber numbers. Asexpected from the previous analysis, the rescaled values become independent of We forlarge Weber numbers. The same holds for the correction terms: B z = 0 . ± .
03 and B u = 0 . ± .
01 defined in equations (3.24)-(3.25) for We & plash wave and crown breakup after disc impact on a liquid surface. Figure 9.
Shape of the splash with all distances rescaled by τ / for four different values of theWeber number We . Each plot contains three instances in time ( τ = 0 . τ = 0 . τ = 0 . −1.5 −1 −0.5 0 0.5 1−1−0.500.511.5 ( x − x p ) / τ ( z − z p ) / τ / τ = 0.1 We = 68 We = 146 We = 3647 We = 91162 Figure 10.
Comparison of the self-similar shape of the splash (i.e., with all length scales rescaledby τ / ) for four different values of the Weber number. Differences in the rescaled shapes areonly observed for the lower values of the Weber number, see for example the shape for We = 68and We = 9 in figure 9. For high Weber numbers ( We & We . Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo We w p / τ x p / τ ( z p + B z τ ) / τ ( u p + B u ) τ τ Figure 11.
The positions x p and z p , the splash width w p , and the vertical velocity u p , fordifferent values of the Weber number. Clearly, all the rescaled values are independent of theWeber number for We & We in this regime (see also figure 10). The inset shows length scalesand velocities in the splash rescaled by τ / and τ − / respectively. After proper rescaling, allthe plotted values become independent of time, which indicates the self-similar behavior. Dueto the boundary condition at the disc, a constant B u is added to the vertical velocity u p and aterm B z times τ is added to the vertical position z p of point P (see text). All data in the insetcorrespond to a simulation at We = 274.
4. Crown breakup transition
We have provided a scaling for the shape of the splash, and have shown that the splashpossesses a self-similar shape for every value of the Weber number. We will now have acloser look at the breakup of the splash into a crown. With the crown breakup, we arereferring to the ejection of droplets from the tip of the splash.In order to determine when the deceleration is large enough to provoke the growth ofperturbations, we define a local Bond number at the tip of the splash using the radius ofcurvature and the acceleration felt by fluid particles in this region of the wave, namely, A tip − g (see equation (1.1)), where A tip = − d U tip d T , (4.1)with − d U tip /d T >
0. However, the determination of the local Bond number at the tipof the splash wave is not an easy task. Indeed, from the numerical simulations Bo tip canbe obtained very accurately, but only in an indirect manner for most Weber numbers.The reason for this is that the tip of the splash is unstable in the numerical simulationsfor We &
30, as is explained in more detail later in this section. Alternatively, the localBond number at the tip can be determined directly from the experiments, but with arelatively large uncertainty. Both the experimental and numerical method of determiningthe local Bond number will be explained now.We obtain experimental values for Bo tip by tracing the tip of the splash, within atime interval T i with a duration of typically 3 ms, in high-speed movies taken at 5400frames per second. A second-order fit to the position data versus time gives us the meandeceleration of the tip A tip . We determine the radius of the rim R C by measuring it plash wave and crown breakup after disc impact on a liquid surface. Figure 12.
Influence of the node density on the tip of the splash, for We = 274. ( a ) Theradius of curvature R C at the tip has a minimum value comparable to the minimum distancebetween the nodes. The radius of curvature increases as the cylindrical rim grows in size, untilthe rim pinches off, and the radius of curvature starts again at its minimum value. This processis sensitive to the node density, but has no influence on the solution away from the tip, and isnot present when We is small enough. ( b ) The cylindrical rim at its maximum size before itpinches off, for a minimum node distance of 0 . c ) The same simulation as ( b ), just after therim has pinched off. The rim is removed from the simulation directly after the pinch-off becauseit does not represent the physical situation (see text). graphically at the beginning and at the end of T i , giving us a minimum and maximumvalue for R C within the time interval. The experimental values for Bo tip are shown infigure 13 in black dots. The error bars are obtained by using the minimum and maximumof R C , because the radius of curvature is the dominant source of uncertainty, due to thesquared appearance of R C in the definition of the Bond number (1.1).In the boundary integral simulations, the radius of curvature and deceleration can bedetermined accurately at any moment of time. The tip of the splash, which is the region inwhich we are interested, is however only stable in a limited range of low Weber numbers.For higher Weber numbers ( We &
30) a neck is formed below the rim, leading to itspinch-off from the rest of the splash wave (see figure 12( b - c )) † . As soon as the pinch-off has occurred, a new neck forms which pinches off a bit later. This induces a series ofpinch-offs in the numerics, which are unphysical for a number of reasons. The first is that,due to axisymmetry, the droplet is tubular, much unlike the droplets that are generatedin experiments. The second is that the details of the pinch-off are strongly dependenton the node density that is used in the BI simulations, as is shown in figure 12, andtherefore needs to be qualified as a numerical artifact associated with neglecting viscouseffects in the computations. Indeed, we analyzed the numerical simulations and realizedthat the driving force leading to the detachment of the toroidal drop from the rest of theliquid mass is precisely the Bernoulli suction effect caused by the liquid acceleration atthe neck region. The series of pinches influence both the length scale and the decelerationof the tip, such that it is not possible to determine the local Bond number at the tip ofthe splash for high Weber numbers in the simulations. Although the tip of the splash isunstable, we can show that the rest of the splash is not influenced by these numericalartifacts. For this reason, we determine the local Bond number at point P, Bo P , usingthe width W p (see figure 2) as the length scale and calculate the deceleration at t he † This pinched-off donut-shaped volume of fluid is subsequently removed from the simulationto prevent it from crossing the free surface in some point which would cause the code to crash. Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo same point, namely, Bo P = ρσ (cid:18) − dU P dT − g (cid:19) W P . (4.2)In figure 13, the different methods of measuring the local Bond number are put together.Using the lower reference point P on the splash shows that for We & We (blue squares, figure 13), and we find that it is constantin time. Comparing Bo P with the experimental values of Bo tip however, shows thatthe local Bond number is overestimated if we use Bo P . But calculating the local Bondnumber in the numerical simulations at different positions on the splash wave resultsonly in a vertical shift of the points in figure 13. Indeed, figure 13 shows that the valueof the local Bond number evaluated at the point Q on the splash wave where the slopeis -2 (thus, point Q is closer to the edge of the rim than point P), is proportional to Bo P defined in equation (4.2). This result is a consequence of the fact shown in figure 11 thatthe dependence of the lengths and velocities in the splash wave are not dependent onthe Weber number for sufficiently large values of this control parameter. This evidenceindicates that we can calculate the local Bond number at different positions by simplymultiplying Bo P by a constant. For Bo tip , we determine this constant to be ∼ . ∼ . Bo tip numerically calculated for low values of the impact Weber number. The local Bondnumber at the tip that we have deduced from Bo P is shown in figure 13 with red circles.Clearly, the transition to breakup into a crown occurs when the local Bond number atthe tip of the splash is of order unity.The proportionality of the local Bond number with the Weber number and its indepen-dence of time is not unexpected, because it results from the self-similarity solution thatbecomes independent of We for large Weber numbers, as shown in figure 11. Using theself-similar scaling for distances, velocities and accelerations, we write the decelerationand typical length scale in dimensional form as: A tip = V D R D C a τ − / (4.3)and R C = R D C R τ / , (4.4)where C a and C R are dimensionless constants, independent of time and Weber numberfor We ≫
1. Substituting (4.3) and (4.4) in (1.1), and using the fact that A tip ≫ g forlarge values of We gives: Bo tip ≃ C a C R ρV D R D σ = C a C R We , (4.5)clearly showing that the local Bond number is independent of time and proportional to We .Using the proportionality Bo tip ∝ We and the crown breakup condition that Bo tip isof order unity around the transition, we can define a condition for the crown breakuptransition based on We that does not involve the local Bond number. Such a condition ispreferable, because whereas the local Bond number is difficult to determine, the Webernumber is directly given by the experimental conditions. There indeed exists one criticalWeber number We crit above which we always observe breakup of the splash into a crown,as can be seen in table 1. The table shows the values of We , Fr and Re at the crown plash wave and crown breakup after disc impact on a liquid surface. log We l og B o l o c a l We crit Bo P , numeric Bo Q , numeric Bo tip , numeric Bo P , numeric Bo tip , experiment w P w Q PQ Figure 13.
The local Bond number measured in different ways as a function of We . The bluesquares are the local Bond numbers determined in the simulations at point P ( Bo P ), the blacktriangles correspond to the local Bond number at the point Q ( Bo Q ) of the inset and red circlesare these same values of Bo P multiplied by 0 .
1. This multiplication factor has been found bymatching to the experimental values of the local Bond number at the tip (black dots). Thegreen diamonds show the local Bond numbers measured at the tip of the splash in the BIsimulations, which is only possible for We .
30 (see text). The shaded area indicates the rangeof experimental conditions in which we observe the crown breakup transition. Note that thistransition occurs when the local Bond number is of order unity (dashed horizontal line). R D (mm) We crit Fr Re
15 145 ± .
77 1 . ·
20 135 ±
19 2 .
51 1 . ·
25 134 ±
11 1 .
59 1 . ·
30 142 ±
19 1 .
18 1 . · Table 1.
Transition to crown breakup of the splash for different disc radii. The value of thecritical Weber number does not differ appreciably with the size of the disc, where the Froudenumber and the Reynolds number at the transition show a clear dependence on R D . breakup transition † . Note that only the Weber number always has approximately thesame value at the transition, showing that the Weber number is indeed the relevantparameter to indicate the transition to crown breakup and droplet ejection. † Clearly, We , Fr and Re can be related after calculating the disc velocity V D using the discradius R D and either of the three dimensionless numbers in table 1. Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo
5. Air cushion effect: force on the disc and drop size
The previous theoretical derivation of the self-similar solution predicts decelerationswithin the splash wave proportional to ∝ t − / , with t the dimensionless time afterimpact. The huge decelerations expected from the analysis for t ≪ constant in time, and not at a diverging velocity proportional to ∝ t − / , as is predictedby the self-similar solution. The origin of the cutoff time or, equivalently, the origin of thecutoff length below which the self similar solution is no longer valid, could in principlelie in any of the following three sources: (i) the fact that the disc edge is rounded, (ii)the effect of liquid viscosity, and (iii) the air layer between the impacting disk and theliquid.First we assess the fact that the disc edge is not perfectly sharp, but possesses a finitecurvature. Since, in view of figure 6 the self-similar solution is only valid for t > t exp ∼ .
03, we can estimate the cut-off length X = x R D ∼ R D t / exp ∼ . ≈ µ m.To evaluate the effect of liquid viscosity we estimate the viscous cutoff time, t v , whosemagnitude is determined by equating the boundary layer thickness, δ , to the width ofthe splash region at the instant of time when the self similar solution begins to be valid,namely, δ ∼ s µR D t v ρ V D ∼ R D t / v → t / v ∼ Re − / D → t v ∼ Re − D ∼ − . (5.1)Clearly, this value is much smaller than t exp and viscosity can therefore not account forthe observed cut-off time.The remaining physical effect is the finiteness of air density, which is known to be re-sponsible for the so called air cushion effect , analysed by Howison et al. (1991); Wilson(1991) for the case of two dimensional geometries. Figure 14a illustrates the disc ap-proaching the free interface forcing the air between disc and liquid surface to flow radi-ally outwards. Using continuity and neglecting compressibility effects, the radial velocityfield of the gas, V g ( R, T ), can be expressed as a function of (dimensional) time T andthe vertical deformation H g ( T ) of the free surface below the disc as − π R (cid:18) V D − d H g d T (cid:19) + 2 π R [ H g + V D ( T s − T )] V g = 0 ⇒ V g = R V D − d H g /dT ) H g + V D ( T s − T ) , (5.2)where we assume that H g in the region below the disc does not depend on R . In (5.2), T s is a constant which will be arbitrarily fixed to T s = R D /V D since its precise valuepossesses a negligible influence on H g ( T ). Using the expression for the flow field givenin (5.2), the gas pressure dependence on R and T can be determined integrating themomentum equation ρ g ∂V g ∂T + ∂∂ R ρ g V g P g ! = 0 , (5.3) plash wave and crown breakup after disc impact on a liquid surface. Figure 14. a) The downward motion of the disc forces the gas to flow radially outwards. Theoverpressure needed to accelerate the gas deforms the interface. b) The linearity of the Laplaceequation permits to express the liquid potential at R = 0 for a generic value of d H/d T asΦ( R = 0 , T ) = − αR D d H/dT , with αR D V D the value of the potential at R = 0 obtained fromthe solution of the Laplace equation subjected to the boundary conditions depicted in this figure.Here, ∂ Φ /∂n indicates the normal velocity at the interface. which yields the following expression for P g ( R, T ): P g − P a = 38 ρ g (cid:0) R D − R (cid:1) (cid:18) V D − d H g /dTH g + V D ( T s − T ) (cid:19) − ρ g R D − R H g + V D ( T s − T ) d H g d T . (5.4)The desired equation for H g ( T ) is determined from the Euler-Bernoulli equation par-ticularized at R = 0, namely, ρ ∂ Φ ∂ T ( R = 0 , T ) + 12 ρ (cid:18) d H g d T (cid:19) + P g ( R = 0) = P a , (5.5)where P a indicates the atmospheric pressure. Indeed, assuming that the normal velocity inthe liquid below the disc does not depend on R , the potential Φ at R = 0 can be expressedas a function of H g since, by virtue of the linearity of the Laplace equation, Φ( R = 0) = − α R D d H g /d T , with α = 2 /π the value of the potential at R = 0 corresponding tothe solution of ∇ Φ = 0 subjected to the boundary conditions sketched in figure 14b(see Iafrati & Korobkin (2011)). Finally, the substitution of equation (5.4) into (5.5)and subsequent non-dimensionalization provides the following differential equation for h = H g /R D : − (cid:18) α + 14 ρ g ρ t s − t + h (cid:19) d hd t + 12 (cid:18) d hd t (cid:19) + 38 ρ g ρ (cid:18) − d h/dth + t s − t (cid:19) = 0 , (5.6)with t = T V D /R D and similarly for t s .Let us point out that the above analysis rests on the assumption that viscous effects inthe gas are unimportant, and this is indeed the case since the characteristic length L e inwhich viscosity transforms a plug flow into a parabolic Poiseuille velocity profile is justthe entrance length in pipe flows, which is such that L e R D ∼ ρ g V g H µ g H R D = Re g H R D , (5.7)0 Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo where V g is given in equation (5.2), Re g = ρ g V D R D /µ g , H = V D ( T S − T ) + H g andwhere use of the continuity equation (5.2) has been made. The estimation in equation(5.7), which is valid while V D − dH g /dT ∼ O ( V D ), suggests to neglect viscous effects when H /R D & Re − g . In our experiments Re g ∼ O (10 ) and, therefore, we can assume thatthe radial velocity profile in the gas is uniform while H /R D & − . The estimationin equation (5.7) also suggests that viscous effects in the gas should be included in theanalytical description if H < − R D . However, when this condition is fulfilled, theliquid normal velocities are already close to that of the impactor, as shown in figure 15.Thus, the gas overpressures associated to gas inertia set the liquid into motion beforegas viscosity needs to be retained into the analytical description of the air cushion effectdescribed above.The integration of equation (5.6), subjected to appropriate initial conditions, provideswith the time evolution for h , dh/dt , the dimensionless gas pressure p ( r = 0) = P g / ( ρ V D )and the height of the gas layer h = t s − t + h shown in figures 15-16. It is depicted infigure 15(a) that the interface is only appreciably accelerated for t ≃ t & t end with t end = 1 .
01, namely, when the disc position is only slightly below thelevel of the undisturbed interface, d h/dt ≃ h ≃
0. Therefore, the liquid normalvelocities below the disc are smaller than the disc impact velocity while t < t end andequal to the disc velocity for t > t end . This fact implies that the speed at which thesplash wave propagates must be, during a time interval t − t s . .
01, smaller than theone predicted by the previous potential flow numerical simulations where the effect of airwas neglected. This is just what figure 4 shows: the splash wave calculated numericallyneeds 0.6 ms less than the real wave to reach the wave radial position at T = 0 . ms ( τ = 0 . t end − t s ≃ .
01, atime interval very similar to the magnitude of the virtual origin in times, t exp = 0 . t exp too large to find a clear trend.Nevertheless, since all the measured values are of the order of 0 .
01, together with thefact that radial position of the splash waves collapse onto a single curve after rescalingwith R D and V D (see figure 3 b ), suggest that t exp is independent of We , in accordancewith our theory.To check the validity of the model presented above, we have carried out two-fluidpotential flow numerical simulations using the numerical code detailed in Gordillo et al. (2007); Gekle & Gordillo (2011). Both the density ratio ρ g /ρ and the impact Webernumber have been varied in wide ranges of values. Figure 15 shows that the time evolutionof the dimensionless height h is almost perfectly reproduced by our theory. The shapeof the entrapped bubble calculated numerically for different values of the density ratio,depicted in figure 16, shows that the interface accelerates beneath the disc uniformlydownwards, except in a very localized region near the disc edge, where the interfaceaccelerates upwards as a consequence of Bernoulli suction (or, equivalently, due to thegrowth of a Kelvin-Helmholtz instability), causing the entrapment of an air pocket. It isalso shown in figure 16 that the maximum height of the gas layer calculated numericallycompares very favorably with the one predicted theoretically in spite of the fact that ourmodel, which was built under the assumption that h does not depend on r , is unableto capture the upwards acceleration of the free interface. Let us point out here that the plash wave and crown breakup after disc impact on a liquid surface. t log h (model)log h (simulation)log dh / dt log p log h t t m D f D t t ( a ) (b) Figure 15. ( a ) Using continuous lines, we plot the time evolutions of the gas layer depth h , theliquid velocity at r = 0, dh/dt , gas pressure at r = 0, p g , and the distance from the disc bottomto the liquid interface, h = t s − t + h for ρ g /ρ = 1 . × − and α = 2 /π . The inset shows thetime evolution of the different variables for instants of time right after the disc bottom is belowthe undisturbed free interface. The discontinuous line represents the function h ( t ) obtainedfrom the two-fluid potential flow code, which is almost perfectly reproduced by the integrationof equation (5.6). ( b ) Time evolution of the functions representing the dimensionless force on thedisc f D (dashed lines) and m D (continuous lines) defined, respectively, in equations (5.9)-(5.10).The insets show that due to the effect of air the force is not infinite any more, but reaches awell-defined maximum slightly after the disc bottom reaches the level of the unperturbed freeinterface and that the maximum value of the dimensionless change in momentum, defined in(5.10), is m D ( t = 1 . ≃ . ∼ O (1). −3.5 −3 −2.5 −2 −1.5−3−2.5−2−1.5−1−0.50 log ρ g / ρ l log h (model)log h (simulation)log dh / dt Figure 16.
Continuous lines indicate the dependence of both h and d h/dt on the density ratio ρ g /ρ at t = 0 .
99 predicted by the integration of equation (5.6). Diamonds indicate the numericalvalues of the air layer thickness at r = 0 when the first contact between the liquid and the discbottom occurs, which occurs near the disc edge, as shown by the different insets. Numericalsimulations reveal that, except close to the disc edge, the air layer depth is mostly uniform. air cushion effect described above shares many similarities with the bubble entrapmentprocesses taking place either during the impact of drops against solid or liquid surfacesThoroddsen et al. (2003); Mandre et al. (2009); Duchemin & Josserand (2011) or in thecollapse of bubbles Gordillo & Fontelos (2007); Gordillo (2008). It should be mentioned2 Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo −3 −2 −1 0 1−2−1.5−1−0.500.51 log t − t l og x p We = 274simulation, t = 0experiment, t = 0 Figure 17.
Comparison between the time evolution of the horizontal position x P of point P for We = 274 in (i) a single-fluid potential-flow numerical simulation in which the liquid normalvelocities beneath the disc is provided by our theoretical prediction (function dh/dt in figure15a, blue circles) and (ii) the experiments (green diamonds). Note that the experimental data isthe bare data of figure 6b without the time shift. The conclusion is that the simulations based onour theoretical analysis excellently reproduces the experiments, with no adjustable constants. here that, while the gas flow is well described by the lubrication approximation in theimpact of drops, viscous effects are negligible in our case.The critical role played by the thin layer of air located beneath the disc on the splashwave speed can be demonstrated in an even more convincing way. To this end, we havecarried out single-fluid potential-flow simulations in which the effect of air is introducedby replacing the constant impact velocity boundary condition by the function dh/dt represented graphically in figure 15a. Figure (17) shows that this new type of numericalsimulations, which indirectly retain the air cushion effect through its influence on thenormal velocity beneath the disc, perfectly reproduces the radial propagation of thesplash wave and, thus, avoids the need of using an artificial time shift to match simulationswith experiments. The good agreement between our analytical model and experimentalresults also confirms that the effect of gas viscosity on the splash wave formation is smallwhen compared to that of gas inertia, as anticipated above.Once the critical role played by air in the impact process has been shown to be verywell described by our theory, we determine the force on the disc, F D ( T ), by integrationof the pressure field given by equation (5.4), to give F D ( T ) = ρ π R D V D f D ( t ) , (5.8) plash wave and crown breakup after disc impact on a liquid surface. f D ( t ) = 116 ρ g ρ " (cid:18) − dh/dth + t s − t (cid:19) − h + t s − t d hdt (5.9)is represented graphically in figure 15(b) for the particular case of ρ g /ρ = 1 . × − .Figure 15(b) shows the remarkable result that f D ( t ) reaches a maximum when the discbottom is below the level of the unperturbed free surface but before the disc bottomtouches the liquid. Let us point out here that the force on an impacting disc was alreadycalculated by Iafrati & Korobkin (2008, 2011), who neglected the effect of air in theanalysis, a starting hypothesis leading to the force on the disc tend to infinite duringthe initial instants of the impact process. The fact that our theory predicts a maximumforce on the disc, indicates that the air cushion effect described above softens the impactprocess, thus avoiding the diverging force deduced from the analysis when air is removed.Our theory also allows to determine the change in the momentum of a disc falling freelyand impacting against a liquid free surface,∆ M D ( T ) = Z T F D ( T ′ ) d T ′ = ρ π R D V D m D ( t ) , (5.10)where the function m D ( t ) defined in equation (5.10) is also represented graphically infigure 15(b). Since m D ( t = 1 . ≃ . V D of a free fallingdisc of density ρ s and thickness H impacting normally onto water-air interface ( ρ g /ρ =1 . × − ) can be estimated as: ρ s π R D H ρ s ∆ V D ≃ . π ρ V D R D → ∆ V D V D ≃ . ρρ S R D H . (5.11)The result expressed by equation (5.11) has been obtained under the assumption thatthe disc velocity is constant and, thus, it can only be strictly applied in the case inwhich the velocity of the solid does not appreciably vary during the impact process.That is, when the increment in the disc momentum given by equation (5.10) is smallwhen compared with the initial momentum. Nevertheless, equation (5.11) can be usedfor qualitative estimation purposes. Doing so, we can conclude that, since for most solids1 < ρ s /ρ <
20, the air cushion effect could reduce the impacting velocity by an amount ∼ O ( V D ) if R D /H > O (10). Let us finally point out that the validity of our study isrestricted to those situations in which the disc impact velocity is sufficiently low so as toneglect gas compressibility effects. 5.1. Drop size
It was demonstrated above that, in the range of times in which the self-similar solutionis valid (i.e., when t exp ≪ t ≪ R C /R D = r tip ∝ ( t − t exp ) / and a tip ∝ ( t − t exp ) − / respectively and, thus, Bo local = ρ V D R D σ a tip r tip ∝ W e. (5.12)Therefore, whenever Bo local > W e > W e crit , it is expected thatwhen the rim decelerates i.e., at any instant of time such that t > t exp , the rim developsa Rayleigh-Taylor instability with a initial length scale d < r tip given by the condition(see Villermaux & Bossa (2011)) ρV D R D a max,tip d σ ≃ → d ≃ W e − / a − / max,tip , (5.13)4 Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo where it has been taken into account that the maximum deceleration at the tip of thesplash wave is such that a max,tip ≫ gR/V D . Equation (5.13) expresses that the charac-teristic length scale of the fingers, d R D , is set by the condition that their weigh per unitlength, namely, ρ V D R D a tip d , equals the surface tension force per unit length, σ . Themaximum deceleration scales as a max,tip ∝ t − / exp ∼ ( t end − t s ) − / and, therefore, theinitial wavelength of maximum growth rate is given by d ∝ t / exp W e − / , (5.14)where use of equation (5.13) has been made. These initial perturbations give rise to theformation of fingers of increasing width and length which break due to the action ofsurface tension forces in a characteristic capillary time given by t cap = V D T cap /R D = V D /R D (cid:0) ρR D d /σ (cid:1) / = W e / d / ∝ W e − / t exp , (5.15)where use of equation (5.14) has been made. Consequently, the breakup time will begiven by t break ∝ ( t exp + t cap ) and the dimensionless size of the drop formed, which isof the order of the width of the rim (see figure the experimental evidence in figure 18a),is given by d drop ∼ w rim ∝ t / break ∝ t / exp (1 + t cap /t exp ) / ∝ t / exp (cid:16) b W e − / (cid:17) / , (5.16)with b an adjustable constant. Figure 18b shows a reasonable agreement between theexperimental drop size measured from the analysis of the images of the type depictedin figure 1 and the drop diameter predicted by equation (5.16). Let us point out thatfigure 18b reveals that the drop size is hardly dependent of the impact Weber number.This fact indicates that the characteristic length scale of the ejected drops is fixed by theair bubble entrapped underneath the disc, being this the physical mechanism underlayingthe experimentally observed cut-off time, t exp .Let us point out that the broad distribution of drop sizes depicted in figure 18b is mostlikely caused by small disturbances of different initial amplitudes distributed unevenlyalong the azimuthal direction. Indeed, a regular periodic structure of drop ejection alongthe rim is expected only when there is an equal amount of noise in all wavelengths because,only in that case, we would observe the growth of just the wavelength selected by theRayleigh-Taylor type of instability, i.e, the one with the largest growth rate. Amongthe reasons that could lead to the initial amplitude of perturbations to depend on thewavelength and on the azimuthal direction are: a slight deviation from axisymmetry,since this will cause the radial component of the gas velocity to vary in the azimuthaldirection, a small tilt of the impacting disc, geometrical imperfections on the disc itselfor the presence of waves before the disc impacts the free surface.
6. Conclusions
We have studied the formation of a splash wave and the transition to the ejection ofdroplets after the impact of a disc on a liquid by using boundary integral simulations,experiments and theoretical analysis. Only by combining these three methods, we havebeen able to analyze the full problem. Although all the different approaches used possesslimitations either in accuracy, stability or mathematical formulation, each of them isable to reveal specific parts of the problem unaccessible by the other methods. Using theoverlapping parts however, we accurately demonstrated the validity of each of the threedifferent approaches used in this study. plash wave and crown breakup after disc impact on a liquid surface. ( a ) log We l og d d r op ( b ) experiment a t (1 + b We −1/4 ) , a = 1, b = 1 a t (1 + b We −1/4 ) , a = 0.25, b = 20 Figure 18. a) Snapshot showing the crown splash and how the edge finding routine identifiesthe ejected drops. Drop diameter is experimentally determined from the analysis of this type ofimages. b) This figure compares the diameter of the ejected drops with the theoretical predictiongiven by equation (5.16) for two different values of b . By approximating the time just after disc impact ( t ≪
1) as a 2-dimensional potentialflow problem and neglecting the effect of air, we have shown that there exist self-similarsolutions of the second kind for any value of the Weber number. These self-similar solu-tions exist because the matching of the inertial terms in the unsteady Bernoulli equationto the far-field boundary condition for the potential gives the same scaling powers as thematching with the surface tension term. When the Weber number is increased to largevalues ( & ∼ t / ), and velocities ( ∼ t − / ) at different positions in the splash. In the ex-periments, the same scaling is found after introducing a small time-shift t exp which wedemonstrated is associated to the finiteness of gas density. Indeed, in spite of its largedensity contrast with water or the solid, the very thin air cushion entrapped between theimpacting solid and the liquid reveals that gas critically affect the velocity of propagationof the splash wave as well as the force on the disc, F D . The presence of air avoids thesingularity in F D predicted when air effects are neglected (see Iafrati & Korobkin (2011))and our theory predicts that the maximum of F D is reached before the disc touches theliquid.We have shown that the transition to droplet ejection (crown breakup transition) iscaused by a Rayleigh-Taylor instability. We experimentally determined that the localBond number based on the tip deceleration and the tip width, is of order unity at thesplash transition. Theoretical analysis predict, and numerical simulations confirm that, Bo tip depends linearly on the Weber number for We & Bo tip and We , we have concluded that the splash transition can be identified bya critical Weber number, which we indeed find in the experiments. We have also foundthat the air cushion effect also sets the maximum deceleration experienced by the rimedge and have used this information to scale the sizes of the drops ejected for those exper-imental conditions in which We > We crit ≃
140 and the impact Weber number is belowthe threshold for which surface sealing occurs. Finally, we have found the interesting re-sult that the thickness of the air layer entrapped beneath the disc sets the characteristic6
Ivo R. Peters, Devaraj van der Meer and Jos´e Manuel Gordillo length scale of the drops ejected when the impact Weber number exceeds the critical one.IRP and DvdM acknowledge financial support by NWO. JMG thanks financial supportby the Spanish Ministry of Education under Project DPI2011-28356-C03-01 and theJunta de Andaluc´ıa under project P08-TEP-03997. The Spanish projects have been partlyfinanced through European funds.
REFERENCESBergmann, R., van der Meer, D., Gekle, S., van der Bos, A. & Lohse, D.
J. Fluid Mech. , 381–409.
Bergmann, Raymond, van der Meer, Devaraj, Stijnman, Mark, Sandtke, Marijn,Prosperetti, Andrea & Lohse, Detlef
Phys. Rev. Lett. (15), 154505. Deegan, R. D., Brunet, P. & Eggers, J.
Nonlinearity (1),C1–C11. Duchemin, L. & Josserand, C.
Phys. Fluids , 091701. Duchemin, Laurent, Popinet, Stephane, Josserand, Christophe & Zaleski, Stephane
Phys. Fluids (9), 3000–3008. Ducleaux, V., Caill´e, F., Duez, C., Ybert, C., Bocquet, L. & Clanet, C.
J. Fluid Mech. , 1–19.
Duez, C., Ybert, C., Clanet, C. & Bocquet, L.
Nature Phys. (3), 180–183. Gaudet, S.
Phys. Fluids (10), 2489. Gekle, S., van der Bos, A., Bergmann, R., van der Meer, D. & Lohse, D.
Phys. Rev.Lett. (8), 084502.
Gekle, Stephan & Gordillo, J. M.
J. Fluid Mech. , 293–330.
Gekle, Stephan & Gordillo, J. M.
Int. J. Numer. Meth. Fluids , 1456–1469. Gekle, Stephan, Gordillo, Jos´e Manuel, van der Meer, Devaraj & Lohse, Detlef
Phys. Rev. Lett. (3), 034502.
Gordillo, J. M.
Phys. Fluids , 112103. Gordillo, J. M. & Fontelos, M. A.
Phys.Rev. Lett. , 144503. Gordillo, J. M. & Gekle, Stephan
J. Fluid Mech. , 331–346.
Gordillo, J. M., Sevilla, A., & Mart´ınez-Baz´an, C.
Phys. Fluids , 077102. Hogrefe, J.E., Peffley, N.L., Goodridge, C.L., Shi, W.T., Hentschel, H.G.E. & Lath-rop, D.P.
Physica D (1-4),183–205.
Howison, S. D., Ockendon, J. R. & Wilson, S. K.
J . Fluid Mech , 215–230.
Iafrati, A. & Korobkin, A.
J. Eng. Math. , 199–224. Iafrati, Alessandro & Korobkin, Alexander A.
Phys. Fluids (7), 2214. Iafrati, Alessandro & Korobkin, Alexander A.
Phys. Fluids (8), 082104. plash wave and crown breakup after disc impact on a liquid surface. Krechetnikov, R.
J. Fluid Mech. , 387.
Krechetnikov, R.
Phys. Fluids (9), 092101. Krechetnikov, Rouslan & Homsy, George M
J. Coll. Int. Sci. (2), 555–9.
Lee, M, Longoria, RG & Wilson, DE
Phys. Fluids (3), 540–550. Lhuissier, H. & Villermaux, E.
J. Fluid Mech. , 5–44.
Lister, John R., Kerr, Ross C., Russell, Nick J. & Crosby, Andrew
J. Fluid Mech. , 1–26.
Longuet-Higgins, M.S. & Oguz, H.
J. FluidMech. , 183–201.
Longuet-Higgins, M. S.
J. Fluid Mech. , 103–121.
Mandre, S., Mani, M. & Brenner, M. P.
Phys. Rev. Lett. , 134502.
May, A.
Int. J.Numer. Meth. Fluids , 1219–1222. Oguz, H. N. & Prosperetti, A.
J. Fluid Mech. , 111–145.
Richardson, E. G.
Proc. Phys. Soc. (4),352–367. Scolan, Y.-M. & Korobkin, A. A.
J. Fluid Mech. , 293–326.
Thoroddsen, S. T., Etoh, T. G. & Takehara, K.
J. Fluid Mech. , 125–134.
Villermaux, E. & Bossa, B.
J. Fluid Mech. , 412–435.
Wagner, Herbrt
Z.Angew. Math. Mech. (4), 193–215. Wilson, S. K.
J. Eng. Mathematics , 265–285. Worthington, A. M.
A study of splashes . London: Longman and Green.
Worthington, A. M. & Cole, R. S.
Philos. T. R. Soc. A , 137–148.
Worthington, A. M. & Cole, R. S.
Philos. T. R. Soc. A , 175–199.
Yakimov, Yu. L.
Fluid Dyn. (5), 679—-682. Yarin, A. L.
Ann.Rev. Fluid. Mech. , 159–192. Zeff, B.W., Kleber, B., Fineberg, J. & Lathrop, D.P.
Nature (6768), 401–404.
Zhang, Li V., Brunet, Philippe, Eggers, Jens & Deegan, Robert D.
Phys. Fluids22