Droplet impact on a thin liquid film: anatomy of the splash

Droplet impact on a thin liquid film: anatomy of the splash

Abstract

We investigate the dynamics of drop impact on a thin liquid film at short times in order to identify the mechanisms of splashing formation. Using numerical simulations and scaling analysis, we show that the splashing formation depends both on the inertial dynamics of the liquid and the cushioning of the gas. Two asymptotic regimes are identified, characterized by a new dimensionless number : when the gas cushioning is weak, the jet is formed after a sequence of bubbles are entrapped and the jet speed is mostly selected by the Reynolds number of the impact. On the other hand, when the air cushioning is important, the lubrication of the gas beneath the drop and the liquid film controls the dynamics, leading to a single bubble entrapment and a weaker jet velocity.

I Introduction

Droplet collision and impact are iconic multiphase flow problems: rain, atomization of liquid jets, ink-jet printing or stalagmite growth all involve impact in one manner or the other (Rein, 1993; Yarin, 2006; Josserand & Thoroddsen, 2016). The droplet may impact on a dry surface, a thin liquid film or a deep liquid bath. In all cases, impact may lead to the spreading of the droplet or to a splash where a myriad of smaller droplets are ejected far away from the zone of impact (Rioboo et al., 2001). Control of the outcome of impact is crucial for applications. Spreading is desirable for coating or ink-jet printing for instance while splashing may improve the efficiency of evaporation and mixing in combustion chambers. Two distinct types of splash, prompt and “ordinary” are now distinguished. The prompt splash is defined as a very early ejection of liquid at a time where is the droplet diameter and its velocity. The second, ordinary type occurs at larger times, through the formation of a vertical corolla ending in a circular rim that destabilizes into fingers and droplets (Deegan et al., 2008). In the prompt splash, a very thin liquid jet first forms, called the ejecta sheet (Thoroddsen, 2002). This sheet is ejected at high velocity, initially almost horizontally and is expected to disintegrate eventually in very small, fast droplets. Very often, these two mechanisms happen in a sequence, the ejecta sheet being a precursor of the corolla, as illustrated by the numerical simulations of droplet impacts on a thin liquid film shown in Figure 1.

a)b)
c)d)
e)f)

Figure 1: Snapshots of a droplet impacting on a thin liquid film for parameters in Tables 1 and 2, i.e. and . For the sake of visualisation, the liquid of the droplet and of the layer have been colored in green and blue respectively, although it is the same liquid. The gas phase is colored in dark blue. The snapshots correspond to the dimensionless time , ,, , and respectively.

The dynamics of droplet impact is complex, involving singular surface deformation and pressure values in the inviscid limit and several instabilities of surface evolution, so that an overall understanding of the whole process is still lacking. In particular, the splash depends on many physical parameters, the most important being the impact velocity. Obviously, high velocities promote the splash while at low velocities the droplet gently spreads. This behavior is mostly characterized by the Reynolds and the Weber numbers defined below. Although droplet impact on solid surfaces or on liquid films show similar output, the physical mechanisms leading to these effects often have different origins. For droplet impact on solids, the surface properties play an important role, through its roughness and the contact line dynamics for instance. A remarkable discovery has been done recently: the surrounding gas (usually air) also plays a crucial role in splash formation (Xu et al., 2005), and understanding in detail the influence of the gas still remains a challenge (Mandre & Brenner, 2012; Riboux & Gordillo, 2014; Klaseboer et al., 2014). In particular, the formation of a thin air layer at the instant of impact smoothes the singularity expected in the absence of any gas and thus “cushions” the impact. It also leads to the entrapment of an air bubble (Thoroddsen et al., 2003; Mehdi-Nejad et al., 2003; Thoroddsen et al., 2005; Mandre et al., 2009; Duchemin & Josserand, 2011),

In this paper we focus on droplet impact on a thin liquid film where splashing and spreading depend mostly on the balance, in the liquid as well as in the gas, between inertial and viscous forces (Stow & Hadfield, 1981; Mundo et al., 1995; Yarin & Weiss, 1995; Josserand & Zaleski, 2003). In the literature about impacts on liquid films, the effect of the surrounding air has not been shown as dramatically as on solid surfaces, although a systematic study of its influence is still lacking. For instance, it has often been noticed experimentally and numerically that air bubbles were entrapped by the impact dynamics (Thoroddsen et al., 2003; Thoraval et al., 2012), and the interplay between these bubbles and the ejecta sheet still needs to be elucidated. In 2003 two of us (Josserand & Zaleski, 2003, which we further refer to as JZ03) have proposed that the splashing/spreading transition observed in experiments and in numerics was controlled by the capillary-inertia balance within the ejecta sheet. The thickness in this theory was selected by a viscous boundary layer. In such a model, scaling laws for the jet thickness and velocity were deduced but the gas dynamics was totally neglected and the existence of an entrapped bubble was not considered at all.

The goal of this paper is therefore to determine the properties of the ejecta sheet using high resolution numerical simulations of axisymetric incompressible newtonian two-phase flow, in order to exhibit the relevant physical mechanisms at the heart of the prompt splash in this framework.

Ii The general problem

ii.1 Geometry and dimensional analysis

We consider a droplet of diameter impacting on a thin liquid film of thickness with a velocity normal to the film interface. Both liquid and gas have densities and , and dynamical viscosities and respectively. The surface tension is . In experiments the droplet is typically produced from some height above the film and falls under gravity . In our case the simulation starts with a small initial air gap between the film and the droplet and a velocity , as shown on Figure 2. The droplet is assumed to be spherical. In order for this assumption to be valid we need a) to have as little effect of the air flow on the droplet shape as possible. This should be verified if the gas Weber number is small or in a viscous regime and b) to assume that the oscillations of the droplet shape caused by the droplet release mechanism (Wang et al., 2012) are as small as possible. The Froude number that quantifies the influence of the gravity during the impact is taken constant and high () for all the simulations, indicating that gravity has only a small effect on the dynamics. We will restrict this study to large liquid Weber numbers .

Figure 2: Sketch of the numerical simulation, a spherical droplet of radius and velocity is considered at a distance from a liquid film of thickness .

The problem is then mostly characterized by the dimensionless numbers

