Initial surface deformations during impact on a liquid pool
Wilco Bouwhuis, Maurice H.W. Hendrix, Devaraj van der Meer, Jacco H. Snoeijer
aa r X i v : . [ phy s i c s . f l u - dyn ] M a r Under consideration for publication in J. Fluid Mech. Initial surface deformations during impacton a liquid pool
Wilco Bouwhuis , Maurice H.W. Hendrix , , Devaraj van der Meer ,and Jacco H. Snoeijer , Physics of Fluids Group, Faculty of Science and Technology, University of Twente, 7500 AEEnschede, The Netherlands, Laboratory for Aero and Hydrodynamics, Delft University of Technology, Leeghwaterstraat21, NL-2628 CA Delft, The Netherlands, Mesoscopic Transport Phenomena, Eindhoven University of Technology, Den Dolech 2, 5612AZ Eindhoven, The Netherlands(Received 30 September 2018)
A tiny air bubble can be entrapped at the bottom of a solid sphere that impacts ontoa liquid pool. The bubble forms due to the deformation of the liquid surface by a localpressure buildup inside the surrounding gas, as also observed during the impact of aliquid drop on a solid wall. Here we perform a perturbation analysis to quantitativelypredict the initial deformations of the free surface of the liquid pool as it is approachedby a solid sphere. We study the natural limits where the gas can be treated as a viscousfluid (Stokes flow) or as an inviscid fluid (potential flow). For both cases we derive thespatio-temporal evolution of the pool surface, and recover some of the recently proposedscaling laws for bubble entrapment. When inserting typical experimental values for theimpact parameters, we find that the bubble volume is mainly determined by the effectof gas viscosity.
1. Introduction
The phenomena resulting from solid-body impacts on liquid surfaces are widely stud-ied because of their omnipresence in nature and industry (Korobkin & Pukhnachov 1988;Howison et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al. et al.
W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer h(r,t)δ(r , t) rz U (a) (b) (c) h (t) R gasliquid entrapped air bubble Figure 1. (a) A solid sphere (radius R ) approaches a liquid surface with velocity U . The gapheight between the bottom of the sphere and the undisturbed water level ( z = 0) is h ( r, t ), r and t being the radial coordinate and time, respectively, with h (0 , t ) = h ( t ). (b) While the spheremoves downwards, the pool deflects by a small amount δ ( r, t ), as a result of the local pressurebuild-up in the air that is squeezed out. In the limit where δ ≪ h , which typically is valid upto very close to the impact time, the profiles are computed analytically. (c) This mechanism willfinally result in air bubble entrapment. et al. et al. et al. H d , and the entrapped bubble volume V b , respectively scale as H d ∼ R St − / , V b ∼ R St − / . (1.1)Here St is the Stokes number, St = ρ l U R/η g , in which ρ l is the density of the liquid, R is the radius of the drop, U is the impact velocity, and η g is the dynamic viscosity ofthe air. This scaling has been confirmed experimentally and numerically (Marston et al. et al. et al. et al. et al. H d ∼ R Ca / , V b ∼ R Ca , (1.2)where Ca = η g U/γ is the capillary number based on the gas properties and surfacetension γ . The crossover in between the two regimes, at which the size of the entrappedair bubble is maximal, is found by equating the predictions for H d from Eqs. (1.1) and(1.2). Then, one finds U ∼ η / g γ / / (cid:16) ρ / l R / (cid:17) , where U is the crossover impactvelocity, leading to maximal bubble entrapment. For an impacting water drop having aradius of 1 mm, this gives 0 .
07 m / s. Indeed, this is of the same order of magnitude as wasobserved experimentally, where the maximum bubble size was found around 0 .
25 m / s(for ethanol drops) (Bouwhuis et al. nitial surface deformations during impact on a liquid pool et al. (2012) it was experimentally found that, in the inertial regime, the bubble volume wasfixed before the rupture of the air film.In this paper, we analytically compute the initial deformations due to sphere impactonto a liquid pool in the inertial regime, where the deflection of the liquid is limited by itsinertia rather than by its surface tension. In experiments, there is generally not enoughresolution to accurately detect these initial deformations, and therefore we use numericalsimulations, bridging also towards larger deformations. By restricting ourselves to smalldeformations of the pool surface, we obtain detailed spatio-temperal information of thedeflection as well as the dependence on experimental parameters. This provides a naturalbridge between scaling theory, which lacks detailed information on the structure of theinterface deflection, and profiles obtained by direct numerical simulations. Similar calcu-lations were previously performed by Yiantsios & Davis (1990) in the capillary regime,recovering the scaling (1.2). Hence, such a small-deformation theory gives an analyticalfoundation to the scaling laws, as well as detailed predictions for the shape of the defor-mation. Although the problem of a cushioning air layer has been solved by Wilson (1991)for an ‘inertial’ air layer, a similar insightful similarity analysis for the inertial (liquid)regime was not yet attempted.The paper is organized as follows. Sec. 2 starts with a dimensional analysis of theproblem and shows the limiting cases when the gas can be described as a potential flowor as a viscous lubrication flow. This section also outlines the formalism based on whichthe interface deformations are computed. In Sec. 3 we present the results for both viscousgas flow and potential gas flow. The analytical results are illustrated for a representativecase of impact on a pool of water, with a sphere of radius R = 1 mm and velocity U = 5 m / s, surrounded by air, as typical in experiments (inertial regime). Here we alsoprovide a detailed comparison of our results with numerical simulations based on theBoundary Integral (BI) method, to validate our analysis and to investigate when theresults start to deviate from the small-deformation regime. In Sec. 4 we conclude on theresults in terms of air bubble entrapment.
2. Formulation
The geometry of the problem is sketched in Figure 1: we consider a solid sphere (radius R ) moving downwards towards the pool with a velocity U (Figure 1a). The velocity ofthe sphere during its fall is assumed to be constant, i.e. we neglect the acceleration ofgravity and the possible deceleration due to the gas flow. The movement of the air inducesan increase of the gas pressure at the bottom of the sphere, which will then deflect thepool surface by a distance δ ( r, t ) (Figure 1b). The deformation δ is defined positive whenthe pool deflects downwards. For as long as the interface deflection is small with respectto the height of the gap, i.e. | δ | ≪ h , the problem can be solved by a perturbationanalysis. In this section we first address the problem by dimensional analysis, and thenprovide the linearized formalism that allows computing the spatio-temporal evolution ofthe deflection δ ( r, t ). 2.1. Dimensional analysis
Let us first consider the gas flow induced by the motion of the sphere. In the regime wherethe height of the gap is much larger than R , the sphere does not experience any influence W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer of the pool. In that case, the Reynolds number of the gas flow is Re g = ρ g U R/η g , where ρ g is the density of air (1 .
204 kg / m ). However, as soon as the gap height becomes small, h /R ≪
1, the air flow will be oriented mainly in the radial direction. As is typical forlubrication flows (Reynolds 1886), one then has to consider a different Reynolds numberthat is obtained from the radial component of the Navier Stokes equation. In termsof scaling laws this gives ρ g u r /L ∼ η g u r /h , where u r is the typical radial gas flowvelocity, and L = √ Rh is the length scale in the radial direction (Hicks et al. et al. et al. U L ∼ u r h , and after elimination of u r one thus finds therelevant Reynolds number Re g, lubr. = ρ g U h /η g . In the thin-gap regime, the relativeinfluence of the viscosity and the inertia of the gas thus involves the gap thickness h instead of the sphere radius R .It is instructive to evaluate these parameters for typical experimental values, such asspheres falling in air ( ρ g = 1 .
204 kg / m , η g = 1 . × − Pa s) with R = 1 mm and U = 5 m / s. The crossover from inertial to viscous gas flow, Re g, lubr. ∼
1, arises when h ∼ µ m. This implies that there exists an “inertial thin-gap regime”, where h /R < g, lubr. > h < µ m,the gas can be described by a purely viscous flow. In the remainder, we therefore considera potential flow analysis during two parts of the trajectory: the large-gap stage h /R ≫ h /R ≪
1. The viscous flow is treated only in the final stages ofimpact, for which h /R ≪ U / ( gR ) is much larger than 1.2.2. From gas pressure to interface deflection
The first step of the analysis is to compute the response of the liquid on a gas pressure P g for the different limiting cases (viscous/inertial gas), as discussed above. Since we setout to compute the initial deformation, we can compute P g assuming the liquid pool isundeformed – the influence of a finite deflection is a correction at higher order in δ/h . Weassume axisymmetry and solve the equations in cylindrical coordinates ( r, z ) (see Figure1). The gas pressure will provide the boundary condition at the liquid pool, generatinga liquid flow as described by the linearized Euler equation: ∂~v∂t = − ρ l ~ ∇ P l , (2.1)where ~v ( r, z, t ) is the velocity field in the liquid, and P l ( r, z, t ) is the pressure inside theliquid. The advection terms in the Euler equation are quadratic in velocity and thereforeof higher order in δ/h , in analogy to the wave analysis (Lamb 1957). Without surfacetension, the gas pressure provides the boundary condition for the liquid pressure P l ( r, z = − δ, t ) ≃ P l ( r, z = 0 , t ) = P g ( r, t ) , (2.2)with the first equality again due to taking into account only leading order terms in δ/h .The resulting deflection is given by the kinematic boundary condition: nitial surface deformations during impact on a liquid pool ∂δ∂t = − v z | z = − δ − v r | z = − δ ∂δ∂r ≃ − v z | z =0 , (2.3)where v z | z =0 is the vertical velocity at the pool surface (to the lowest order in δ/h ).Substituting condition (2.3) into the vertical component of Eq. (2.1) gives ∂ δ∂t = 1 ρ l ∂P l ∂z (cid:12)(cid:12)(cid:12)(cid:12) z =0 . (2.4)The above equation shows that in order to compute δ ( r, t ), one requires a spatialderivative ∂P l /∂z . Hence, we need to find the pressure distribution inside the liquid thatis induced by P g at the free surface. For an incompressible liquid this can be achievedby taking the divergence of Eq. (2.1), which owing to ~ ∇ · ~v = 0 reduces to ∇ P l = 0.As the boundary condition is axisymmetric, it is natural to express the pressure as theaxisymmetric solution of the Laplace equation: P l ( r, z, t ) = ∞ Z c P g ( k, t )J ( kr )e kz k d k, (2.5)where the integration variable k is the wave number, and J ( kr ) is the Bessel function ofthe first kind with order ν = 0. The amplitude of the ‘modes’ J ( kr )e kz is given by theHankel transform of order 0 of the gas pressure P g ( r, t ), c P g ( k, t ) = ∞ Z P g ( r, t )J ( kr ) r d r. (2.6)Substituting this expression for the pressure into Eq. (2.4) gives ∂ δ∂t ( r, t ) = ∞ Z c P g ( k, t ) ρ l J ( kr ) k d k, (2.7)where we note an additional factor k coming from the derivative of ∂P l /∂z .The basic procedure for determining ∂ δ/∂t from the gas pressure is now clear: oneneeds to find the Hankel transform of the gas pressure (Eq. (2.6)), subsequently take thederivative of the result in the z -direction and evaluate the expression at z = 0, and finallytake the inverse Hankel transform (Eq. (2.7)). In the following section we will performthese steps for the gas pressures computed in the limits of Stokes gas flow and inviscidgas flow.
3. Results
Stokes gas flow
We now turn to the Stokes flow in the lubrication limit, which is valid for Re g, lubr. ≪ h /R ≪
1. In case of vanishing interface deformation, the gas pressure building upbelow an impacting sphere becomes (Davis et al. P g ( r, t ) = 3 η g U Rh (cid:16) r Rh (cid:17) = 3 η g UR (cid:18) RL (cid:19) F ( u ) . (3.1) W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer
Here we factorized the result in dimensional parameters determining the magnitude ofthe pressure and a dimensionless function F ( u ) that contains the spatial informationof the pressure profile. For this, we introduced L ( t ) = p Rh ( t ) as the relevant radiallength scale, while the geometrical function reads F ( u ) = 1 (cid:0) u (cid:1) ; u ( t ) = rL ( t ) . (3.2)Note that in the limit of vanishing thickness h , the pressure tends to diverge, P g ∼ h − , while the width of the peak becomes increasingly small, L ∼ h / . These singulartendencies are regularized when the deformations of the surface become comparable to h , but yet, set the characteristic scales for the enclosed bubble volume.We continue the analysis by inserting the gas pressure profile in Eq. (2.7), and find aclosed form expression: ∂ δ∂t ( r, t ) = 3 η g Uρ l RL (cid:18) RL (cid:19) G ( u ) . (3.3)Once more we recognize a dimensional prefactor that determines the scale of the acceler-ation, while the time-dependence follows from L ( t ) and u ( t ), and the spatial dependencethrough G ( u ). The additional factor 1 /L appearing in (3.3) originates from the scaling u = r/L . The spatial similarity profile is G ( u ) = ∞ R c F J ( ku ) k d k , where c F ( k ) is theHankel transform of F ( u ). The analytical expression for c F ( k ) is found to be c F ( k ) = √ k K ( √ k ) , (3.4)where K ( k ) is the modified Bessel function of the second kind with order ν = 1, andthe analytical expression for G ( u ) is G ( u ) = − (cid:16) u √ u +2 (cid:17) − E (cid:16) u √ u +2 (cid:17) + 14E (cid:16) u √ u +2 (cid:17) ( u + 2) / . (3.5)K and E are the complete elliptic integrals of the first and second kind, respectively.To illustrate and validate our analysis, we compare the predicted profiles with the re-sults obtained by Boundary Integral (BI) simulations (Pozrikidis 1997; Oguz & Prosperetti1993; Bergmann et al. et al. (2012, 2013): the liquid within the pool is described as a potential flow, while the pres-sure along the pool surface is explicitly calculated from the viscous lubrication equationfor the gas flow. To be able to confirm our theoretical predictions in the inertial regimewithout the influences of surface tension and hydrostatics (which are both very small, asmentioned in the Introduction), γ and g are equal to zero in our simulations. In the limitof small deflection, the simulations should thus recover Eq. (3.3).Figure 2a shows the configuration on the scale of the sphere, for typical impact pa-rameters for a sphere in air ( R = 1 mm, U = 5 m / s). The interface deflection δ is shownin Figure 2b, at the moment when the sphere is at a height h = 100 µ m. At this time, δ ≪ h ≪ R , for which we expect agreement between the BI results and our predictionfrom Eq. (3.3). Figure 2c shows the acceleration ∂ δ/∂t versus r . The solid line is theresult from the BI simulations and indeed gives perfect agreement with the prediction,represented by the dots. nitial surface deformations during impact on a liquid pool δ ( μ m ) −3 −2 −1 0 1 2 3020406080100120 r (mm) ∂ δ / ∂ t ( m / s ) TheoryBI (a)(b)(c) z ( mm ) Figure 2.
Deflection of the pool interface for Stokes gas flow; R = 1 mm, U = 5 m / s, startingheight of the (bottom of the) sphere h s = 0 . h = 0 . δ as a function of r , and (c) ∂ δ/∂t as a function of r . The solid red lines result from the Boundary Integral (BI) simulation. Thetheoretical result from Eq. (3.3) has been superimposed in panel c (blue dots). Note the differencein scales on the vertical axes of panel a and b. The BI results agree perfectly with the theoreticalpredictions, as long as | δ | ≪ h . −9 −7 −5 −3 −10 −8 −6 −4 h = h s − U t (m) δ | r = ( m ) TheoryBI δ = h Figure 3.
Deflection of the pool interface on the axis, δ r =0 , plotted against h ( t ), for Stokesgas flow; R = 1 mm, U = 5 m / s, starting height h s = 0 . h , the deflection δ | r =0 converges towards a -1/2 power-law. The BIresults perfectly agree with the theoretical predictions, until δ and h become of comparablemagnitude, pointed out by the crossing with the solid gray line δ | r =0 = h . At that moment δ r =0 saturates to a constant value, which is the ‘dimple height’ H d of Bouwhuis et al. (2012). The actual deflection profile δ ( r, t ) can not be integrated explicitly from (3.3), due tothe time-dependence through L and u . However, we can derive δ | r =0 , the deflection ofthe pool surface on the axis, which does not involve L ( t ). Using that ∂/∂t = − U ∂/∂h ,we find W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer ∂ δ | r =0 ∂h = 3 η g G (0) ρ l U R (cid:18) RL (cid:19) = 3 η g G (0) ρ l U R (cid:18) Rh (cid:19) / , (3.6)where (3.5) implies G (0) = √ π . The solution of Eq. (3.6) for δ | r =0 is subject tostart-up effects as long as h ∼ h s , where h s is the initial height of the gap. If we let theinitial height h s → ∞ , we find δ | r =0 ≃ √ π η g ρ l U (cid:18) Rh (cid:19) / . (3.7)This predicts that the central height increases dramatically when h decreases, as δ ∼ h − / . Figure 3 shows the BI result for δ | r =0 against h (solid line), superimposed withthe theoretical predictions (dashed line, taking into account the finite initial height h s ).Indeed, as soon as h ≪ h s , δ | r =0 converges to a − / δ ∼ h (indicated by thesolid gray line) and the lubrication approximation ceases to be valid. At this point, thedeflection converges to a constant, which will be the final dimple height H d . As statedin the Introduction, this will determine the dimple volume, and thus the entrapped airbubble volume, independently of the air film rupture process.The current analysis provides a rigorous foundation for the scaling results obtained pre-viously in Marston et al. (2011); Hicks & Purvis (2011); Hicks et al. (2012); Mandre & Brenner(2012); Bouwhuis et al. (2012). There, the ‘dimple height’ H d was observed to approacha constant value during the final stages of the impact. Figure 3 shows that this heightcan be estimated from δ r =0 ∼ h ∼ H d . Using (3.7), this immediately gives H d ∼ η g R / ρ l U H / d ∼ R St − / , (3.8)where St = ρ l U R/η g is the Stokes number. The corresponding volume of the entrappedbubble volume then scales as V b ∼ L H d ∼ RH d ∼ R St − / , (3.9)where we use the common estimate that L sets the lateral scale of the bubble. Theseare precisely the scaling predictions for the inertial regime (for Stokes gas flow), wherethe assumptions H d ∼ δ and L ∼ ( H d R ) / were further validated (Marston et al. et al. et al. et al. Potential gas flow
As motivated in Sec. 2.1, the inertial phase of the impacting sphere consists of twodistinct stages: the large-gap regime h ≫ R and the thin-gap regime h ≪ R . Belowwe separately treat both limiting cases analytically. We furthermore perform a numericalpotential flow calculation for the full range of h /R , to validate the analysis and to showhow the two stages are connected.3.2.1. Large-gap regime: h ≫ R When the sphere is very far from the pool surface, the flow field can be described bythe well-known potential flow field around a moving sphere of radius R . The introductionof the (undeformed) pool surface, however, requires that the gas velocity has no vertical nitial surface deformations during impact on a liquid pool v z | z =0 = 0. This boundary condition can be satisfied using the ‘methodof images’, corresponding to two approaching spheres having radius R with approachingvelocity U towards a mirroring horizontal line ( z = 0). Applying the superposition of thepotentials for the two moving spheres, one obtains the potential φ ( r, z, t ) = U R ( z − R − h ) (cid:16) r + ( z − R − h ) (cid:17) / − z + R + h (cid:16) r + ( z + R + h ) (cid:17) / . (3.10)It is important to realize that the introduction of the second moving sphere not onlyinfluences the flow around z = 0, but also gives a small, unwanted velocity on theboundary of the original sphere. In the limit of very large gaps, R/h ≪
1, this correctionbecomes negligible and (3.10) gives the asymptotically correct potential.We now extract the gas pressure profile on the level of the pool surface z = 0, byapplying the unsteady Bernoulli equation: P g ( r, t ) = ρ g U "(cid:18) Rζ (cid:19) F ( u ) + 92 (cid:18) Rζ (cid:19) F ( u ) ≃ ρ g U (cid:18) Rζ (cid:19) F ( u ) . (3.11)Here, ζ ( t ) = R + h ( t ) = R + h s − U t , the radial direction is scaled as u ( t ) = r/ζ , whilethe spatial profiles are F ( u ) = 2 − u (1 + u ) / ; (3.12) F ( u ) = − u (1 + u ) . (3.13)Since (3.10,3.11) are only valid for h ≫ R , we only keep the dominant first term in(3.11). Note that the width of the pressure peak is now set by the scale ζ = h + R . Thiscan be contrasted with the width in the thin-gap limit, L = √ Rh , which becomes verynarrow.Next, from (3.11) we can compute the induced acceleration profile using (2.7): ∂ δ∂t ( r, t ) = ρ g U ρ l ζ (cid:18) Rζ (cid:19) G ( u ) . (3.14)One recognizes a dimensional prefactor that is separated from the spatio-temporal de-pendence. The function G ( u ) = ∞ R c F J ( ku ) k d k is the spatial similarity profile, where c F ( k ) is the Hankel-transform of F ( u ). For G ( u ) we did not find any analytical expres-sion, but one can numerically calculate the given integral (cf. Figure 4).Once again, we can analytically compute the behavior of the central deflection, δ | r =0 : ∂ δ | r =0 ∂h = ρ g G (0) ρ l R (cid:18) Rζ (cid:19) . (3.15)Recalling that ∂/∂h = ∂/∂ζ and ζ → R for h → R , this implies that the final δ r =0 scales as ρ g R/ρ l . In contrast to the result for viscous flow, the typical deformation versus h depends only on the density ratio ρ g /ρ l , but not on the impact velocity. While the0 W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer density ratio is typically small, we anticipate that the resulting deflection for a millimeter-sized sphere can be a few microns. This is actually comparable to typical deflections inthe viscous lubrication phase. However, the pool is not deformed locally over a smallwidth √ Rh , but over the scale of the entire sphere, and therefore it will be of littleconsequence for the formation of the dimple and the size of the entrapped air bubble.3.2.2. Thin-gap regime: h ≪ R In the inertial thin-gap limit, the gas is squeezed out mainly in the radial direction.To predict the pressure profile for this stage of the impact, we use the depth-integratedcontinuity equation (Snoeijer et al. et al. ∂h∂t + 1 r ∂∂r ( rhu r ) = 0 , (3.16)where u r ( r, t ) is the height-averaged radial gas velocity in the gap. Assuming a plug flowthat does not depend on the z -coordinate, this average simply gives u r ( r, t ) = u r ( r, t ).This analytical description is similar to what has been done by Wilson (1991), who alsostudied cushioning air-layers at solid-liquid impact in the inertial thin-gap regime, thoughin 2D Cartesian coordinates, for general shapes of the impacting solid. In the present case,the bottom of the impacting solid sphere can be described as h = h ( t ) + r / (2 R ), andthus, ∂h/∂t = ∂h /∂t = − U . Hence, we can integrate (3.16) to find u r = u r = U r h (cid:16) r Rh (cid:17) . (3.17)The velocity profile (3.17) has a local maximum at r = √ Rh , and vanishes for r = 0and r = ∞ . Substituting the profile into the radial component of the Euler equation andintegrating over r gives the gas pressure: P g ( r, t ) = ρ g U R h r Rh (cid:16) r Rh (cid:17) = ρ g U (cid:18) RL (cid:19) F ( u ) , (3.18)with L ( t ) = p Rh ( t ), u ( t ) = r/L , and F ( u ) = 1 + u (cid:0) u (cid:1) . (3.19)Note that the geometry of the thin-gap again gives rise to a highly localized pressureprofile of a width √ Rh . The gas pressure again tends to diverge as h →
0, but moreslowly than in the viscous case: the inertial gas pressure in the thin-gap-limit is propor-tional to 1 /h , in contrast to the more singular scaling for the viscous gas flow scenario,1 /h .From (2.7) we deduce the pool surface acceleration ∂ δ∂t ( r, t ) = ρ g U ρ l L (cid:18) RL (cid:19) G ( u ) , (3.20)where G ( u ) = ∞ R c F J ( ku ) k d k , with c F ( k ) the Hankel-transform of F ( u ). At the origin r = 0, this reduces to nitial surface deformations during impact on a liquid pool ∂ δ | r =0 ∂h = ρ g G (0)2 ρ l R (cid:18) RL (cid:19) . (3.21)Just like in case of the large-gap regime, the central deflection has no dependence onimpact velocity. Solving gives δ r =0 ∼ h / + integration constants. From this we con-clude that in the inertial thin-gap limit, the pressure tends to diverge for h →
0, butthe deflection δ converges . Contrarily to the final stages in the case of viscous gas flow,the inertial gas pressure is not sufficiently singular to induce a strongly enhanced deflec-tion. The integration constants depend on the full history of the impact process, whichthus involves the dynamics during the preceding large-gap regime. To predict the actualdeflection during the final stages of sphere impact, it is thus not sufficient to consider thelarge-gap or thin-gap regime of the potential gas flow problem, but requires numericalsimulation of the full impact process over all h /R .3.2.3. Numerical simulations
Simulating the potential gas flow impact process using the BI technique calls for adifferent approach with respect to the case of Stokes gas flow. The reason is that werequire the gas pressure over the full range of gap thickness, including h ∼ R , for whichno analytical solution for the gas pressure is available that can serve as a boundarycondition for the liquid pool. As a consequence, the gas phase must be also computednumerically, which we achieve using the Boundary Integral code. We thus need to runtwo separate simulations. The process is started by a BI simulation of a solid sphereimpacting towards an undeformed surface, with in between a potential gas flow. Fromthis simulation, the gas pressure profile along the pool surface ( z = 0) is extracted. In thesecond BI simulation, this pressure is applied on a deformable pool surface, from whichwe eventually determine the resulting pool deflections. This is again a valid method aslong as δ/h ≪
1. The pressure data is transmitted from the first simulation to the secondsimulation through an extensive data file. Note that by doing two separate simulations,one needs to take into account the different length scales during the impact process(for h = 10 mm →
100 nm), implying very sensitive local node spacings and timedependencies. This was achieved by adapting the node spacing and time steps to ensureconvergence of the numerical results.Figure 4a and b respectively show the configuration on the length-scale of the sphere,and the interface deflection δ ( r, t ) for R = 1 mm, U = 5 m / s, and h = h s = 10 mm(i.e., the large-gap regime). Figure 4c shows the acceleration profile at the correspond-ing time, and is observed to agree very well with the asymptotic result of Eq. (3.14)(blue dots). The very small difference between the BI result and the theoretical pre-dictions can be explained by the fact that h s /R = 10, implying an expected differenceof about 10% between the theory and the numerical simulations. We remark that thecorresponding deformation (Figure 4b) is very small, as we look to the very initial defor-mations in the start-up regime. The demand h /R ≫ r = 0. We confirmed that this instability has a numerical origin and that it does notinfluence the result on the scale of the deformations. The thin-gap regime is analyzed inFigure 5. We again find very good agreement between the analytical gas velocity profile(panel a) and the pressure profile (panel b) and the BI results (here h = 100 nm).2 W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer δ ( μ m ) −25 −12.5 0 12.5 25051015 x 10 −3 r (mm) ∂ δ / ∂ t ( m / s ) TheoryBI z ( mm ) (a)(b)(c) (a) x 10 −8 Figure 4.
Deflection of the pool interface for potential gas flow in the limit h ≫ R ; R = 1 mm, U = 5 m / s, h = h s = 10 mm (thus, h /R = 10). (a) Global view plot of the sphere and poolcontours, (b) δ against r , and (c) ∂ δ/∂t against r . The solid red lines result from the BIsimulation. The theoretical result from Eq. (3.14) has been superimposed in panel c (blue dots).Note the difference in scales on the vertical axes of panel a and b. The BI results are nicelyagreeing with the theoretical predictions, until h /R becomes of order 1. u r ( m / s ) BITheory 0 0.02 0.04 0.06 0.08 0.1012 x 10 r (mm) P ( P a ) BITheory0 0.02 0.04 0.06 0.08 0.1 r (mm) (a) (b) Figure 5.
Inertial gas flow in the thin-gap limit. Theoretical prediction (red dashed line) andBI gas flow simulation result (blue solid line) of the velocity (panel a) and pressure (panel b)profile within the gas; R = 1 mm, U = 5 m / s, h s = h = 100 nm (thus, h /R = 10 − ). We findvery good agreement between the theoretical predictions and the BI results. The crossover between the large-gap and thin-gap limits is illustrated in Figure 6,showing the gas pressure on the symmetry axis r =0. As predicted, in the limit h /R ≫ ρ g U ( R/ζ ) (red dashed line), and inthe limit h /R ≪ ρ g U R/ (2 h ) (green dashed line). This confirmsthe validity of the analytical approaches. Finally, we investigate the deflection of the poolthat is induced by the numerically obtained gas pressure. Figure 7 shows the deflectionat r =0, the inertial (gas) counterpart of Figure 3. As expected, δ r =0 deviates from the nitial surface deformations during impact on a liquid pool −7 −5 -3 -1 −2 h (m) P ( P a ) BI P r=0 = 2 ρ g U ( R / ζ ) = ½ ρ g U R / h P r=0 h = R Figure 6.
Behavior of the gas pressure on the axis, P r =0 , plotted against the gap height h ( t ),for potential gas flow. R = 1 mm, U = 5 m / s, h s = 10 mm. The red dashed line is the theoreticalprediction in the regime h ≫ R ; The green dashed line is the theoretical prediction in the regime h ≪ R . The dashed line points out the crossover h = R . The BI gas flow simulation result(blue solid line) indeed follows these predicted behaviors in the corresponding regimes, with acrossover at h ∼ R . −7 −5 −3 −1 −10 −8 −6 h = h s − U t (m) δ | r = ( m ) δ = h h = R Theory (large gap limit)BI
Figure 7.
Deflection of the pool interface on the axis, δ r =0 , plotted against the gap height h ( t ),for potential gas flow; R = 1 mm, U = 5 m / s, h s = 10 mm. The solid red line is the resultfrom the BI solid sphere on liquid pool simulations (Sec. 3.2.3). The theoretical result from Eq.(3.15) for the large-gap regime has been superimposed (blue dashed line). δ | r =0 saturates to aconstant. The BI results perfectly agree with the large-gap predictions in the regime h ≫ R .In the regime h ≪ R , δ r =0 deviates from this prediction, but the difference is relatively small.The dashed gray line points out the crossover h = R ; The solid gray line points out δ | r =0 = h . large-gap prediction in the small-gap regime, though the deviation is not very large.This means that, despite the fact that the gas pressure tends to diverge for h → ∂ δ/∂t , which in principlecould again be directly calculated from the pressure profile † . † A second reason is that, in the numerical simulations, the very small gap height of 100 nmneeds a very high local node density on both the pool surface and the sphere surface; The dif-ference in length scales of R and h is four decades, which is very challenging. This induces verysmall time steps to be able to calculate a fair second derivative of the deflection profile in time. Inaddition, the pressure along the pool surface needs to be extracted from a prior solid-sphere-on–solid-surface simulation (through an extensive data file), which makes the discretization morecomplicated. W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer
The large-gap prediction for the final δ r =0 is thus satisfactory, and we conclude withthe following scaling law for the resulting dimple height H d for the inertial gas scenarioas was concluded from Eq. (3.15): H d ∼ R ρ g ρ l . (3.22)This dimple height is independent of the impact velocity of the sphere. Since the surfacedeformation is the sum of the deformations in the both the large-gap and the thin-gaplimit, it is unclear what the correct radial and axial length scales are that lead to thevolume of the pinched bubble.
4. Conclusion
We performed a perturbation analysis to investigate the initial deflections of a liquidsurface, induced by the approach of an impacting solid sphere. The analysis assumes thedeflection is limited by the inertia of the liquid pool (i.e., not by its surface tension) and weconsider two natural limits for the surrounding medium: Stokes gas flow and potentialgas flow. We obtained a quantitative prediction for the pool surface deflection, whichwas validated numerically, and recovered previously proposed scaling laws for bubbleentrapment.While the ‘cushioning’ of an inertial gas layer had been analyzed before (Wilson 1991),most recent work on liquid or solid impact assumes a viscous gas layer. Surprisingly, ouranalysis reveals that inertial and viscous cushioning both lead to a pool deflection ofthe order 1 µ m, for typical experimental conditions. However, the Stokes gas pressurestrongly tends to diverge for h →
0, much more strongly than during the inertial gasphase. In addition, this viscous lubrication pressure profile is very localized, while mostof the inertial deflection is generated during the initial phase where the pool deflection isspread over the entire width of the sphere. This explains why the experimental results onbubble entrapment are in close agreement with the scaling law (3.9) (Tran et al. et al. et al. et al. et al. R , the dynamics willexhibit two different types of crossover: a geometric crossover based on the relative thick-ness of the gap h/R , and a crossover from inertial to viscous gas flow. The order in whichthese crossovers occur depends on the parameters of the problem. In our numerical ex-amples we assumed one first reaches the thin-gap regime, before the lubrication Reynoldsnumber (based on the gap thickness h ), becomes smaller than unity. This order can bereversed for impact at smaller velocities or for a sphere sinking in a more viscous medium.In that case, however, one needs to bear in mind that the influence of the pool surfacetension will become more important, corresponding to the capillary impact regime. Inthis case, the thin film potentially has time to drain out before a bubble is formed, makingthe entrapment process more complex (Klaseboer et al. et al. nitial surface deformations during impact on a liquid pool Acknowledgments
We gratefully acknowledge Stephen Wilson and Hanneke Gelderblom for insightfuldiscussion. This work was supported by STW and NWO through a VIDI Grant No.11304.
REFERENCESBergmann, R. P. H. M., van der Meer, D., Gekle, S., van der Bos, J. & Lohse, D.
J. Fluid Mech. ,381.
Bouwhuis, W., van der Veen, R. C. A., Tran, T., Keij, D. L., Winkels, K. G., Peters,I. R., van der Meer, D., Sun, C., Snoeijer, J. H. & Lohse, D.
Phys. Rev. Lett , 264501.
Bouwhuis, W., Winkels, K. G., Peters, I. R., Brunet, P., van der Meer, D. & Snoeijer,J. H.
Phys. Rev. E ,023017. van Dam, D. & Le Clerc, C. Phys. Fluids , 3403. Davis, R. H., Serayssol, J.-M. & Hinch, E. J.
J. Fluid Mech. , 479–497.
Deng, Q., Anilkumar, A. V. & Wang, T. G.
J. Coll. and Interf. Sc. , 523–532.
Do-Quang, M. & Amberg, G.
Phys. Fluids , 022102.
Driscoll, M. M. & Nagel, S. R.
Phys. Rev. Lett. , 154502.
Hicks, P. D., Ermanyuk, E. V., Gavrilov, N. V. & Purvis, R.
J. Fluid Mech. , 310–320.
Hicks, P. D. & Purvis, R.
Phys. Fluids , 062104. Howison, S. D., Ockendon, J. R. & Wilson, S. K.
J. Fluid Mech. , 215–230.
Klaseboer, E., Chevaillier, J. P., Gourdon, C. & Masbernat, O.
J. Colloid.Int. Sc. , 274–285.
Klaseboer, E., Manica, R. & Chan, D. Y. C.
Phys. Rev. Lett , 194501.
Korobkin, A. A., Ellis, A. S. & Smith, F. T.
J. Fluid Mech. , 365–394.
Korobkin, A. A. & Pukhnachov, V. V.
Ann. Rev. FluidMech. , 159–185. Lamb, H.
Hydrodynamics , 6th edn. Cambridge University Press.
Mandre, S. & Brenner, M. P.
J.Fluid Mech. , 148–172.
Marston, J. O., Vakarelski, I. U. & Thoroddsen, S. T.
J. Fluid Mech. , 660–670.
Moore, M. R. & Oliver, J. M.
IMA J.Appl. Math. , 661–680. Oguz, H. N. & Prosperetti, A.
J. Fluid Mech. , 111–145.
Pozrikidis, C.
Introduction to theoretical and computational fluid dynamics , 1st edn.Oxford University Press.
Reynolds, O.
Phil.Trans. R. Soc. Lond. , 157–234. W. Bouwhuis, M.H.W. Hendrix, D. van der Meer, and J.H. Snoeijer
Smith, F. T., Li, L. & Wu, G. X.
J.Fluid Mech. , 291–318.
Snoeijer, J. H., Brunet, P. & Eggers, J.
Phys. Rev. E , 036307. Thoroddsen, S. T., Etoh, T. G., Takehara, K., Ootsuka, N. & Hatsuki, A.
J. Fluid Mech. , 203.
Thoroddsen, S. T., Thoraval, M. J., Takehara, K. & Etoh, T. G.
J. Fluid Mech. , 469–479.
Tran, T., de Maleprade, H., Sun, C. & Lohse, D.
J. Fluid Mech. , R3.
Wilson, S. K.
J. Eng. Math. , 265–285. Yiantsios, S. G. & Davis, R. H.
J. Fluid Mech. , 547–573.
Yoon, Y., Borrell, M., Park, C. C. & Leal, L. G.
J. Fluid Mech.525