Maximum drop radius and critical Weber number for splashing in the dynamical Leidenfrost regime

Depósito de investigación de la Universidad de Sevilla https://idus.us.es/ "This document is the Accepted Manuscript version of a Published Work that appeared in final form in Journal of Fluid Mechanics copyright © Cambridge University Press after peer review and technical editing by the publisher. To access the final edited and published work see 10.1017/jfm.2016.496. Under consideration for publication in J. Fluid Mech. 1 Maximum drop radius and critical Weber number for splashing in the dynamical Leidenfrost regime Guillaume Riboux†& Jos´ e Manuel Gordillo ´ Area de Mec´anica de Fluidos, Departamento de Ingenier´ıa Aeroespacial y Mec´anica de Fluidos, Universidad de Sevilla, Avenida de los Descubrimientos s/n 41092, Sevilla, Spain. (Received ?? and in revised form ??) At room temperature, when a drop impacts against a smooth solid surface at a velocity above the so called critical velocity for splashing, the drop loses its integrity and fragments into tiny droplets violently ejected radially outwards. Below this critical velocity, the drop simply spreads over the substrate. Splashing is also reported to occur for solid substrate temperatures above the Leidenfrost temperature, TL, for which a vapor layer prevents the drop from touching the solid. In this case, the splashing morphology differs from the one reported at room temperature because, thanks to the presence of the gas layer, the shear stresses acting on the liquid can be neglected. Our purpose here is to predict, for wall temperatures above TL, the critical Weber number for splashing as well as the maximum spreading radius. First, making use of Boundary Integral Simulations, we calculate both the time evolution of the liquid velocity as well as the height of the sheet which is ejected tangentially to the substrate. These results are then used as boundary conditions for the one dimensional mass and momentum equations describing the dynamics of the rim limiting the expanding liquid sheet. Our predictions for both the maximum spreading radius and for the critical Weber number for splashing, are in good agreement with experimental observations. 1. Introduction The understanding of the spreading or the break up processes of a drop impacting onto a solid substrate is an area of current active research because of its relevance in a number of technological applications such as coating, cleaning, cooling, and combustion (Josserand & Thoroddsen 2016). It is known that, at room temperature, drop splashing does not only depend on the radius R, on the impacting velocity Vand on the density ρ, the viscosity µand the interfacial tension σof the liquid, but also on the material properties of the gas, on the gas pressure (Xu et al. 2005) and on the physicochemical properties of the substrate (Duez et al. 2007). For the case of drops impacting a smooth and dry solid substrate at room temperature illustrated in figure 1, it can be observed that once the drop touches the solid at T= 0, a small bubble -which has no influence on the subsequent dynamicsis entrapped near the axis of symmetry. Subsequently, at the ejection time Te, a thin sheet of liquid is expelled from the radial position A(Te), with π A2(Te) the area of the wetted region. Then, if the impact velocity is below the critical velocity for splashing, V < V ∗, the ejected lamella spreads tangentially along the solid; however, if V > V ∗, the edge of the liquid sheet dewets the solid as a consequence of the lift forces exerted by the surrounding gas (Riboux & Gordillo (2014), from now on RG14). Then, capillary and Rayleigh–Taylor disturbances grow in the azimuthal direction of the †Email address for correspondence: grib[email protected] 2Guillaume Riboux & Jos´e Manuel Gordillo rim, causing its disintegration into drops. At room temperature, splashing will occur only if the liquid front is able to dewet the solid. However, when the temperature of the substrate is above TL, with TLthe Leidenfrost temperature, the liquid never touches the substrate, with independence of the value of the impact velocity because, under these conditions, the drop levitates on its own vapor (Shirota et al. 2016). Consequently, the critical velocity for splashing will strongly depend on whether the solid temperature is below or above TL, as it has been recently reported by Staat et al. (2015). In addition, as a consequence of the frictionless motion of the liquid with the substrate, the maximum spreading diameter of impacting drops increases for temperatures above TLwith respect to the case of drop impact at room temperature (Lastakowski et al. 2014; Wilderman et al. 2016). In this contribution we aim to predict, for temperatures of the substrate above TL, the maximum spreading radius of the drop as well as the critical impact velocity above which the drop disintegrates into smaller droplets. Moreover, for impact velocities larger than the critical one, we will also provide results for the velocities and for the diameters of the tiny droplets ejected. For this purpose, the flow field within the ejected lamella is modeled using a ballistic approximation and next, making use of integral balances of mass and momentum, the dynamics of the rim is described. A similar study was carried out by Riboux & Gordillo (2015) (from now on RG15), where the analytical expressions for the liquid velocity and the height of the liquid sheet deduced in RG14 were used as boundary conditions for the equations describing the flow in the expanding lamella. One of our main contributions here is that these analytical expressions are substituted by universal time-dependent functions calculated numerically using a Boundary Integral Method (BIM). We find that, for short times after impact, the simulations are in excellent agreement with the theoretical predictions in RG14 but, for times T∼R/V , these numerical results depart from those in RG14 and also from the equations provided in Roisman (2009); Eggers et al. (2010). The differences found are essential to correctly predict the experiments. Potential flow numerical simulations will be carried out using the Boundary Integral method described in Rodr´ıguez–Rodr´ıguez et al. (2006); Gordillo & Gekle (2010). This type of code is adequate to simulate the splashing of drops for substrate temperatures above TLsince: i) the vorticity in the falling drop is initially zero and ii) the shear stresses acting in the circular region of radius A(T) below the drop -see figure 1are negligible because, in the dynamic Leidenfrost regime, a thin vapor layer prevents the liquid from touching the wall. Therefore, the production of vorticity is restricted to regions near the free surface where the interfacial curvature is highest, e.g., the region where the impacting drop meets the ejected liquid sheet. Consequently, the flow field after the impact is mostly irrotational except in very narrow regions localized near the free surface. Motivated by this fact, in this contribution, we approximate the velocity field within the drop as v=∇Φ, with Φ the velocity potential satisfying the impenetrability condition at Z= 0, ∂Φ/∂Z = 0. In addition, thanks to the symmetry of the flow with respect to the plane Z= 0, vis calculated here as the result of the head on collision of two drops of identical radii Rmoving with respective velocities Vezand −Vez. Here, ez indicates the unit vector pointing in the opposite direction to that of the falling drop. The comparison of the experimental images corresponding to an ethanol drop impacting a solid substrate at room temperature with the numerical profiles is provided in figure 1. Except in the region near the edge of the lamella, the agreement between experiments and numerical results, is remarkable. The discrepancies in the position of the rim observed in figure 1 are due to the fact that, at room temperature, the viscous shear stresses at the wall contribute to further decelerate the advancing front. We used our own experiments, Maximum drop radius and splashing criterium in the dynamical Leidenfrost regime 3 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 ( a ) ( b ) ( c ) ( d ) ( e ) A ( T ) ( f ) ( g ) ( h ) ( i ) ( a ) ( b ) ( c ) ( d ) ( e ) ( f ) ( g ) Figure 1. Top: The sequence of images illustrate an ethanol drop of initial radius R=1.03 mm impacting against a dry smooth solid surface at room temperature with a velocity V=1.69 m·s−1 -We = ρV 2R/σ = 100 and Re = ρV R/µ = 1360at times (a)t≃0, (b)t= 0.10, (c)t= 0.19, (d) t= 0.29, (e)t= 0.38, (f)t= 0.48, (g)t= 0.57, (h)t= 0.76, (i)t= 0.95 with t=TV/R. The bright line represents the result of the Boundary Integral Simulation, which is carried out using a number of nodes Nthat increases dynamically in time (N(T= 0) = 501). Nodes are clustered in the regions with the highest curvature. The radius of the wetted area A(T)/R =a(t) = √3t, (Riboux & Gordillo 2014) is indicated in (e). Bottom: The right column is a zoom of the region where the lamella is ejected for different instants of time, namely, (a)t= 0.016, (b) 0.049, (c) 0.082, (d) 0.164, (e) 0.246, (f) 0.328 and (g) 0.491. For reference, the left column illustrates the numerical results for the same values of t. The absence of viscous shear stresses causes the rim to flow faster in the numerical simulations than in the experiments. performed at room temperature, to compare with the numerical results since we could not find in the literature experimental images of drops impacting a hot substrate in the Leidenfrost regime with enough spatio-temporal resolution. 4Guillaume Riboux & Jos´e Manuel Gordillo Figure 2. The potential flow numerical results start imposing the initial wetted radius, a= sin θ0≈π/36 with θ0= 5◦. The inset (a) shows the re-gridding procedure used to stabilize the numerical simulation, consisting in placing nodes in between the nodes of the previous time step. The inset (b) is a zoom of the region from which the lamella is ejected at three different instants of time: te,t= 0.15 and t= 0.35. The figure also illustrates the wetted radius a(t) as well as the radial position of the rim, rt(t), and its velocity, vt(t). The inset also shows the definition of height of the lamella, h(r, t), not to be confused with the width of the rim, ht(t). 0 0.03 0.06 0.09 0.12 0.15 0.18 t 0 0.2 0.4 0.6 0.8 1 te ( a ) rt ( t ) a ( t )  3 t ˙ rt ( t )= vt ( te ) 101102103 We 10-2 10-1 100 te ∝ We − 2 / 3 ( b ) Figure 3. (a) The figure shows the numerical results for both a(t) and the radial position of the edge lamella rt(t) in the case of We = 100. The vertical line indicates the ejection time predicted by equation (2.2) and the dashed line represents the predicted position of the rim, ejected from r=√3tewith an initial velocity vt(te) = 1/2p3/te, with tegiven in equation (2.2). (b) The ejection times predicted by the Boundary Integral Simulations are such that te∝We−2/3, in agreement with equation (2.2). 2. Potential flow simulations Following the notation in RG14, lower case variables will be used in what follows to refer to dimensionless variables, which are constructed here using as scales for velocity, length and pressure V,Rand ρV 2, respectively. Since the Froude number Fr = V2/gR 1 and the shear stresses acting on the liquid are negligible, the only relevant dimensionless parameter characterizing the splashing of droplets in the Leidenfrost regime is the Weber number, defined here as We = ρV 2R/σ. The most relevant event taking place after the impact of a drop in the dynamic Leiden- Maximum drop radius and splashing criterium in the dynamical Leidenfrost regime 5 frost regime is that at te, with tethe ejection time -see figure 2-, an extremely thin liquid sheet of initial thickness ht(te) is expelled tangentially to the substrate with an initial velocity vt(te). In RG14 an algebraic equation for tewas deduced based on the following facts: i) prior to the ejection of the lamella, the time evolution of the radius of the wetted area is a(t) = √3t†, ii) at te, the velocity of the tip of the lamella is equal to the velocity of the wetted radius, i.e., vt(te) = ˙a(te) with dots denoting time derivatives and iii) the lamella can only be ejected if its tip advances faster than the radius of the wetted area. We concluded in RG14 that tecan be calculated solving the algebraic equation √3 2Re−1t−1/2 e+ We−1= ¨a h2 t=c2t3/2 ewith c= 1.1,(2.1) which expresses the fact that the ejection time is the instant at which the deceleration of the edge of the lamella coincides with the deceleration of the wetted area, ¨a. For wall temperatures above TL, the edge of the lamella is not decelerated by the action of viscous shear stresses and, therefore, setting to infinity the Reynolds number Re in equation (2.1) yields, te= (c2We)−2/3.(2.2) Using the analytical expressions for both vt(te) and ht(te) derived in RG14 as well as the result in equation (2.2), it can be concluded that vt(te) = 1/2p3/te∝We1/3and ht(te) = √12 t3/2 e/π ∝We−1. Figure 3 confirms the result in equation (2.2), te∝We−2/3. Once the lamella is ejected, in RG14 we also deduced that, for t > teand t1, both the liquid velocity and the height of the liquid layer at r=a(t) = √3t, which is the radial position where the drop meets the lamella, are respectively given by va= 2˙a=p3/t and ha=√12 t3/2/(3π). Figure 4 confirm these theoretical predictions for We ⩾100. More precisely, figure 4 shows that, while the analytical expression for vais valid for arbitrary times, the corresponding analytical expression for hadeparts from the numerical results for t⩾0.1 (see figure 4c). Therefore, the equations describing the liquid velocity in the lamella u(r, t) and its height h(r, t) at r=√3ti.e., at the radial position from which the liquid sheet is ejected for times t > te, are respectively given by u(r=a, t) = va=p3/t , and (h(r=a, t) = ha=√12 t3/2/(3π) for t < 0.1 h(r=a, t) = P(t) for t⩾0.1, (2.3) with P(t) = 5 X i=0 (0.1pi)tiand p0=−2.453 ×10−3, p1= 1.321, p2=−1.176, p3= 0.4943, p4=−0.1047, p5= 8.89 ×10−3 (2.4) a polynomial which is fitted to the numerical results. The ad hoc radial velocity field within the drop proposed by both Roisman (2009) and Eggers et al. (2010) is u(r, t) = r/(t+τ), with τan adjustable order unity constant. While the stagnation point type of flow within the impacting drop hypothesized by Roisman (2009) and Eggers et al. (2010) is a good approximation to the real flow field for t∼O(1), for times such that t > te,t1, the velocity field at r≃√3tdoes not correspond to a †This result was derived for the very first time by Riboux and Gordillo in RG14 using the linearization of the boundary conditions for the potential flow (Wagner 1932). 6Guillaume Riboux & Jos´e Manuel Gordillo 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 t 0 1 2 3 4 5 6 7 8 9 va ( a ) We=100 ∗ We=100 We=300  3 /t 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 t 0 1 2 3 4 5 6 7 8 9 ha ( b ) × 10 − 2 (  12 / 3 π ) t 3 / 2 P ( t ) [1] [2] 0 0.05 0.1 0.15 0.2 t 0 4 8 12 ha × 10 − 3 ( c ) We=100 ∗ (  12 / 3 π ) t 3 / 2 P ( t ) Figure 4. (a) Comparison between the theoretical expression va= 2˙a=p3/t deduced in RG14 and the numerical result. The good agreement between theory and experiments depicted in this figure is independent of the Weber number whenever We ⩾100. Two of the simulations shown are carried out using N(T= 0) = 501 whereas in the case marked with an asterisk [∗], N(T= 0) = 1001. (b) Comparison of the height of the lamella at r=√3tpredicted by the models in Eggers et al. (2010) ([1]) and Roisman (2009) ([2]) with the numerical result. (c) Comparison between the theoretical expression ha=h(a(t), t) = (√12/3π)t3/2deduced in RG14, valid for t1 and the numerical result. The vertical line indicates the instant at which we have set the transition between the result predicted by the potential flow theory and the fitting polynomial P(t) in equation (2.4). stagnation point flow (see RG14 for details). Interestingly, our theory predicts a radial velocity at r=a(t) = √3t,u(r=a(t), t) = a(t)/t =p3/t for t1 which is, by coincidence, the radial velocity corresponding to a stagnation point of flow of the type u=r/t particularized at r=√3t. This is why figure 4ashows a fairly good between the analytical expression va=p3/t and the numerical results for arbitrary values of t. Figure 4balso shows the values of ha(t) = h(r=√3t, t) predicted by the model in Roisman (2009), h(r, t) = 8 η (t+τ)2exp −6η r2 (t+τ)2,(2.5) with η= 0.39, τ= 2 ×0.25 and with a constant equal to 8 because here we define dimensionless lengths using the drop radius instead of its diameter, as well as the results predicted by the model in Eggers et al. (2010), h(r, t) = 1 (t+τ)2H(x) H(x) = 3.19 (1 + C x2)6with x(r, t) = r t+τ        (2.6) with C= 0.604 and τ= 1 (Lastakowski et al. 2014). Clearly, neither the model by Roisman (2009) not that by Eggers et al. (2010) is in agreement with the numerical results. Maximum drop radius and splashing criterium in the dynamical Leidenfrost regime 7 0.1 0.2 z ( a ) 0.1 0.2 z ( b ) 0.1 0.2 z ( c ) 0.1 0.2 z ( d ) a ( t ) 0.1 0.2 z ( e ) 0.1 0.2 z ( f ) ha ( t ) 0.1 0.2 z ( g ) 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 r 0.1 0.2 z ( h ) Figure 5. Comparison, for We = 100, between the height of the lamella predicted by the model -ballistic equations (3.1) subjected to the boundary conditions given in equation (2.3)- and the numerical result. The vertical line illustrates the radial position r=a(t) = √3tfor several instants of time as well as the height of the lamella at r=√3t,ha(t) = h(r=√3t, t); vais the liquid velocity at r=a(t) = √3a, namely va=u(r=√3t, t). (a)t= 0.19, (b) 0.25, (c) 0.32, (d) 0.45, (e) 0.55, (f) 0.70, (g) 0.90 and (h)t= 1.10. 3. Prediction of the critical Weber number for splashing and of the maximum spreading radius for droplets impacting in the Leidenfrost regime For a given instant of time t > te, the lamella extends from the position where the drop meets the liquid sheet i.e., r=√3t, down to r=rt(t), with rt(t) indicating the radial position of the rim. In RG15, we showed that both the liquid velocity u(r, t) and the height of the liquid layer within the lamella, h(r, t), can be calculated using the pair of equations (Gordillo & Gekle 2010; Villermaux & Bossa 2011), Du Dt= 0 and D ln (rh) Dt=−∂ u ∂ r ,(3.1) with D/Dt≡∂/∂ t +u ∂/∂r the material derivative. The former of the two equations in (3.1) expresses that fluid particles conserve their velocities within the lamella from r=√3tdown to the radial position rt(t) where the rim is located; the latter, is the continuity equation. The pair of equations in (3.1) representing the ballistic motion of fluid particles along the lamella are solved subjected to the initial conditions given in equation (2.3) by means of the Lagrangian Numerical Method already used in RG15. The results of the model for h(r, t), shown in figure 5, are in excellent agreement with the numerical simulations. The radial position of the edge of the lamella, rt(t), as well as its thickness, ht(t), are deduced applying the integral balances of mass and momentum at the rim (Taylor 1959; 8Guillaume Riboux & Jos´e Manuel Gordillo Culick 1960), drt dt=vt, π 4 dh2 t dt= [u(rt, t)−vt]h(rt, t), π h2 t 4 dvt dt= [u(rt, t)−vt]2h(rt, t)−2 We−1.                (3.2) The system of equations (3.2) is solved subjected to the following initial conditions, rt=a(te) = √3te, vt=vt(te) = 1/2p3/te, ht=ht(te) = √12 t3/2 e/π ,        (3.3) with tegiven by equation (2.2). Figure 6ashow the solution of the system of equations (3.1)–(3.3) for a particular value of the Weber number, whereas figure 6billustrates the comparison between the predicted maximum spreading radius and the experimental data in Tran et al. (2012). As it was pointed out in Wilderman et al. (2016), there are two clearly differentiated regions in the plot rmax vs We of figure 6b: the impact is elastic for We .10 since the initial kinetic energy is transformed into surface energy, as it apparent from the capillary waves developing at the interface of the impacting drop whereas, for We &10, the drop spreading process is dominated by inertia. It is discussed in Wilderman et al. (2016) that approximately one-half of the initial kinetic energy is dissipated in the sudden expansion connecting the liquid sheet with the rim for We ⩾10, a fact which is implicitly taken into account in the mass and momentum conservation equations (3.2). In spite of the fact that the system of equations (3.2) is analogous to that used in Roisman (2009); Eggers et al. (2010); Villermaux & Bossa (2011); Lastakowski et al. (2014), the crucial difference between the present study and previous ones which is that, in our case, we do not impose a specific form of the velocity field within the lamella. Instead, since fluid particles conserve their velocities within the liquid sheet, our system of equations (3.2) subjected to the real initial conditions given in (2.3), with tecalculated self-consistently through equation (2.2), provides with the correct expressions for both the liquid velocity u(rt) and the height of the lamella h(rt) upstream of the toroidal rim, as it can be inferred from the good agreement existing between experiments and predictions depicted in figure 6b. Another of the advantages of not imposing an ad-hoc velocity field within the lamella rests on the fact that both the critical Weber number and the diameters and the velocities of the droplets ejected for impact velocities V > V ∗, which strongly depend on h(rt), can be calculated self-consistently. Indeed, in RG15, following the ideas in Villermaux & Bossa (2011); Agbaglah et al. (2013), we deduced a criterium for the disintegration of the edge of the liquid sheet based on the fact that the droplets composing the spray result from the amplification in the azimuthal direction of Rayleigh–Taylor and capillary instabilities. The growth rates of the disturbances developing in the azimuthal direction of the toroidal rim are highly attenuated as a consequence of the simultaneous growth of its thickness. In consequence, we concluded that drops will only be ejected when the time characterizing the radial growth of the rim, Th= (R/V )th= (1/HtdHt/dT)−1, is substantially larger than the capillary time Tc= (R/V )tc=ρ H3 t/8σ1/2. The breakup time tbis fixed at the instant at which tc/th≃0.085 for the reasons explained in RG15, which are reproduced here for the sake of clarity: i) the characteristic time of growth of a capillary