which are the liquid Reynolds and Stokes numbers of the impact, and the density ratio. These numbers compare the droplet inertia with viscous effects in the liquid and gas respectively, and compare the liquid to the gas inertia. We note that has to be very large if one wants simultaneously to have and . An additional dimensionless number is the aspect ratio between the liquid film and the droplet which we keep relatively small and constant in this study. We expect the initial gas layer aspect ratio to be irrelevant if the conditions described above () are satisfied and is larger than the characteristic thickness defined below. Compressibility effects are characterized by the Mach numbers , where (resp. ) is the speed of sound in the liquid (resp. gas) and in all our simulations, these Mach numbers remained small enough so that compressibility effects could be neglected (Lesser & Field, 1983).

We note that the axisymmetric flow assumption is not valid when digitations and splash droplets form. However, at short time, before these instabilities can develop and particularly before the jet is created, we can consider this assumption correct citer et s’inspirer des resultats experimentaux. Otherwise the general 3D problem remains a grand numerical challenge because of the large range of scales involved (Rieber & Frohn, 1998; Gueyffier & Zaleski, 1998). Despite some recent numerical results (Fuster et al., 2009) realistic 3D numerical simulations of droplet impact at short times are yet hard to attain.

Moreover a solid basis for the analysis of the scaling of 3D flow may only be attained when the scaling of 2D flow has been uncovered. We thus postpone a detailed 3D study of droplet impact to future work.

ii.2 Scaling analysis

We analyze now the different mechanisms at play during droplet impact using simple scaling arguments. Recall that surface tension and gravity can be neglected, and in a first step, we will consider also that the surrounding gas has negligible effects. We then quite naturally define as the time at which an undeformed, spherical droplet at uniform velocity would touch the undeformed, planar liquid surface. With this definition of the time, the origin initial time is

We will now use an important geometrical argument first suggested by Wagner (1932): considering the intersection of the falling sphere with the impacted film, it is straightforward to define the vertical lengthscale as and the horizontal one . These apparently simple scalings arising purely from geometry are in fact very robust and relevant to the description of impact at short times: for instance, it has been shown that gives a correct estimate of the so-called spreading radius defined as the radius where the pressure is maximal. Remarkably, the geometrical velocity of this intersection

(1)

diverges at , questioning the incompressible assumption. However, although formally such a geometrical velocity diverges, fluid velocities remain much smaller and compressibility can be safely neglected for small Mach numbers. To make the article self-contained, we will now recall rapidly the results obtained by JZ03. The key point is the numerical observation that the pressure field and the velocity field are perturbed over the length scale so that a kind of inner-outer asymptotic analysis can be performed, in which the flow is uniform at scales larger that , potential at scales of order and viscous at scales much smaller than (a more rigorous asymptotic analysis has been developed later in Howison et al. (2005)). In this analysis the ejecta speed is obtained using a mass conservation argument between the impacting droplet and the ejected sheet, assuming that the thickness of the jet is selected by a viscous length. More precisely, one can compute in this framework first the mass flux from the falling undeformed sphere through the undeformed film surface

(2)

This flux can be absorbed either by surface deformation of the droplet and of the film or by the formation of an ejecta sheet. In JZ03, we have assumed that the thickness of such a sheet or jet is given by a viscous boundary layer formed at the basis of the jet leading to a viscous length scale

(3)

Then conservation of the volume flux through this jet implies that the jet velocity has to satisfy

(4)

Remarkably, this gives a nonlinear relationship between the jet and the impacting droplet velocity since then . Such a law has obviously some physical restrictions: first of all, the flux formula (2) is valid only for since our whole analysis is for short times, and makes no sense for times of order . Furthermore, the jet velocity has to be larger than the geometrical velocity . Indeed if one would expect the ejecta sheet to be overrun by the falling droplet. This condition together with (1) yields a “geometric” limiting time

(5)

When the ejecta forms, a bulge or rim appears at its tip according to the mechanism of Taylor (1959) and Culick (1960). This rim moves backwards at the Taylor-Culick velocity

(6)

where is the thickness of the ejecta sheet. The ejecta sheet cannot form if its velocity is smaller than the Taylor-Culick velocity constructed with the thickness of the ejecta sheet. We thus obtain that the ejecta can form only when which from equations (4, 3, 6) yields

(7)

Both conditions (5) and (7) must be satisfied at short times because as stated above the whole theory does not make sense for larger times. Then we must have which yields the condition

We restrict the present study to the dynamics where splashing is always present. More precisely since we consider situations such that and , then of the two conditions for splashing given above, is always more restrictive than . Therefore, in our configuration, the jet appears only when its velocity is bigger than the geometrical one.

Finally, the inertial pressure of the impact can be computed using the rate of change of vertical momentum in the droplet, following JZ03:

(8)

In this equation, the vertical momentum in the droplet is affected only in a half-sphere of radius . Eq. (8) gives the impact pressure

(9)

which leads to

(10)

as observed in numerical simulations (see JZ03). Note that a detailed analysis of the potential flow for a droplet falling on a solid surface has been performed by Philippi et al. (2015). The reasoning in the latter paper may be straightforwardly transposed to the impact on a liquid surface to obtain for the pressure field in the neighborhood of

(11)

which is very similar at to the scaling in (10). However the pressure field of (11) has an additional singularity for , not predicted by Equation (10) at . This singularity is indeed observed in our numerical simulations as well as in JZ03 and in Philippi et al. (2015); Duchemin & Josserand (2011).

In the theory above, contact occurs at , the vertical length scales are and , the flow pressure and the geometric velocity are singular with an infinite limit at . However, the effect of the gas layer was not taken into account so far. When instead the gas layer is taken into account, the above analysis is an approximation valid at length scales where is the thickness of the gas layer. The scale of the impact pressure is thus

(12)

The singularity of velocity and pressure is regularized and it can be said that the gas “cushions” the shock of the impact.

Moreover it is observed in experiments by Thoroddsen et al. (2003) and in numerical simulations (Josserand & Zaleski, 2003; Mehdi-Nejad et al., 2003; Korobkin et al., 2008) that contact does not occur on the symmetry axis but on a circle of radius so that a bubble is entrapped, as observed also for droplet impact on solid surfaces by Thoroddsen et al. (2005); Kolinski et al. (2012). As before, horizontal and vertical length scales are related at short times by . Thus the thickness of the gas layer or bubble at the time of contact relates to the contact radius by

(13)

These short-time asymptotics have to match the initial conditions at negative time . Let , be the positions of the drop and film surfaces respectively, and be the thickness of the gas layer. To fix ideas, let us consider initial conditions such that , close to the impact time so that on the axis . Then the gas layer is thin, there is a separation of horizontal and vertical length scales so that a lubrication approximation is valid over distances . As long as the lubrication pressure (estimated below) thus obtained in the gas layer is much smaller than the impact pressure the liquid advances almost undeformed while expelling the gas. In this regime

(14)

When the lubrication pressure becomes large enough to deform the liquid and slow down the thinning of the gas layer, the time is of order of a so called air-cushioning time scale and the thickness reaches the air cushion length scale . Matching with the initial velocity, eq. (14) then gives the relationship which together with the separation of scales condition (13) links the time and position of the contact through .

In order to determine these air cushioning space and time scales, we find the dominant balance in the lubrication equation, following in part recent works on impacts on solid surfaces (Korobkin et al., 2008; Mandre et al., 2009; Duchemin & Josserand, 2011; Klaseboer et al., 2014). Our theory starts from the incompressible lubrication equation in cylindrical geometry:

(15)

where is the pressure in the gas layer. The factor in front of the lubrication pressure comes from the Poiseuille velocity profile valid for laminar flows, obtained with a zero radial velocity at and , which is assumed because of the small horizontal velocity in the liquid before splashing. Using the above geometrical argument for the pressure term and , we obtain the following scaling for the lubrication pressure in the gas film of thickness :

The usual lubrication scaling for the bubble entrapment is then obtained by writing that during this “cushioning phase”, the air pressure for balances the impact pressure, that is yielding:

Here we have used the relation , where is the time of bubble entrapment. This relation gives the following scaling for the bubble entrapment:

(16)

The cushioning phase starts when is negative and of order and ends when first contact occurs at a positive time . This leads to two remarks: one is that the time of “cushioning” is both the time scale of the duration of this phase, and the time coordinate of the two instants of the beginning of the cushioning phase and the end of it at the first contact. We do not have strong arguments or data to show that this two instants are symmetric around . However, interestingly, our numerical simulations show that a kind of droplet/film symmetry holds, so that at , we have to a high degree of accuracy (see Figure 4 below).

In this approach, in agreement with previous works  (Mandre et al., 2009; Duchemin & Josserand, 2011), gas inertia effects have been totally neglected but such an assumption is questionable, as suggested in a recent study (Riboux & Gordillo, 2014). We propose here to take into account inertial pressure losses for the pressure in the gas layer using simple scaling arguments. Implementing inertial effects in thin film equation remains a difficult problem which is still a question of scientific debate (Luchini & Charru, 2010) and here we will only estimate scaling laws for the contribution of inertia.

The variation of the film height is given by the divergence of the horizontal gas flux in the layer

where is the averaged horizontal velocity between and . We can determine the scale for by considering the momentum balance in a thin gas layer

where is a constant depending on the profile of the flow in the gas layer and is a constant characterizing turbulent friction. This equation will hold if the flow remains thin () and does not separate. Since in incompressible flow pressure is defined up to a constant and the pressure at the exit of the thin gas layer flow (taken here for ) is the pressure at infinity, it is convenient to set this pressure at the exit to zero. Then equals the pressure difference and can be estimated at the bubble entrapment yielding

(17)

Here, we have taken for the first (dominant) term the lubrication pressure already computed above. The second and third terms result from two kinds of inertial effects, the turbulent friction term and a possible singular head loss due to flow separation. It is readily seen that the ratio between the first two terms is the local Reynolds number of the gas layer . The third term is a singular head loss. The constant and depend on the precise geometry of the flow and are difficult to estimate. However, it can be seen that the singular head loss is much smaller for our problem by a factor than the turbulent friction term so that we will neglect the singular head loss in the following developments.

The film pressure may be finally obtained by estimating the horizontal velocity scale as which together with eq. (13) yields

(18)

Remarkably, the neglected singular head loss term would give an additional contribution in the form .

Finally, equating and yields now an implicit equation for

(19)

which can be written in terms of the dimensionless variables , and

(20)

The above equation is cubic in and its solution gives the dimensionless height of the film as a function of and with two asymptotic regimes separated by a critical Stokes number

(21)

The first regime, for corresponds to the case computed above with lubrication only (16) and is of the form

as also stated by Mandre et al. (2009); Mani et al. (2010). In the other regime, for we get

(22)

The estimates for the time at which the bubble is entrapped results from the estimates for through . For we recover eq. (16)

(23)

where the dimensionless time is noted , and all the other scalings are obtained straightforwardly from the scalings for the height given above. Similarly from Equation (10) the pressure in and on top of the gas layer, neglecting surface tension effects is (eq. 16):

(24)

This theory will now serve as a framework to interpret the numerical simulations reported below. The main prediction is that air delays contact by a time of order and that a bubble of typical radius is entrapped.

The above considerations however do not say at what time the liquid sheet is ejected. It sets contraints on the air cushioning effect and we can only turn to numerical simulations to see how the air layer dynamics interacts with jet formation.

To conclude this scaling analysis, it is interesting to notice that an alternative theory has been proposed recently (Klaseboer et al., 2014). There a different scaling for the entrapment has been obtained (leading to instead of ) based on the balance between the lubrication pressure in the gas and the Bernoulli pression in the drop . Although the pressure amplitude in numerical simulations of drop impact has been shown to obey the singular law (11), experimental studies have yet to distinguish between these two predictions.

Iii Numerical method

In our continuum-mechanics modelling approach, fluid dynamics is Newtonian, incompressible, with constant surface tension. In the “one-fluid approach” (Tryggvason et al., 2011) one considers a single fluid with variable viscosity and density, and a singular surface tension force, yielding the Navier-Stokes equations that read:

(25)
(26)

where , is the pressure, denotes the unit normal to the interface and is the two-dimensional Dirac distribution restricted to the interface, and are the space-dependent fluid densities and viscosities equal to their respective values and in each phase. This set of equations can be written using dimensionless variables, rescaling lengths by , velocities by , times by , densities by and pressures by so that the Navier-Stokes equations become:

(27)

where now in the liquid phase while and in the gas. Equation (27) is solved using the methods described in (Popinet, 2003, 2009; Lagrée et al., 2011; Tryggvason et al., 2011), that is by discretizing the fields on an adaptive quadtree grid, using a projection method for the pressure, the time stepping and the incompressibility condition. The advection of the velocity fields is performed using the second-order Bell-Collela-Glaz scheme, and momentum diffusion is treated partially implicitly. The interface is tracked using a Volume of Fluid (VOF) method with a Mixed Youngs-Centered Scheme (Tryggvason et al., 2011) for the determination of the normal vector and a Lagrangian-Explicit scheme for VOF advection. Curvature is computed using the height-function method. Surface tension is computed from curvature by a well-balanced Continuous-Surface-Force method. Density and viscosity are computed from the VOF fraction by an arithmetic mean. This arithmetic mean is followed by three steps of iteration of an elementary filtering. This whole set of methods is programmed either in the Gerris flow solver (Popinet, 2014), or in the Gerris scripts that were designed to launch these computations.

Four refinement criteria are used as follows 1) the local value of the vorticity, 2) the presence of the interface as measured by the value of the gradient of the VOF “color function” 3) a measure of the error in the discretisation of the various fields based on an a posteriori error estimate of a given field as a cost function for adaptation. This a posteriori error is estimated by computing the norm of the Hessian matrix of the components of the velocity field, estimated using third-order-accurate discretisation operators, 4) when near the interface, the curvature is used as the adaptation criterion. To measure the degree of refinement so obtained, recall that on a quadtree grid, a level of refinement means that the grid cell is times smaller than the reference domain or “box”. When adaptively refining, a predefined maximum level is used for the adaptation on curvature (4), moreover adaptation on vorticity (1) and on the error (3) may lead to a maximum level of refinement and finally cells near the interface (2) are always refined to level at the least. Note that criterion (3) is generally more efficient than (1) so the latter could have been dropped altogether.

Iv Results of simulations

iv.1 Impact dynamics

We perform series of simulations with parameters set as in Table 1. For water-like fluids, these constant Weber and Froude numbers correspond to a drop of radius mm falling at velocity . The numerical simulations are performed for different liquid and gas viscosities characterized by the Reynolds and Stokes numbers varying from to and from to respectively. For a mm diameter drop impacting at , it would typically cover the range between one eighth to twenty times the water viscosity, and one fourth to ten times the air viscosity. In particular, we have done simulations for three Stokes numbers (, and ) a large range of Reynolds numbers (, , , , , , and ) that will be used to analyse the dynamical properties of the impact.

We Fr
500 815 826.4 0.2 1/30
Table 1: Dimensionless values of the parameters for all simulations reported

In an initial phase the droplet falls undeformed until air cushioning effects set in. Then at some point in time the jet forms and a reconnection of the interfaces on the droplet and film occurs. Figure 1 shows a droplet impact with physical parameters approximating a glycerinated water droplet falling in air. The main dimensional parameters are given in Table 2, other parameters approximate air at ambient temperature, leading to the dimensionless numbers and in complement to the dimensionless numbers of Table 1.

4 m s m kg ms kg s m
Table 2: Dimensional values of the main parameters for Figure 1

The grid is refined based on the four criteria above so that the smallest cell has size . Figure 3 shows two views of the grid refinement for a case where the liquid viscosity is twice higher than in Table 2 the gas one ten times larger, so that and , all the other parameters being the same than on figure 1

a)b)

Figure 3: The grid refinement used in a) Figure 5a and b) Figure 5d. Here, the viscosity is twice that of Table 2, so that with a gas viscosity such that .

We have checked that higher does not change the results significantly. This can also be verified from Figure 4 where it can be seen directly that the simulation starts at a time that is much larger than the apparent time scale of the air cushioning effect. In fact, the velocity of the south pole of the droplet in unperturbed until very short times around .

Figure 4: The positions of the bottom of the droplet and the top of the film on the axis as a function of time for , . The dashed line represents the mean position that decreases like .

In Figures 1 and 4, one can observe that a bubble is indeed entrapped by the impact due to the cushioning of the gas beneath the droplet. A very thin ejecta sheet is formed, followed by the growth of a thicker corolla. These figures do not however give a full account of the level of accuracy reached in the calculation, as shown on Figure 5 where successive zooms of the interface are shown around the instant when the droplet contacts the liquid film. The large range of scales between the droplet diameter and the small features in Figure 5d is apparent, Figure 3 showing the corresponding grid.

a) b)
c) d)

Figure 5: Successive zooms of the interface at the same time of impact for and . a) shows the general view; b) shows the zoom of the interface corresponding to the square traced on figure a) and so on from b) to c) and c) to d). In the first and the latter figure, the physical scales are shown based on a mm diameter drop. The scale shown in the latter figure corresponds in fact to twenty times the smallest mesh size of the numerical simulation.

In particular, small bubbles (which are actually toroidal because of the axial symmetry) can be seen prior to the ejection of the thin liquid sheet. In this case, the liquid of the droplet has already made contact with the liquid layer before the jet formation. Such small toroidal bubble entrapment might correspond to those observed in experiments recently (Thoroddsen et al., 2012). However it is also strongly controlled by the size of the grid, as reconnection in the VOF methods depends on the grid size. Here we note that with the grid size is one order of magnitude away from the molecular length scales, so the VOF reconnection although not physically realistic, may, in the future, approach the length scales at which molecular forces trigger reconnection in the real world. Finally, the mechanism of jet formation can be observed in Figure 6, where the vorticity field both in the liquid and gas phase is shown prior to the ejection corresponding to the zoom of Figure 5 d). It exhibits a vortex dipole at the origin of the jet, as already described in JZ03.

a)b)

Figure 6: The vorticity field corresponding to Figure 5d : a) in the liquid and b) in the gas.

We can now characterize the splashing dynamics as the liquid and gas viscosities vary.

iv.2 Spreading radius

One of the crucial quantities involved in the scaling analysis is the geometrical radius that acts as the horizontal length scale. In order to verify that the horizontal scale behaves like , we investigate the evolution with time of the spreading radius, defined as the point in the liquid where the velocity is maximal. Figure 7 shows the evolution of the spreading radius with time for all the simulations performed. The square-root scaling () is observed over a large range of time with the same prefactor for all simulations. Remarkably, the figure shows that this geometrical argument for the horizontal characteristic length is particularly robust and that the liquid properties (viscosities, densities) only influence the dynamics at short times.

Figure 7: The dimensionless spreading radius defined as the radius where the velocity is maximal in the liquid, as function of the dimensionless time , in a log-log plot for all the simulations performed in this study. The straight line indicates the slope corresponding to the geometrical law shifted below for eye guiding. As suggested by this geometrical relation, the different curves collapse all on a single one almost parallel to the expected law, showing that viscous, capillary and lubrication effects only alter slightly this dynamics. Differences can however be seen at short times.

iv.3 Initial gas sheet formation

We now study the formation of the gas sheet and its scaling. The transition from free fall to air cushioning can be seen on the time history of the heights of the film and the bubble and . Figure 4 shows both heights and on the axis as a function of time for and . This corresponds to a value of the air viscosity times its ordinary value, while the liquid is times more viscous than water. It is seen that the heights behave linearly until some time near (approximately . The linear behavior of before impact is an indication that the bottom of the drop falls at the free fall velocity (the gravity correction is even smaller than due to the short time of observation) almost unperturbed from its initial value. At on the other hand the cushioning dynamics have fully set in. After time 0, and remains approximately constant. This half velocity linear decrease of both and can be understood simply by momentum conservation as already suggested by Tran et al. (2013). The scale of the gas layer may thus conveniently be defined as .

To determine the scaling of two series of simulations have been performed at and for variable Stokes number . Together with the numbers in Table 1 these completely define the simulation parameters. The dimensionless height is plotted on Figure 8 together with relation (20). The unknown turbulent friction coefficient has been fitted by trial and error to .

Figure 8: Height scale of the gas layer as a function of the Stokes number for two values of the Reynolds number (red square) and (black circle). The fit using equation (20) is shown (solid line) as well as the two scaling laws (dotted line) and (dashed line).

While the numerical data points are not exactly on top of the fit the hypothesis of a transition from a limit at small to a behavior at larger (but still small ) is compatible with the data, with in some range around . However, it is worth to remark that the rightmost part of the graph is closer to a power law behavior, suggesting that the alternative scenario proposed by Klaseboer et al. (2014) might be valid here. In fact, as the analysis on the pressure will demonstrate it below, this alternative scaling is not valid here and the has to be seen as a best fit scaling in an intermediate regime (Jian et al. (2015)). In order to test the scaling of at very low , when the effect of is most marked, we perform a series of simulations at the smallest value of in Figure 8 and variable , keeping all other numbers constant. The results are plotted on Figure 9. We observe a linear increase of with , in agreement with the linear relation (22), with a constant in the limit that depends on the Stokes number.

Figure 9: Height scale of the gas layer as a function of for and and all other parameters as in Table 1.

iv.4 Impact pressure

In order to investigate quantitatively the various mechanisms involved in the impact dynamics and the jet foramtion, we follow the evolution in time of the maximum pressure on the axis in the gas layer as shown on Figure 10 for different Stokes number at constant . We see a very large pressure peak compared to the Bernoulli pressure . The insert of figure 10 shows a good agreement with the scaling predicted in (24). This dependence of the pressure field with the Stokes number is clearly in disagreement with the Bernoulli argument proposed in Klaseboer et al. (2014). Together with the variations of found compatible with a behavior, it indicates that the lubrication is the dominant regime in the air cushioning, by opposition to the alternative scenario of Klaseboer et al. (2014). We should also emphasize here that this high pressure in the gas layer can lead to the compression of the gas. Indeed, taking for instance the typical values for water drop impact and , we obtain a pressure un the air of the order of one half of the atmospheric pressure.

Figure 10: The pressure in the gas layer in the cushioning regime defined as the maximum gas pressure on the axis down the drop. The dimensionless pressure (using ) is plotted as a function of the dimensionless time for a fixed and for varying Stokes number , , , , , , and , from top to bottom. The insert shows the maximum pressure over time of these curves as function of the Stokes number. The dashed line draws the expected scaling following the prediction (24)

iv.5 Ejecta sheet velocity

We study now the evolution of the velocity maximum in the liquid. This quantity is indeed an interesting proxy for several measurements. It is directly related to estimates of the velocity of the jet itself. It is easier and less ambiguous to measure than the ejecta thickness, which varies widely at its base. Finally it has a telltale “spike” at the instant of jet ejection marking , the time of emergence of the ejecta sheet. On Figure 11 we show the maximum dimensionless velocity as a function of the dimensionless time shifted to the beginning of the simulation (remind that the origin of time corresponds to the time for which the geometrical falling sphere would interact with the liquid layer), for and , and other parameters in Table 1. It is seen that the maximum velocity deviates from the initial velocity around , that is when the droplet approaches the liquid film, after which it increases rapidly, then reaches a maximum and decreases slowly. The maximum is often remarkably spiked.

a

Figure 11: The maximum velocity as a function of the dimensionless time for and . Other parameters as in Table 1.
Figure 12: Left: the maximum velocity as a function of the radius where it is located for different times (indicated by the blue signs). Right: the interface profiles corresponding to the same time. The position at which the velocity is maximum is again marked by blue marks. The location of the maximum velocity at time of reversal of the interface curvature, marking the beginning of jet ejection, is marked with a red sign.

Figure 12 shows the value of the maximum velocity in the liquid and the location at which it is reached as the time varies for the same parameters. It is clear in that case that the sharp peak corresponding to the maximum velocity also corresponds to the time of reversal of the curvature of the interface, marking the beginning of the ejection of the jet. Detailed investigations show that in all cases investigated in this paper this spike corresponds exactly to the time of formation of the ejecta sheet, and that the maximum is located at its base.

However, zooming on the base of the jet as done in Figure 5d shows that a set of tiny bubbles is already formed, meaning that first contact between the drop and the sheet has already occurred before the jetting time in Figure 12. In other words in that case contact happens markedly before jet formation.

Defining the speed of the jet as this spike velocity, we can investigate how the jet velocity depends on the Reynolds numbers for the different Stokes numbers simulated. This velocity is shown on Figure 13 as function of the Reynolds number in order to check the validity of the scaling law (4): . The results are somehow puzzling: indeed, as the predicted law exhibits a reasonably good agreement for “low” Reynolds number (below ), important deviations appear at larger Reynolds where another scaling is apparently at play, consistent with a fit with . However, it is interesting to notice that the jet velocity shows almost no dependence on the Stokes number below , as suggested by the viscous length theory of JZ03, while a small dependance can be identified in the higher Reynolds regime.

Figure 13: The jet velocities as function of the Reynolds number for all the simulations performed in this study, in a log-log plot.The predicted power law relation 4 is plotted on showing only a reasonable agreement at Reynolds number lower than (solid line). For higher Reynolds numbers another scaling appears, consistent with a law (dashed line). For , where many Stokes numbers have been considered, remark a slight dependence of the velocity with the Stokes number, following the higher the Stokes number, the lower is the jet velocity.

In order to better understand this discrepancy between the predicted law and the numerical results, the dimensionless time of the jet formation needs to be investigated. Two scaling laws for this time are in competition: on the one hand the viscous length theory without accounting for the gas lubrication effect suggests that (relation 5); on the other hand, the cushioning of the gas suggests that this effect is delayed by the time (relation 23).

Figure 14: as a function of the inverse of the Reynolds number () for different Stokes number: black circle , red square and green diamond . Parallel dashed straight lines are drawn for each Stokes number for guiding the eyes.

First of all, Figure 14 exhibits parallel straight lines when plotting as a function of for three different Stokes numbers, indicating on the one hand that evolves linearly with . On the other hand, the straight lines are different for each Stokes number suggesting that a time delay depending on the Stokes number has to be considered. This can be seen on Figure 15 where is plotted, for , as a function of as suggested by the theoretical law obtained for the bubble entrapment  (23).

Figure 15: as a function of the Stokes number to the power predicted by the theory (23) , for constant Reynolds number . The dashed line indicates a linear relationship between these two quantities.

Again a nice straight line is observed, demonstrating eventually that the obeys the following relation:

(28)

which is equivalent by elementary algebra to

(29)

with and fitting parameters. Figure 16 confirms this relation, by plotting as a function of , together with the latter relation (29) for and , showing a good collapse of the data on this master curve.

Figure 16: The data (black circles) for the jet formation time plotted as the cloud of points for all the simulations performed in this study. Formula 29 is shown (solid curve) for and .

From Figure 14 one can argue that the fitting constants and may eventually depend slightly on the Reynolds and Stokes numbers respectively. However, no clear and strong trend could be extracted from the data available so far and we will consider and as constant as a first approximation. The above study suggests that the impact dynamics can be decomposed in two dynamical stages: a first one dominated by the cushioning dynamics involving a time scale dependence. Then, the ejection mechanism of the liquid sheet arises after a time delay proportional to .

a) b)

Figure 17: a) Maximum velocity as function of time for and for various Reynolds number (, , , , , , and for the curves from bottom to top). b) Same curves where the velocity has been rescaled by the theoretical prediction . In this case the order of the curves for increasing Reynolds numbers are reversed, ranging from top to bottom. The initial velocity has been subtracted at short time for clarity.

These different regimes can now be investigated through the evolution of the maximum velocity as function of time for all the parameters simulated here. Firstly, we show on figure 17 a), the maximum dimensionless velocity as function of the dimensionless time for different Reynolds number for a fixed Stokes number , where it can be observed that the higher the Reynolds number, the higher is the velocity as expected by the JZ03 prediction (4). This is investigated in figure 17 b), where these dimensionless velocities rescaled by the predicted scaling are plotted as function of time. If a reasonable collapse of the curve is obtained firstly at short times and for Reynolds numbers below , the curves for higher Reynolds numbers deviates from the master curve starting at the dimensionless time of the jet formation. As expected however, the time decreases as the Reynolds number increases, but for the high Reynolds numbers, the velocity peak appears during the velocity rise indicating that eventually the two mechanisms of air cushioning and jet formation interact. Somehow, the air cushioning effect is interrupted by the ejection of the liquid sheet. This could explain why the jet velocity at high Reynolds number does not follow the prediction (4). On the other hand, for lower Reynolds number, one can see that the two mechanisms of air cushioning and the jet formation are well separated in time.

a) b)

Figure 18: a) Maximum velocity as function of time for and for various Stokes numbers (, , , , and for the curves from left to right). b) For three Stokes numbers , and the maximum velocity rescaled by are shown for different Reynolds numbers ranging from to , in a log-linear plot. Three sets of curves can be identified, one for each Stokes number, from left to right as the Stokes number increases. Here, the initial velocity has been subtracted at short time for clarity.

The dependence of the dynamics on the Stokes numbers can be observed on figure 18 a) where the maximum velocity is shown at for different Stokes numbers. As expected by the lubrication theory, the higher the Stokes number the slower is the rise of the velocity. On the other hand, since the formation of the jet is delayed by this cushioning dynamics, we observe that the velocity peak is also delayed and is slightly decreasing as the Stokes number increases. This is in contradiction with the JZ03 initial prediction (4) that was obtained neglecting the gas cushioning, explaining why this relation is not verified. The velocity curves are rescaled on figure 18 b) by as suggested by the prediction (4) for three Stokes numbers and different Reynolds numbers (up to for a given Stokes number). The curves arrange in three sets (one for each Stokes number for all the Reynolds numbers) in time showing clearly that the Stokes number influences mostly the accelerating regime. The predicted scaling for the jet velocity is seen through the maximum of these curves that are very close one from each other. However, as observed in the jet velocity curve shown in figure 13 and on figures rescaleRe, many curves do not reach this maximum because of the cushioning dynamics.

a) b)

Figure 19: a) Rescaled velocities by for three Stokes number and up to eight Reynolds numbers shown as function of the rescaled time following the lubrication scaling 23. b) Same curves but with the rescaled time involving a scaling The initial velocity has been subtracted at short time for clarity.

The curves for the three Stokes numbers on Figure 18 result in three different curves which are roughly translated on the logarithmic time axis, suggesting a scaling dependence on the Stokes number as proposed by the relation (30). Using exactly this scaling to rescale the time as , figure 19 a) shows the maximum velocity for more than twenty different Stokes and Reynolds numbers. The collapse of all these curves is reasonably good, although a better collapse is obtained when rescaling the time by as shown on figure 19 b). This better collapse is somehow reminiscent of the best fit obtained for in figure 8. Therefore, the best collapse is obtained using a intermediate regime scaling where the time is best fit by . These results suggest that in this first dynamical stage, where the air cushioning is dominant, the velocity obeys the following relation:

(30)

valid a priori before the jet formation, where can be taken between and . is the universal function that describes this cushioning regime. Immediately after the jet formation, and in a kind of plateau region, the velocity at the base of the jet does scale with for intermediate Reynolds numbers as predicted by Equation (4), while it departs slightly from this scaling for high Reynolds numbers. Although this evolution is in agreement with the initial theory of JZ03, it is important to notice that the velocity is not constant and is in fact decreasing with time after the jet formation time . This second stage is determined by the jet dynamics and is a priori not influenced by the surrounding gas. Recalling that the timescale for the jet formation in the absence of gas cushioning is , it is tempting to investigate the evolution with time of the maximum velocity as a function of the rescaled time , following:

(31)

where is the universal function describing this second dynamical stage. The rescaled velocities are shown on figure 20 as function of the rescaled time for all the simulations performed in this study. Remarkably a nice collapse of the curves is observed for the large time dynamics i. e. for the time after the jet formation (). The dashed line in this log-log plot shows a good fit of the data in this regime, indicating a power law, so that the maximum velocity obeys eventually for :

(32)
Figure 20: The rescaled maximum velocity shown as a function of the rescaled time based on the geometrical timescale for jet ejection (see relation 5) in a log-log plot. The dotted-dashed line indicates the power law scaling . The initial velocity has been subtracted at short time for clarity.

This behavior suggests an explanation for the selection of the maximum jet velocity, considering that this latter regime starts at the jetting time . Recall firstly the asymptotic scalings observed for the jet maximum velocity shown on figure 13, namely for and otherwise. Interestingly, these two asymptotics are consistent with the function (32) when considering the two asymptotics for the jet formation time . Indeed, for we have so that, using 32, the first scaling is. On the other hand, in the other regime, , for which we obtain:

giving remarkably the velocity scaling for large Reynolds numbers. Note also that it predicts a power law dependence on the Stokes number that we have not tested so far, although it can be qualitatively seen on figure 13. Finally, using the formula (29) for :

we obtain an effective formula for the maximum jet velocity:

(33)
Figure 21: Black circles: the numerically measured jet velocity shown as function of the effective jet velocity eq. 33 for all the cases studied here. The best fit is obtained for while and have been computed before. The dashed line shows the theoretical line .

This formula is tested on figure 21 where the measured jet velocity for all the numerical simulations performed here is compared to eq. (33), where the best fit is obtained for . If the results are correctly distributed on the line for the jet velocities up to , it departs slightly from it for the higher velocities that correspond to the highest Reynolds number studied here.

V Discussion and Conclusion

In our numerical experiments, we have focused on the jet velocity and the effect of the gas layer. The formation and the scaling of the gas layer were analyzed. At the same time that the gas layer forms, a large peak of pressure is observed which scales as as predicted by the lubrication theory that we have outlined in Section II, leading eventually to a bubble entrapment. The formation of the gas layer itself shows interesting symmetry properties between the dynamics of the droplet surface and the dynamics of the liquid layer . The two are equally deformed at the time defined geometrically as described above. The time during which deformation sets up is very short compared to the characteristic time and also to the time of free fall . The corresponding layer thickness is consistent with a scaling as outlined by the lubrication theory, although an intermediate inertial regime is also observed. Moreover, Figure 10 demonstrates that the pressure is determined by the lubrication scaling, following , confirming that the dominant mechanism for the air cushioning is determined by the lubrication regime. The results for the jetting time exhibit an interplay between the gas cushioning (whose time scale is determined by ) and the liquid viscous boundary layer mechanism described in JZ03 (with a time scale in ). However, it is interesting to remark that the air cushioning dynamics could almost equivalently be fit by a time scale (instead of ), indicating that the intermediate inertial regime is also observed. The dependence of the layer thickness shows in addition a linear relationship with the density ratio (Figure 9) as predicted by the theory in equation (22).

These results evidence a transition from a regime where the gas thickness and cushioning effect are insignificant on the jet dynamics to a regime where the air cushioning controls the jet dynamics. From equation (29), we find that the air cushioning regime occurs for provided is small enough. The corresponding number varies only with for a given liquid-gas pair, so it fixes a limiting for which the regime changes for a given pair of fluids. For an air-water system, taking it gives a transition at or for a mmm droplet and . The jet regime is thus characterized by a new dimensionless number that is related to the ratio between the time for bubble entrapment to the geometrical time for the jet formation :

(34)

We have observed in our simulations that in most cases the jet emergence is simultaneous with the connection of the two interfaces, that is the bubble becomes trapped at the time of jet formation. In fact, the situation is more complex and can be analyzed using this jet number that quantifies the transition between a regime where the air cushioning is insignificant () to a regime when it is dominant (). When the air cushioning is insignificant, we find that the jet forms at the geometrical time . with a jet velocity of the order of as predicted earlier in JZ03. Moreover, we find that this velocity scale is present before jet formation indicating the existence of a large velocity in the droplet prior to the formation of the jet. This large velocity (asymptotically infinitely larger than ) is indicative of the focusing of the liquid velocity in a small region inside the droplet prior to the emergence of the jet. In the regime where the gas thickness and cushioning effect are insignificant at small the air layer has to close before jet formation. A trace of that is seen in the presence of small bubbles on Figure 5d before jet formation. This is even clearer on figure 22 that shows the details of the interface dynamics between the time of the first contact between the drop and the liquid film and the formation of the jet, for and , leading to .

a)b)
c)d)

Figure 22: Zoom on the connexion of the upper interface of the drop with the liquid layer for and , the other parameters being that of table 1. The white lines represent the interfaces and the colors scale in the liquid phase is realted to the velocity.The interfaces are shown on time a), b), c) and d). The jet number and we thus observe that after the first connexion of the two interface on figure a), small bubbles are created at the front of the connexion by the retraction of the thin gas layer that is still between the falling drop and the liquid film, figure b) and c) until the jet emerges on figure d). On the other side of the connected region, the thin air film forming the entrapped large bubble is also retracting by capillarity.

We observe that a delay exists between the first connexion of the interfaces (figure 22 a) and the jet formation (igure 22 d), during which small bubble are entrapped by the thin gas film dynamics.

On the other hand, in the regime where is large (and similarly it is not clear how jet formation and air bubble closing interact but one expects that no small bubbles are entrapped and that the jet is formed simultaneously to the bubble entrapment. This is illustrated on figure 23 for and , so that .

a)b)
c)d)

Figure 23: Zoom on the connexion of the upper interface of the drop with the liquid layer for and , the other parameters being that of table 1. The interface is shown on time a), b), c) and d). The jet number and the first connexion of the two interface occurs between figures b) and c). No additional bubbles are entrapped and the jet is directly emerging from this connexion, figure d).

The splashing dynamics appears thus as a complex combination between the inertial dynamics of the liquid with the cushioning of the gas.

References

  1. Culick, F. E. C. 1960 Comments on a ruptured soap film. J. Appl. Phys. 31, 1128–1129.
  2. Deegan, R., Brunet, P. & Eggers, J. 2008 Complexities of splashing. Nonlinearity 21, C1.
  3. Duchemin, L. & Josserand, C. 2011 Curvature singularity and film-skating during drop impact. Phys. Fluids 23, 091701.
  4. Fuster, D., Agbaglah, G., Josserand, C., Popinet, S. & Zaleski, S. 2009 Numerical simulation of droplets, bubbles and waves: state of the art. Fluid Dyn. Res. 41, 065001.
  5. Gueyffier, D. & Zaleski, S. 1998 Formation de digitations lors de l’impact d’une goutte sur un film liquide. C. R. Acad. Sci. IIb 326, 839–844.
  6. Howison, S., Ockendon, J., Oliver, J., Purvis, R. & Smith, F. 2005 Droplet impact on a thin fluid layer. J. Fluid. Mech. 542, 1–23.
  7. Jian, Z., Josserand, C., Ray, P., Duchemin, L., Popinet, S. & Zaleski, S. 2015 Modelling the thickness of the air layer in droplet impact. ICLASS 2015 .
  8. Josserand, C. & Thoroddsen, S. 2016 Drop impact on a solid surface. Annu. Rev. Fluid Mech. 48, 365–391.
  9. Josserand, C. & Zaleski, S. 2003 Droplet splashing on a thin liquid film. Phys. Fluids 15, 1650.
  10. Klaseboer, E., Manica, R. & Chan, D. Y. 2014 Universal behavior of the initial stage of drop impact. Physical review letters 113 (19), 194501.
  11. Kolinski, J. M., Rubinstein, S. M., Mandre, S., Brenner, M. P., Weitz, D. A. & Mahadevan, L. 2012 Skating on a film of air: drops impacting on a surface. Phys. Rev. Lett. 108, 074503.
  12. Korobkin, A., Ellis, A. & Smith, F. 2008 Trapping of air in impact between a body and shallow water. J. Fluid Mech. 611, 365–394.
  13. Lagrée, P. Y., Staron, L. & Popinet, S. 2011 The granular column collapse as a continuum: validity of a two-dimensional Navier–Stokes model with a (I)-rheology. Journal of Fluid Mechanics 686, 378–408.
  14. Lesser, M. & Field, J. 1983 The impact of compressible liquids. Annu. Rev. Fluid Mech. 15, 97.
  15. Luchini, P. & Charru, F. 2010 Consistent section-averaged equations of quasi-one-dimensional laminar flow. J. Fluid. Mech. 565, 337–341.
  16. Mandre, S. & Brenner, M. 2012 The mechanism of a splash on a dry solid surface. J. Fluid Mech. 690, 148–172.
  17. Mandre, S., Mani, M. & Brenner, M. 2009 Precursors to splashing of liquid droplets on a solid surface. Phys. Rev. Lett. 102, 134502.
  18. Mani, M., Mandre, S. & Brenner, M. 2010 Events before droplet splashing on a solid surface. J. Fluid Mech. 647, 163–185.
  19. Mehdi-Nejad, V., Mostaghimi, J. & Chandra, S. 2003 Air bubble entrapment under an impacting droplet. Phys. Fluids 15 (1), 173–183.
  20. Mundo, C., Sommerfeld, M. & Tropea, C. 1995 Droplet-wall collisions: Experimental studies of the deformation and breakup process. Int. J. Multiphase Flow 21, 151.
  21. Philippi, J., Lagrée, P.-Y. & Antkowiak, A. 2015 Drop impact on solid surface: short time self-similarity. In preparation.
  22. Popinet, S. 2003 Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries. J. Comp. Phys. 190 (2), 572–600.
  23. Popinet, S. 2009 An accurate adaptive solver for surface-tension-driven interfacial flows. J. Comput. Phys. 228, 5838–5866.
  24. Popinet, S. 2014 Gerris flow solver-http://gfs.sourceforge.net/ .
  25. Rein, M. 1993 Phenomena of liquid drop impact on solid and liquid surfaces. Fluid Dyn. Res. 12, 61.
  26. Riboux, G. & Gordillo, J. 2014 Experiments of drops impacting a smooth solid surface: A model of the critical impact speed for drop splashing. Phys. Rev. Lett. 113, 024507.
  27. Rieber, M. & Frohn, A. 1998 Numerical simulation of splashing drops. Academic Press, proceedings of ILASS98, Manchester.
  28. Rioboo, R., Marengo, M. & Tropea, C. 2001 Outcomes from a drop impact on solid surfaces. Atomization and Sprays 11, 155–165.
  29. Stow, C. & Hadfield, M. 1981 An experimental investigation of fluid flow resulting from the impact of a water drop with an unyielding dry surface. Proc. R. Soc. London, Ser. A 373, 419.
  30. Taylor, G. I. 1959 The dynamics of thin sheets of fluid III. Disintegration of fluid sheets. Proc. Roy. Soc. London A 253, 313–321.
  31. Thoraval, M.-J., Takehara, K., Etoh, T., Popinet, S., Ray, P., Josserand, C., Zaleski, S. & Thoroddsen, S. 2012 von kármán vortex street within an impacting drop. Phys. Rev. Lett. 108, 264506.
  32. Thoroddsen, S. 2002 The ejecta sheet generated by the impact of a drop. J. Fluid Mech. 451, 373.
  33. Thoroddsen, S., Thoraval, M.-J., Takehara, K. & Etoh, T. 2012 Micro-bubble morphologies following drop impacts onto a pool surface. J. Fluid Mech. 708, 469–479.
  34. Thoroddsen, S. T., Etoh, T. G. & Takehara, K. 2003 Air entrapment under an impacting drop. Journal of Fluid Mechanics 478, 125–134.
  35. Thoroddsen, S. T., Etoh, T. G., Takehara, K., Ootsuka, N. & Hatsuki, A. 2005 The air bubble entrapped under a drop impacting on a solid surface. Journal of Fluid Mechanics 545, 203–212.
  36. Tran, T., de Maleprade, H., Sun, C. & Lohse, D. 2013 Air entrainment during impact of droplets on liquid surfaces. Journal of Fluid Mechanics 726.
  37. Tryggvason, G., Scardovelli, R. & Zaleski, S. 2011 Direct Numerical Simulations of Gas-Liquid Multiphase Flows. Cambridge University Press.
  38. Wagner, H. 1932 über stoss und gleitvorgänge und der oberfläshe von flüssigkeiten. Zeit. Angewandte Math. Mech. 12 (4), 193–215.
  39. Wang, A.-B., Kurniawan, T. & Tsai, P.-H. 2012 About oscillation effect on the penetration depth of drop induced vortex ring. In XXIII ICTAM, 19-24 August 2012, Beijing, China.
  40. Xu, L., Zhang, W. & Nagel, S. 2005 Drop splashing on a dry smooth surface. Phys. Rev. Lett. 94, 184505.
  41. Yarin, A. & Weiss, D. 1995 Impact of drops on solid surfaces: self-similar capillary waves, and splashing as a new type of kinematic discontinuity. J. Fluid Mech. 283, 141–173.
  42. Yarin, A. L. 2006 Drop impact dynamics: Splashing, spreading, receding, bouncing… Annu. Rev. Fluid Mech. 38, 159.