Hostname: page-component-6898cfcb55-cbk24 Total loading time: 0 Render date: 2026-11-11T04:45:08.736Z Has data issue: false hasContentIssue false
  • English
  • Français

Initial surface deformations during impact on a liquid pool

Published online by Cambridge University Press:  20 April 2015

Wilco Bouwhuis,
Maurice H. W. Hendrix,
Devaraj van der Meer and
Jacco H. Snoeijer
Wilco Bouwhuis*
Affiliation:
Physics of Fluids Group, Faculty of Science and Technology, University of Twente, 7500 AE Enschede, The Netherlands
Maurice H. W. Hendrix
Affiliation:
Physics of Fluids Group, Faculty of Science and Technology, University of Twente, 7500 AE Enschede, The NetherlandsLaboratory for Aero and Hydrodynamics, Delft University of Technology, Leeghwaterstraat 21, NL-2628 CA Delft, The Netherlands
Devaraj van der Meer
Affiliation:
Physics of Fluids Group, Faculty of Science and Technology, University of Twente, 7500 AE Enschede, The Netherlands
Jacco H. Snoeijer
Affiliation:
Physics of Fluids Group, Faculty of Science and Technology, University of Twente, 7500 AE Enschede, The NetherlandsMesoscopic Transport Phenomena, Eindhoven University of Technology, Den Dolech 2, 5612 AZ Eindhoven, The Netherlands
*
†Email address for correspondence: w.bouwhuis@ziggo.nl
Rights & Permissions [Opens in a new window]

Abstract

A tiny air bubble can be entrapped at the bottom of a solid sphere that impacts onto a liquid pool. The bubble forms due to the deformation of the liquid surface by a local pressure buildup inside the surrounding gas, as also observed during the impact of a liquid drop on a solid wall. Here, we perform a perturbation analysis to quantitatively predict the initial deformations of the free surface of a liquid pool as it is approached by a solid sphere. We study the natural limits where the gas can be treated as a viscous fluid (Stokes flow) or as an inviscid fluid (potential flow). For both cases we derive the spatiotemporal evolution of the pool surface, and recover some of the recently proposed scaling laws for bubble entrapment. On inserting typical experimental values for the impact parameters, we find that the bubble volume is mainly determined by the effect of gas viscosity.

JFM classification

Drops and Bubbles: DropsDrops and Bubbles: Drops and Bubbles

Information

Type
Papers
Information
Journal of Fluid Mechanics , Volume 771 , 25 May 2015, pp. 503 - 519
Check for updates
Copyright
© 2015 Cambridge University Press 

1. Introduction

The phenomena resulting from solid-body impacts on liquid surfaces are widely studied because of their omnipresence in nature and industry (Korobkin & Pukhnachov Reference Korobkin and Pukhnachov1988; Howison, Ockendon & Wilson Reference Howison, Ockendon and Wilson1991; Korobkin, Ellis & Smith Reference Korobkin, Ellis and Smith2008; Deng, Anilkumar & Wang Reference Deng, Anilkumar and Wang2009; Do-Quang & Amberg Reference Do-Quang and Amberg2009; Marston, Vakarelski & Thoroddsen Reference Marston, Vakarelski and Thoroddsen2011; Hicks et al. Reference Hicks, Ermanyuk, Gavrilov and Purvis2012; Moore & Oliver Reference Moore and Oliver2014). These involve splashing, jet formation, cavity formation and air bubble entrapment. The entrapment of tiny micrometre-sized air bubbles between a solid object and a pool is due to a mechanism similar to that of the impact of a liquid drop on a solid surface (Smith, Li & Wu Reference Smith, Li and Wu2003; van Dam & Le Clerc Reference van Dam and Le Clerc2004; Thoroddsen et al. Reference Thoroddsen, Etoh, Takehara, Ootsuka and Hatsuki2005; Driscoll & Nagel Reference Driscoll and Nagel2011; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Mandre & Brenner Reference Mandre and Brenner2012; Klaseboer, Manica & Chan Reference Klaseboer, Manica and Chan2014) or of a drop onto a liquid pool (Yiantsios & Davis Reference Yiantsios and Davis1990; Hicks & Purvis Reference Hicks and Purvis2011; Thoroddsen et al. Reference Thoroddsen, Thoraval, Takehara and Etoh2012; Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013). The air that surrounds the falling object is squeezed out between the solid and the pool surface during the final stages of impact, resulting in a local pressure buildup in the gas. This pressure will induce a small deformation of the liquid surface (figure 1 b), which will finally result in the entrapment of a tiny air bubble by the rupture of the enclosed air film (figure 1 c). For many applications these air bubbles are undesirable, and, hence, the prediction of their size is of great importance.

Figure 1.

(a) A solid sphere (radius $R$ ) approaches a liquid surface with velocity $U$ . The gap height between the bottom of the sphere and the undisturbed water level ( $z=0$ ) is $h(r,t)$ , where $r$ and $t$ are the radial coordinate and time respectively, with $h(0,t)=h_{0}(t)$ . (b) While the sphere moves downwards, the pool deflects by a small amount ${\it\delta}(r,t)$ , as a result of the local pressure buildup in the air that is squeezed out. In the limit where ${\it\delta}\ll h_{0}$ , which typically is valid up to very close to the impact time, the profiles are computed analytically. (c) This mechanism will finally result in air bubble entrapment.

There are two main types of theoretical approach to determine the deformations of the liquid surface and predict the size of the entrapped air bubble, namely full numerical solutions of the problem and scaling arguments (Wilson Reference Wilson1991; Hicks & Purvis Reference Hicks and Purvis2011; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Hicks et al. Reference Hicks, Ermanyuk, Gavrilov and Purvis2012; Mandre & Brenner Reference Mandre and Brenner2012; Klaseboer et al. Reference Klaseboer, Manica and Chan2014). Combined with experiments, these have led to the observation that the size of the air bubble for impact of a liquid drop on a flat solid surface is determined by either the inertia of the liquid or the surface tension (Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012). For increasingly high impact velocities, liquid inertia dominates and reduces the size of the air film at impact (‘inertial regime’), while surface tension dominates for lower velocities or smaller drop sizes (‘capillary regime’). The case of a solid sphere impacting on a pool leads to similar bubble entrapment, and, moreover, in the inertial regime the same scaling law (including the multiplicative prefactor) as for the impact of a drop on a solid has been observed (Marston et al. Reference Marston, Vakarelski and Thoroddsen2011; Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013). Here, the final centre height difference between the two surfaces, which is called the dimple height $H_{d}$ , and the entrapped bubble volume $V_{b}$ respectively scale as

(1.1a,b ) $$\begin{eqnarray}H_{d}\sim RSt^{-2/3},\quad V_{b}\sim R^{3}St^{-4/3}.\end{eqnarray}$$

Here, $St$ is the Stokes number, $St={\it\rho}_{l}UR/{\it\eta}_{g}$ , in which ${\it\rho}_{l}$ is the density of the liquid, $R$ is the radius of the drop, $U$ is the impact velocity and ${\it\eta}_{g}$ is the dynamic viscosity of the air. This scaling has been confirmed experimentally and numerically (Hicks & Purvis Reference Hicks and Purvis2011; Marston et al. Reference Marston, Vakarelski and Thoroddsen2011; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Hicks et al. Reference Hicks, Ermanyuk, Gavrilov and Purvis2012; Mandre & Brenner Reference Mandre and Brenner2012; Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013). On the other hand, in the capillary regime (small velocities or small drops), the scaling analysis predicts (Yiantsios & Davis Reference Yiantsios and Davis1990; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012)

(1.2a,b ) $$\begin{eqnarray}H_{d}\sim RCa^{1/2},\quad V_{b}\sim R^{3}Ca,\end{eqnarray}$$

where $Ca={\it\eta}_{g}U/{\it\gamma}$ is the capillary number based on the gas properties and surface tension ${\it\gamma}$ . The crossover between the two regimes, at which the size of the entrapped air bubble is maximal, is found by equating the predictions for $H_{d}$ from (1.1a,b ) and (1.2a,b ). Then, one finds $U_{0}\sim {\it\eta}_{g}^{1/7}{\it\gamma}^{3/7}/({\it\rho}_{l}^{4/7}R^{4/7})$ , where $U_{0}$ is the crossover impact velocity, leading to maximal bubble entrapment. For an impacting water drop having a radius of $1~\text{mm}$ , this gives $0.07~\text{m}~\text{s}^{-1}$ . Indeed, this is of the same order of magnitude as was observed experimentally, where the maximum bubble size was found at approximately $0.25~\text{m}~\text{s}^{-1}$ (for ethanol drops) (Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012). Generically, for drops or spheres falling at their terminal velocity of a few metres per second, the impact thus takes place in the inertial regime, where the effect of surface tension can thus be neglected when focusing on the air bubble entrapment. It should be noted that surface tension will enter during the rupture of the air film, which, however, appears to be on a different time scale. In Bouwhuis et al. (Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012) it was experimentally found that, in the inertial regime, the bubble volume was fixed before the rupture of the air film.

In this paper, we analytically compute the initial deformations due to sphere impact onto a liquid pool in the inertial regime, where the deflection of the liquid is limited by its inertia rather than by its surface tension. In experiments, there is generally not enough resolution to accurately detect these initial deformations, and therefore we use numerical simulations, which also provide a bridge towards larger deformations. By restricting ourselves to small deformations of the pool surface, we obtain detailed spatiotemporal information on the deflection as well as the dependence on experimental parameters. This provides a natural bridge between scaling theory, which lacks detailed information on the structure of the interface deflection, and profiles obtained by direct numerical simulations. Similar calculations were previously performed by Yiantsios & Davis (Reference Yiantsios and Davis1990) in the capillary regime, recovering the scaling (1.2a,b ). Hence, such a small-deformation theory gives an analytical foundation to the scaling laws, as well as detailed predictions of the shape of the deformation. Although the problem of a cushioning air layer has been solved by Wilson (Reference Wilson1991) for an ‘inertial’ air layer, a similar insightful similarity analysis for the inertial (liquid) regime has not yet been attempted.

This paper is organized as follows. Section 2 starts with a dimensional analysis of the problem and shows the limiting cases when the gas can be described as a potential flow or as a viscous lubrication flow. This section also outlines the formalism based on which the interface deformations are computed. In § 3 we present the results for both viscous gas flow and potential gas flow. The analytical results are illustrated for a representative case of impact on a pool of water, with a sphere of radius $R=1~\text{mm}$ and velocity $U=5~\text{m}~\text{s}^{-1}$ , surrounded by air, as is typical in experiments (inertial regime). Here, we also provide a detailed comparison of our results with numerical simulations based on the boundary integral (BI) method, to validate our analysis and to investigate when the results start to deviate from the small-deformation regime. In § 4 we conclude on the results 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 a pool with a velocity $U$ (figure 1 a). The velocity of the sphere during its fall is assumed to be constant, i.e. we neglect the acceleration of gravity and the possible deceleration due to the gas flow. The movement of the air induces an increase of the gas pressure at the bottom of the sphere, which will then deflect the pool surface by a distance ${\it\delta}(r,t)$ (figure 1 b). The deformation ${\it\delta}$ is defined as positive when the pool deflects downwards. For as long as the interface deflection is small with respect to the height of the gap, i.e. $|{\it\delta}|\ll h$ , the problem can be solved by a perturbation analysis. In this section we first address the problem by dimensional analysis, and then provide the linearized formalism that allows computation of the spatiotemporal evolution of the deflection ${\it\delta}(r,t)$ .

2.1. Dimensional analysis

Let us first consider the gas flow induced by the motion of the sphere. In the regime where the height of the gap is much larger than $R$ , the sphere does not experience any influence of the pool. In this case, the Reynolds number of the gas flow is $Re_{g}={\it\rho}_{g}UR/{\it\eta}_{g}$ , where ${\it\rho}_{g}$ is the density of air ( $1.204~\text{kg}~\text{m}^{-3}$ ). However, as soon as the gap height becomes small, $h_{0}/R\ll 1$ , the air flow will be oriented mainly in the radial direction. As is typical for lubrication flows (Reynolds Reference Reynolds1886), one then has to consider a different Reynolds number which is obtained from the radial component of the Navier Stokes equation. In terms of scaling laws this gives ${\it\rho}_{g}u_{r}^{2}/L\sim {\it\eta}_{g}u_{r}/h_{0}^{2}$ , where $u_{r}$ is the typical radial gas flow velocity and $L=\sqrt{Rh_{0}}$ is the length scale in the radial direction (Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Hicks et al. Reference Hicks, Ermanyuk, Gavrilov and Purvis2012; Mandre & Brenner Reference Mandre and Brenner2012; Klaseboer et al. Reference Klaseboer, Manica and Chan2014). Application of mass conservation on the air gives $UL\sim u_{r}h_{0}$ , and after elimination of $u_{r}$ one thus finds the relevant Reynolds number $Re_{g,lubr.}={\it\rho}_{g}Uh_{0}/{\it\eta}_{g}$ . In the thin-gap regime, the relative influence of the viscosity and the inertia of the gas thus involves the gap thickness $h_{0}$ instead of the sphere radius $R$ .

It is instructive to evaluate these parameters for typical experimental values, such as spheres falling in air ( ${\it\rho}_{g}=1.204~\text{kg}~\text{m}^{-3}$ , ${\it\eta}_{g}=1.82\times 10^{-5}~\text{Pa}~\text{s}$ ) with $R=1~\text{mm}$ and $U=5~\text{m}~\text{s}^{-1}$ . The crossover from inertial to viscous gas flow, $Re_{g,lubr.}\sim 1$ , arises when $h_{0}\sim 3~{\rm\mu}\text{m}$ . This implies that there exists an ‘inertial thin-gap regime’, where $h_{0}/R<1$ and $Re_{g,lubr.}>1$ at the same time. Only for the final stages of the impact, $h_{0}<3~{\rm\mu}\text{m}$ , can the gas be described by a purely viscous flow. In the remainder of this paper, we therefore consider a potential flow analysis during two parts of the trajectory: the large-gap stage $h_{0}/R\gg 1$ and the thin-gap stage $h_{0}/R\ll 1$ . The viscous flow is treated only in the final stages of impact, for which $h_{0}/R\ll 1$ and it is thus justified to resort to lubrication theory. The various limits will be worked out separately in § 3.

The liquid pool is assumed to be a low-viscosity liquid and is treated for small-amplitude deformations. These are essentially the same assumptions as for the propagation of linear surface waves (Lamb Reference Lamb1957), where the amplitude is small with respect to the length scales of the problem. We focus on the ‘inertial regime’ of impact, where the deformation is limited by the acceleration of the liquid and not by the surface tension of the liquid–air interface. The influence of gravity will also be neglected in the theory; the Froude number based on the impact parameters, $Fr=U^{2}/(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 set out to compute the initial deformation, we can compute $P_{g}$ assuming that the liquid pool is undeformed – the influence of a finite deflection is a correction at higher order in ${\it\delta}/h$ . We assume axisymmetry and solve the equations in cylindrical coordinates $(r,z)$ (see figure 1). The gas pressure will provide the boundary condition at the liquid pool, generating a liquid flow as described by the linearized Euler equation:

(2.1) $$\begin{eqnarray}\frac{\partial \boldsymbol{v}}{\partial t}=-\frac{1}{{\it\rho}_{l}}\boldsymbol{{\rm\nabla}}P_{l},\end{eqnarray}$$

where $\boldsymbol{v}(r,z,t)$ is the velocity field in the liquid and $P_{l}(r,z,t)$ is the pressure inside the liquid. The advection terms in the Euler equation are quadratic in velocity and therefore of higher order in ${\it\delta}/h$ , in analogy to the wave analysis (Lamb Reference Lamb1957). Without surface tension, the gas pressure provides the boundary condition for the liquid pressure,

(2.2) $$\begin{eqnarray}P_{l}(r,z=-{\it\delta},t)\simeq P_{l}(r,z=0,t)=P_{g}(r,t),\end{eqnarray}$$

with the first equality again due to taking into account only leading order terms in ${\it\delta}/h$ . The resulting deflection is given by the kinematic boundary condition:

(2.3) $$\begin{eqnarray}\frac{\partial {\it\delta}}{\partial t}=-v_{z}|_{z=-{\it\delta}}-v_{r}|_{z=-{\it\delta}}\frac{\partial {\it\delta}}{\partial r}\simeq -v_{z}|_{z=0},\end{eqnarray}$$

where $v_{z}|_{z=0}$ is the vertical velocity at the pool surface (to the lowest order in ${\it\delta}/h$ ). Substitution of condition (2.3) into the vertical component of (2.1) gives

(2.4) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}}{\partial t^{2}}=\frac{1}{{\it\rho}_{l}}\left.\frac{\partial P_{l}}{\partial z}\right|_{z=0}.\end{eqnarray}$$

The above equation shows that in order to compute ${\it\delta}(r,t)$ , one requires a spatial derivative $\partial P_{l}/\partial z$ . Hence, we need to find the pressure distribution inside the liquid that is induced by $P_{g}$ at the free surface. For an incompressible liquid this can be achieved by taking the divergence of (2.1), which due to $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}=0$ reduces to ${\rm\nabla}^{2}P_{l}=0$ . As the boundary condition is axisymmetric, it is natural to express the pressure as the axisymmetric solution of the Laplace equation:

(2.5) $$\begin{eqnarray}P_{l}(r,z,t)=\int _{0}^{\infty }\widehat{P_{g}}(k,t)\text{J}_{0}(kr)\text{e}^{kz}k\,\text{d}k,\end{eqnarray}$$

where the integration variable $k$ is the wavenumber, and $\text{J}_{0}(kr)$ is the Bessel function of the first kind with order ${\it\nu}=0$ . The amplitude of the ‘modes’ $\text{J}_{0}(kr)\text{e}^{kz}$ is given by the Hankel transform of order 0 of the gas pressure $P_{g}(r,t)$ ,

(2.6) $$\begin{eqnarray}\widehat{P_{g}}(k,t)=\int _{0}^{\infty }P_{g}(r,t)\text{J}_{0}(kr)r\,\text{d}r.\end{eqnarray}$$

Substitution of this expression for the pressure into (2.4) gives

(2.7) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}}{\partial t^{2}}(r,t)=\int _{0}^{\infty }\frac{\widehat{P_{g}}(k,t)}{{\it\rho}_{l}}\text{J}_{0}(kr)k^{2}\,\text{d}k,\end{eqnarray}$$

where we note an additional factor $k$ coming from the derivative of $\partial P_{l}/\partial z$ .

The basic procedure for determining $\partial ^{2}{\it\delta}/\partial t^{2}$ from the gas pressure is now clear: one needs to find the Hankel transform of the gas pressure (2.6), subsequently take the derivative of the result in the $z$ direction and evaluate the expression at $z=0$ , and finally take the inverse Hankel transform (2.7). In the following section we will perform these steps for the gas pressure computed in the limits of Stokes gas flow and inviscid gas flow.

3. Results

3.1. Stokes gas flow

We now turn to the Stokes flow in the lubrication limit, which is valid for $Re_{g,lubr.}\ll 1$ and $h_{0}/R\ll 1$ . In the case of vanishing interface deformation, the gas pressure building up below an impacting sphere becomes (Davis, Serayssol & Hinch Reference Davis, Serayssol and Hinch1986; Yiantsios & Davis Reference Yiantsios and Davis1990)

(3.1) $$\begin{eqnarray}P_{g}(r,t)=\frac{3{\it\eta}_{g}UR}{h_{0}^{2}\left(1+\displaystyle \frac{r^{2}}{2Rh_{0}}\right)^{2}}=\frac{3{\it\eta}_{g}U}{R}\left(\frac{R}{L}\right)^{4}F_{1}(u).\end{eqnarray}$$

Here, we factorized the result in dimensional parameters, determining the magnitude of the pressure and a dimensionless function $F_{1}(u)$ that contains the spatial information on the pressure profile. For this, we introduced $L(t)=\sqrt{Rh_{0}(t)}$ as the relevant radial length scale, while the geometrical function reads

(3.2a,b ) $$\begin{eqnarray}F_{1}(u)=\frac{1}{\left(1+{\textstyle \frac{1}{2}}u^{2}\right)^{2}};\quad u(t)=\frac{r}{L(t)}.\end{eqnarray}$$

It should be noted that in the limit of vanishing thickness $h_{0}$ , the pressure tends to diverge, $P_{g}\sim h_{0}^{-2}$ , while the width of the peak becomes increasingly small, $L\sim h_{0}^{1/2}$ . These singular tendencies are regularized when the deformations of the surface become comparable to $h_{0}$ , but still set the characteristic scales for the enclosed bubble volume.

We continue the analysis by inserting the gas pressure profile in (2.7), and find a closed form expression:

(3.3) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}}{\partial t^{2}}(r,t)=\frac{3{\it\eta}_{g}U}{{\it\rho}_{l}RL}\left(\frac{R}{L}\right)^{4}G_{1}(u).\end{eqnarray}$$

Once more we recognize a dimensional prefactor that determines the scale of the acceleration, while the time dependence follows from $L(t)$ and $u(t)$ , and the spatial dependence through $G_{1}(u)$ . The additional factor $1/L$ appearing in (3.3) originates from the scaling $u=r/L$ . The spatial similarity profile is $G_{1}(u)=\int _{0}^{\infty }\widehat{F_{1}}\text{J}_{0}(ku)k^{2}\,\text{d}k$ , where $\widehat{F_{1}}(k)$ is the Hankel transform of $F_{1}(u)$ . The analytical expression for $\widehat{F_{1}}(k)$ is found to be

(3.4) $$\begin{eqnarray}\displaystyle \widehat{F_{1}}(k)=\sqrt{2}k\text{K}_{1}(\sqrt{2}k), & & \displaystyle\end{eqnarray}$$

where $\text{K}_{1}(k)$ is the modified Bessel function of the second kind with order ${\it\nu}=1$ , and the analytical expression for $G_{1}(u)$ is

(3.5) $$\begin{eqnarray}G_{1}(u)=\frac{-8\text{K}\left(\displaystyle \frac{u}{\sqrt{u^{2}+2}}\right)-\text{E}\left(\displaystyle \frac{u}{\sqrt{u^{2}+2}}\right)+14\text{E}\left(\displaystyle \frac{u}{\sqrt{u^{2}+2}}\right)}{\left(u^{2}+2\right)^{5/2}}.\end{eqnarray}$$

K and E are the complete elliptic integrals of the first and second kind respectively.

Figure 2.

Deflection of the pool interface for Stokes gas flow; $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ , starting height of the (bottom of the) sphere $h_{s}=0.5~\text{mm}$ , current height $h_{0}=0.1~\text{mm}$ . (a) Global view of the sphere and pool contours, (b) the pool deflection ${\it\delta}$ as a function of $r$ and (c) $\partial ^{2}{\it\delta}/\partial t^{2}$ as a function of $r$ . The solid red lines result from the BI simulation. The theoretical result from (3.3) has been superimposed in (c) (blue dots). Note the difference in scales on the vertical axes of (a) and (b). The BI results agree perfectly with the theoretical predictions, as long as $|{\it\delta}|\ll h$ .

Figure 3.

Deflection of the pool interface on the axis, ${\it\delta}_{r=0}$ , plotted against $h_{0}(t)$ , for Stokes gas flow; $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ , starting height $h_{s}=0.5~\text{mm}$ . The solid red line is the result from the BI simulation. The theoretical result from (3.6) has been superimposed. After a start-up regime for large $h_{0}$ , the deflection ${\it\delta}|_{r=0}$ converges towards a $-1/2$ power law. The BI results agree perfectly with the theoretical predictions, until ${\it\delta}$ and $h_{0}$ become of comparable magnitude, pointed out by the crossing with the solid grey line ${\it\delta}|_{r=0}=h_{0}$ . At that moment ${\it\delta}_{r=0}$ saturates to a constant value, which is the ‘dimple height’ $H_{d}$ of Bouwhuis et al. (Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012).

To illustrate and validate our analysis, we compare the predicted profiles with the results obtained by BI simulations (Oguz & Prosperetti Reference Oguz and Prosperetti1993; Pozrikidis Reference Pozrikidis1997; Bergmann et al. Reference Bergmann, van der Meer, Gekle, van der Bos and Lohse2009). The simulation method is the same as in Bouwhuis et al. (Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012, Reference Bouwhuis, Winkels, Peters, Brunet, van der Meer and Snoeijer2013): the liquid within the pool is described as a potential flow, while the pressure along the pool surface is explicitly calculated from the viscous lubrication equation for the gas flow. To be able to confirm our theoretical predictions in the inertial regime without the influences of surface tension and hydrostatics (which are both very small, as mentioned in the introduction), ${\it\gamma}$ and $g$ are equal to zero in our simulations. In the limit of small deflection, the simulations should thus recover (3.3).

Figure 2(a) shows the configuration on the scale of the sphere, for typical impact parameters for a sphere in air ( $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ ). The interface deflection ${\it\delta}$ is shown in figure 2(b), at the moment when the sphere is at a height $h_{0}=100~{\rm\mu}\text{m}$ . At this time, ${\it\delta}\ll h_{0}\ll R$ , for which we expect agreement between the BI results and our prediction from (3.3). Figure 2(c) shows the acceleration $\partial ^{2}{\it\delta}/\partial t^{2}$ versus $r$ . The solid line is the result from the BI simulations and indeed gives perfect agreement with the prediction, represented by the dots.

The actual deflection profile ${\it\delta}(r,t)$ cannot be integrated explicitly from (3.3), due to the time dependence through $L$ and $u$ . However, we can derive ${\it\delta}|_{r=0}$ , the deflection of the pool surface on the axis, which does not involve $L(t)$ . Using that $\partial /\partial t=-U\partial /\partial h_{0}$ , we find

(3.6) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}|_{r=0}}{\partial h_{0}^{2}}=\frac{3{\it\eta}_{g}G_{1}(0)}{{\it\rho}_{l}UR^{2}}\left(\frac{R}{L}\right)^{5}=\frac{3{\it\eta}_{g}G_{1}(0)}{{\it\rho}_{l}UR^{2}}\left(\frac{R}{h_{0}}\right)^{5/2},\end{eqnarray}$$

where (3.5) implies $G_{1}(0)=(3/8)\sqrt{2}{\rm\pi}$ . The solution of (3.6) for ${\it\delta}|_{r=0}$ is subject to start-up effects as long as $h_{0}\sim h_{s}$ , where $h_{s}$ is the initial height of the gap. If we let the initial height $h_{s}\rightarrow \infty$ , we find

(3.7) $$\begin{eqnarray}{\it\delta}|_{r=0}\simeq \frac{3}{2}\sqrt{2}{\rm\pi}\frac{{\it\eta}_{g}}{{\it\rho}_{l}U}\left(\frac{R}{h_{0}}\right)^{1/2}.\end{eqnarray}$$

This predicts that the central height increases dramatically when $h_{0}$ decreases, as ${\it\delta}\sim h_{0}^{-1/2}$ . Figure 3 shows the BI result for ${\it\delta}|_{r=0}$ against $h_{0}$ (solid line), superimposed with the theoretical predictions (dashed line, taking into account the finite initial height $h_{s}$ ). Indeed, as soon as $h_{0}\ll h_{s}$ , ${\it\delta}|_{r=0}$ converges to a $-1/2$ power law. As expected, the simulation results depart from the analytical prediction when ${\it\delta}\sim h_{0}$ (indicated by the solid grey line) and the lubrication approximation ceases to be valid. At this point, the deflection converges to a constant, which will be the final dimple height $H_{d}$ . As stated in the introduction, this will determine the dimple volume, and thus the entrapped air bubble volume, independently of the air film rupture process.

The current analysis provides a rigorous foundation for the scaling results obtained previously in Hicks & Purvis (Reference Hicks and Purvis2011), Marston et al. (Reference Marston, Vakarelski and Thoroddsen2011), Bouwhuis et al. (Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012), Hicks et al. (Reference Hicks, Ermanyuk, Gavrilov and Purvis2012) and Mandre & Brenner (Reference Mandre and Brenner2012). There, the ‘dimple height’ $H_{d}$ was observed to approach a constant value during the final stages of the impact. Figure 3 shows that this height can be estimated from ${\it\delta}_{r=0}\sim h_{0}\sim H_{d}$ . Using (3.7), this immediately gives

(3.8) $$\begin{eqnarray}H_{d}\sim \frac{{\it\eta}_{g}R^{1/2}}{{\it\rho}_{l}UH_{d}^{1/2}}\sim RSt^{-2/3},\end{eqnarray}$$

where $St={\it\rho}_{l}UR/{\it\eta}_{g}$ is the Stokes number. The corresponding volume of the entrapped bubble volume then scales as

(3.9) $$\begin{eqnarray}V_{b}\sim L^{2}H_{d}\sim RH_{d}^{2}\sim R^{3}St^{-4/3},\end{eqnarray}$$

where we use the common estimate that $L$ sets the lateral scale of the bubble. These are precisely the scaling predictions for the inertial regime (for Stokes gas flow), where the assumptions $H_{d}\sim {\it\delta}$ and $L\sim (H_{d}R)^{1/2}$ were further validated (Hicks & Purvis Reference Hicks and Purvis2011; Marston et al. Reference Marston, Vakarelski and Thoroddsen2011; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Hicks et al. Reference Hicks, Ermanyuk, Gavrilov and Purvis2012; Mandre & Brenner Reference Mandre and Brenner2012; Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013).

3.2. Potential gas flow

As described in § 2.1, the inertial phase of the impacting sphere consists of two distinct stages: the large-gap regime $h_{0}\gg R$ and the thin-gap regime $h_{0}\ll R$ . Below, we separately treat both limiting cases analytically. We furthermore perform a numerical potential flow calculation for the full range of $h_{0}/R$ , to validate the analysis and to show how the two stages are connected.

3.2.1. Large-gap regime: $h_{0}\gg R$

When the sphere is very far from the pool surface, the flow field can be described by the well-known potential flow field around a moving sphere of radius $R$ . The introduction of the (undeformed) pool surface, however, requires that the gas velocity has no vertical component, or $v_{z}|_{z=0}=0$ . This boundary condition can be satisfied using the ‘method of images’, corresponding to two approaching spheres having radius $R$ with approaching velocity $U$ towards a mirroring horizontal line ( $z=0$ ). By applying the superposition of the potentials for the two moving spheres, one obtains the potential

(3.10) $$\begin{eqnarray}{\it\phi}(r,z,t)=\frac{UR^{3}}{2}\left[\frac{\left(z-R-h_{0}\right)}{\left(r^{2}+\left(z-R-h_{0}\right)^{2}\right)^{3/2}}-\frac{z+R+h_{0}}{\left(r^{2}+\left(z+R+h_{0}\right)^{2}\right)^{3/2}}\right].\end{eqnarray}$$

It is important to realize that the introduction of the second moving sphere not only influences the flow around $z=0$ , but also gives a small unwanted velocity on the boundary of the original sphere. In the limit of very large gaps, $R/h_{0}\ll 1$ , this correction becomes 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$ , by applying the unsteady Bernoulli equation:

(3.11) $$\begin{eqnarray}P_{g}(r,t)={\it\rho}_{g}U^{2}\left[\left(\frac{R}{{\it\zeta}}\right)^{3}F_{2}(u)+\frac{9}{2}\left(\frac{R}{{\it\zeta}}\right)^{6}F_{3}(u)\right]\simeq {\it\rho}_{g}U^{2}\left(\frac{R}{{\it\zeta}}\right)^{3}F_{2}(u).\end{eqnarray}$$

Here, ${\it\zeta}(t)=R+h_{0}(t)=R+h_{s}-Ut$ , the radial direction is scaled as $u(t)=r/{\it\zeta}$ , while the spatial profiles are

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle F_{2}(u)=\frac{2-u^{2}}{(1+u^{2})^{5/2}}, & \displaystyle\end{eqnarray}$$
(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle F_{3}(u)=\frac{-u^{2}}{(1+u^{2})^{5}}. & \displaystyle\end{eqnarray}$$
Since (3.10) and (3.11) are only valid for $h_{0}\gg R$ , we only keep the dominant first term in (3.11). It should be noted that the width of the pressure peak is now set by the scale ${\it\zeta}=h_{0}+R$ . This can be contrasted with the width in the thin-gap limit, $L=\sqrt{Rh_{0}}$ , which becomes very narrow.

Figure 4.

Deflection of the pool interface for potential gas flow in the limit $h_{0}\gg R$ ; $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ , $h_{0}=h_{s}=10~\text{mm}$ (thus, $h_{0}/R=10$ ). (a) Global view plot of the sphere and pool contours, (b)  ${\it\delta}$ against $r$ and (c)  $\partial ^{2}{\it\delta}/\partial t^{2}$ against $r$ . The solid red lines result from the BI simulation. The theoretical result from (3.14) has been superimposed in (c) (blue dots). Note the difference in scales on the vertical axes of (a) and (b). The BI results are in good agreement with the theoretical predictions, until $h_{0}/R$ becomes of order 1.

Next, from (3.11) we can compute the induced acceleration profile using (2.7):

(3.14) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}}{\partial t^{2}}(r,t)=\frac{{\it\rho}_{g}U^{2}}{{\it\rho}_{l}{\it\zeta}}\left(\frac{R}{{\it\zeta}}\right)^{3}G_{2}(u).\end{eqnarray}$$

One recognizes a dimensional prefactor that is separated from the spatiotemporal dependence. The function $G_{2}(u)=\int _{0}^{\infty }\widehat{F_{2}}\text{J}_{0}(ku)k^{2}\,\text{d}k$ is the spatial similarity profile, where $\widehat{F_{2}}(k)$ is the Hankel transform of $F_{2}(u)$ . For $G_{2}(u)$ we did not find any analytical expression, but one can numerically calculate the given integral (cf. figure 4).

Once again, we can analytically compute the behaviour of the central deflection, ${\it\delta}|_{r=0}$ :

(3.15) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}|_{r=0}}{\partial h_{0}^{2}}=\frac{{\it\rho}_{g}G_{2}(0)}{{\it\rho}_{l}R}\left(\frac{R}{{\it\zeta}}\right)^{4}.\end{eqnarray}$$

Recalling that $\partial /\partial h_{0}=\partial /\partial {\it\zeta}$ and ${\it\zeta}\rightarrow 2R$ for $h_{0}\rightarrow R$ , this implies that the final ${\it\delta}_{r=0}$ scales as ${\it\rho}_{g}R/{\it\rho}_{l}$ . In contrast to the result for viscous flow, the typical deformation versus $h_{0}$ depends only on the density ratio ${\it\rho}_{g}/{\it\rho}_{l}$ , and not on the impact velocity. While the density ratio is typically small, we anticipate that the resulting deflection for a millimetre-sized sphere can be a few microns. This is actually comparable to typical deflections in the viscous lubrication phase. However, the pool is not deformed locally over a small width $\sqrt{Rh_{0}}$ , but over the scale of the entire sphere, and therefore it will be of little consequence for the formation of the dimple and the size of the entrapped air bubble.

3.2.2. Thin-gap regime: $h_{0}\ll 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-integrated continuity equation (Snoeijer, Brunet & Eggers Reference Snoeijer, Brunet and Eggers2009; Bouwhuis et al. Reference Bouwhuis, Winkels, Peters, Brunet, van der Meer and Snoeijer2013)

(3.16) $$\begin{eqnarray}\frac{\partial h}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left(rh\overline{u}_{r}\right)=0,\end{eqnarray}$$

where $\overline{u}_{r}(r,t)$ is the height-averaged radial gas velocity in the gap. Assuming a plug flow that does not depend on the $z$ coordinate, this average simply gives $\overline{u}_{r}(r,t)=u_{r}(r,t)$ . This analytical description is similar to that of Wilson (Reference Wilson1991), who also studied cushioning air layers at solid–liquid impact in the inertial thin-gap regime, although in 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_{0}(t)+r^{2}/(2R)$ , and thus $\partial h/\partial t=\partial h_{0}/\partial t=-U$ . Hence, we can integrate (3.16) to find

(3.17) $$\begin{eqnarray}u_{r}=\overline{u}_{r}=\frac{Ur}{2h_{0}\left(1+\displaystyle \frac{r^{2}}{2Rh_{0}}\right)}.\end{eqnarray}$$

The velocity profile (3.17) has a local maximum at $r=\sqrt{2Rh_{0}}$ , and vanishes for $r=0$ and $r=\infty$ . Substitution of the profile into the radial component of the Euler equation and integration over $r$ gives the gas pressure:

(3.18) $$\begin{eqnarray}P_{g}(r,t)=\frac{{\it\rho}_{g}U^{2}R}{2h_{0}}\left(\frac{1+\displaystyle \frac{r^{2}}{4Rh_{0}}}{\left(1+\displaystyle \frac{r^{2}}{2Rh_{0}}\right)^{2}}\right)=\frac{{\it\rho}_{g}U^{2}}{2}\left(\frac{R}{L}\right)^{2}F_{4}(u),\end{eqnarray}$$

with $L(t)=\sqrt{Rh_{0}(t)}$ , $u(t)=r/L$ and

(3.19) $$\begin{eqnarray}F_{4}(u)=\frac{1+{\textstyle \frac{1}{4}}u^{2}}{\left(1+{\textstyle \frac{1}{2}}u^{2}\right)^{2}}.\end{eqnarray}$$

It should be noted that the geometry of the thin gap again gives rise to a highly localized pressure profile of width $\sqrt{Rh_{0}}$ . The gas pressure again tends to diverge as $h_{0}\rightarrow 0$ , but more slowly than in the viscous case: the inertial gas pressure in the thin-gap limit is proportional to $1/h_{0}$ , in contrast to the more singular scaling for the viscous gas flow scenario, $1/h_{0}^{2}$ .

From (2.7) we deduce the pool surface acceleration

(3.20) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}}{\partial t^{2}}(r,t)=\frac{{\it\rho}_{g}U^{2}}{2{\it\rho}_{l}L}\left(\frac{R}{L}\right)^{2}G_{4}(u),\end{eqnarray}$$

where $G_{4}(u)=\int _{0}^{\infty }\widehat{F_{4}}\text{J}_{0}(ku)k^{2}\,\text{d}k$ , with $\widehat{F_{4}}(k)$ the Hankel transform of $F_{4}(u)$ . At the origin $r=0$ , this reduces to

(3.21) $$\begin{eqnarray}\frac{\partial ^{2}{\it\delta}|_{r=0}}{\partial h_{0}^{2}}=\frac{{\it\rho}_{g}G_{4}(0)}{2{\it\rho}_{l}R}\left(\frac{R}{L}\right)^{3}.\end{eqnarray}$$

Just like in the case of the large-gap regime, the central deflection has no dependence on the impact velocity. Solving gives ${\it\delta}_{r=0}\sim h_{0}^{1/2}+\text{integration constants}$ . From this we conclude that in the inertial thin-gap limit, the pressure tends to diverge for $h_{0}\rightarrow 0$ , but the deflection ${\it\delta}$ 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 deflection. The integration constants depend on the full history of the impact process, which thus involves the dynamics during the preceding large-gap regime. To predict the actual deflection during the final stages of sphere impact, it is thus not sufficient to consider the large-gap or thin-gap regime of the potential gas flow problem; numerical simulation of the full impact process over all $h_{0}/R$ is required.

3.2.3. Numerical simulations

Simulation of the potential gas flow impact process using the BI technique calls for a different approach with respect to the case of Stokes gas flow. The reason is that we require the gas pressure over the full range of gap thickness, including $h_{0}\sim R$ , for which no analytical solution for the gas pressure is available that can serve as a boundary condition for the liquid pool. As a consequence, the gas phase must also be computed numerically, which we achieve using the BI code. We thus need to run two separate simulations. The process is started by a BI simulation of a solid sphere impacting towards an undeformed surface, with a potential gas flow in between. From this simulation, the gas pressure profile along the pool surface ( $z=0$ ) is extracted. In the second BI simulation, this pressure is applied on a deformable pool surface, from which we eventually determine the resulting pool deflections. This is again a valid method as long as ${\it\delta}/h\ll 1$ . The pressure data are transmitted from the first simulation to the second simulation through an extensive data file. It should be noted that in performing two separate simulations, one needs to take into account the different length scales during the impact process (for $h_{0}=10~\text{mm}\rightarrow 100~\text{nm}$ ), implying very sensitive local node spacings and time dependences. This was achieved by adapting the node spacing and time steps to ensure convergence of the numerical results.

Figure 5.

Inertial gas flow in the thin-gap limit. Theoretical prediction (red dashed line) and BI gas flow simulation result (blue solid line) of the velocity (a) and pressure (b) profiles within the gas; $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ , $h_{s}=h_{0}=100~\text{nm}$ (thus, $h_{0}/R=10^{-4}$ ). We find very good agreement between the theoretical predictions and the BI results.

Figures 4(a) and 4(b) show the configuration on the length scale of the sphere and the interface deflection ${\it\delta}(r,t)$ respectively, for $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ and $h_{0}=h_{s}=10~\text{mm}$ (i.e. the large-gap regime). Figure 4(c) shows the acceleration profile at the corresponding time, and it is observed to agree very well with the asymptotic result of (3.14) (blue dots). The very small difference between the BI result and the theoretical predictions can be explained by the fact that $h_{s}/R=10$ , implying an expected difference of approximately $10\%$ between the theory and the numerical simulations. We remark that the corresponding deformation (figure 4 b) is very small, as we look at the very initial deformations in the start-up regime. The requirement $h_{0}/R\gg 1$ implies a large initial gap height, which, for the parameter values chosen in Fig. 4 to validate the asymptotics, corresponds to non-physically small deflections. The sensitiveness of the very small pool deflection gave rise to switch-off of the smoothing procedure normally used within the simulations (Oguz & Prosperetti Reference Oguz and Prosperetti1993), such that a tiny instability remained visible around the axis, $r=0$ . We confirmed that this instability has a numerical origin and that it does not influence the result on the scale of the deformations. The thin-gap regime is analysed in figure 5. We again find very good agreement between the analytical gas velocity profile (a) and the pressure profile (b) and the BI results (here, $h_{0}=100~\text{nm}$ ).

Figure 6.

Behaviour of the gas pressure on the axis, $P_{r=0}$ , plotted against the gap height, $h_{0}(t)$ , for potential gas flow; $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ , $h_{s}=10~\text{mm}$ . The red dashed line is the theoretical prediction in the regime $h_{0}\gg R$ ; the green dashed line is the theoretical prediction in the regime $h_{0}\ll R$ . The black dashed line points out the crossover $h_{0}=R$ . The BI gas flow simulation result (blue solid line) indeed follows these predicted behaviours in the corresponding regimes, with a crossover at $h_{0}\sim R$ .

Figure 7.

Deflection of the pool interface on the axis, ${\it\delta}_{r=0}$ , plotted against the gap height, $h_{0}(t)$ , for potential gas flow; $R=1~\text{mm}$ , $U=5~\text{m}~\text{s}^{-1}$ , $h_{s}=10~\text{mm}$ . The solid red line is the result from the BI solid sphere on liquid pool simulations (§ 3.2.3). The theoretical result from (3.15) for the large-gap regime has been superimposed (blue dashed line); ${\it\delta}|_{r=0}$ saturates to a constant. The BI results agree perfectly with the large-gap predictions in the regime $h_{0}\gg R$ . In the regime $h_{0}\ll R$ , ${\it\delta}_{r=0}$ deviates from this prediction, but the difference is relatively small. The dashed grey line points out the crossover $h_{0}=R$ ; the solid grey line points out ${\it\delta}|_{r=0}=h_{0}$ .

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_{0}/R\gg 1$ the pressure calculated by BI (blue line) equals $2{\it\rho}_{g}U^{2}(R/{\it\zeta})^{3}$ (red dashed line), and in the limit $h_{0}/R\ll 1$ the pressure equals ${\it\rho}_{g}U^{2}R/(2h_{0})$ (green dashed line). This confirms the validity of the analytical approaches. Finally, we investigate the deflection of the pool that is induced by the numerically obtained gas pressure. Figure 7 shows the deflection at $r=0$ , the inertial (gas) counterpart of figure 3. As expected, ${\it\delta}_{r=0}$ deviates from the 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_{0}\rightarrow 0$ , the influence of the inertial thin-gap limit remains relatively small. For this particular example, it enhances the deflection by less than a factor of $2$ . This is also one of the reasons why we do not show the corresponding theoretical profile for $\partial ^{2}{\it\delta}/\partial t^{2}$ , which in principle could again be directly calculated from the pressure profile. A second reason is that, in the numerical simulations, the very small gap height of $100~\text{nm}$ needs a very high local node density on both the pool surface and the sphere surface; the difference in length scales of $R$ and $h_{0}$ is four decades, which is very challenging. This necessitates very small time steps to be able to calculate a fair second derivative of the deflection profile in time. In addition, 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 more complicated.

The large-gap prediction for the final ${\it\delta}_{r=0}$ is thus satisfactory, and we conclude with the following scaling law for the resulting dimple height $H_{d}$ for the inertial gas scenario as was concluded from (3.15):

(3.22) $$\begin{eqnarray}H_{d}\sim R\frac{{\it\rho}_{g}}{{\it\rho}_{l}}.\end{eqnarray}$$

This dimple height is independent of the impact velocity of the sphere. Since the surface deformation is the sum of the deformations in both the large-gap and the thin-gap limit, it is unclear what the correct radial and axial length scales are that lead to the volume of the pinched bubble.

4. Conclusion

We performed a perturbation analysis to investigate the initial deflections of a liquid surface, induced by the approach of an impacting solid sphere. The analysis assumed that the deflection is limited by the inertia of the liquid pool (i.e. not by its surface tension) and we considered two natural limits for the surrounding medium: Stokes gas flow and potential gas flow. We obtained a quantitative prediction for the pool surface deflection, which was validated numerically, and recovered previously proposed scaling laws for bubble entrapment.

While the ‘cushioning’ of an inertial gas layer had been analysed before (Wilson Reference Wilson1991), most recent work on liquid or solid impact assumes a viscous gas layer. Surprisingly, our analysis reveals that inertial and viscous cushioning both lead to a pool deflection of the order of $1~{\rm\mu}\text{m}$ , for typical experimental conditions. However, the Stokes gas pressure tends to diverge strongly for $h_{0}\rightarrow 0$ , much more strongly than during the inertial gas phase. In addition, this viscous lubrication pressure profile is very localized, while most of the inertial deflection is generated during the initial phase where the pool deflection is spread over the entire width of the sphere. This explains why the experimental results on bubble entrapment are in close agreement with the scaling law (3.9) (Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013), while in addition (3.8) was validated for the case of a liquid drop impact on a solid (Hicks & Purvis Reference Hicks and Purvis2011; Marston et al. Reference Marston, Vakarelski and Thoroddsen2011; Bouwhuis et al. Reference Bouwhuis, van der Veen, Tran, Keij, Winkels, Peters, van der Meer, Sun, Snoeijer and Lohse2012; Hicks et al. Reference Hicks, Ermanyuk, Gavrilov and Purvis2012; Mandre & Brenner Reference Mandre and Brenner2012; Tran et al. Reference Tran, de Maleprade, Sun and Lohse2013); all these results are based on the viscous lubrication regime.

For completeness, we will summarize the possible scenarios for impact of a sphere onto a pool, which can be achieved for different experimental parameters. Assuming an initially high Reynolds number based on the size of the impacting object $R$ , the dynamics will exhibit two different types of crossover: a geometric crossover based on the relative thickness of the gap, $h/R$ , and a crossover from inertial to viscous gas flow. The order in which these crossovers occur depends on the parameters of the problem. In our numerical examples we assumed that one first reaches the thin-gap regime, before the lubrication Reynolds number (based on the gap thickness $h$ ) becomes smaller than unity. This order can be reversed 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 surface tension will become more important, corresponding to the capillary impact regime. In this case, the thin film potentially has time to drain out before a bubble is formed, making the entrapment process more complex (Klaseboer et al. Reference Klaseboer, Chevaillier, Gourdon and Masbernat2000; Yoon et al. Reference Yoon, Borrell, Park and Leal2005).

In this work, we have elaborated on the impact of a solid sphere on a liquid surface. Similar perturbation analysis can be performed for drop impact on a solid, or drop impact on a pool, although the details will be different. This explains why the same scaling laws are observed in all these cases.

Acknowledgement

We gratefully acknowledge S. Wilson and H. Gelderblom for insightful discussion. This work was supported by STW and NWO through a VIDI grant no. 11304.

References

Bergmann, R. P. H. M., van der Meer, D., Gekle, S., van der Bos, J. & Lohse, D. 2009 Controlled impact of a disk on a water surface: cavity dynamics. J. Fluid Mech. 633, 381–409.Google Scholar
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. 2012 Maximal air bubble entrainment at liquid-drop impact. Phys. Rev. Lett. 109, 264501.Google Scholar
Bouwhuis, W., Winkels, K. G., Peters, I. R., Brunet, P., van der Meer, D. & Snoeijer, J. H. 2013 Oscillating and star-shaped drops levitated by an airflow. Phys. Rev. E 88, 023017.Google Scholar
van Dam, D. & Le Clerc, C. 2004 Experimental study of an ink-jet printed droplet on a solid substrate. Phys. Fluids 16, 3403–3414.CrossRefGoogle Scholar
Davis, R. H., Serayssol, J.-M. & Hinch, E. J. 1986 The elastohydrodynamic collision of two spheres. J. Fluid Mech. 163, 479–497.CrossRefGoogle Scholar
Deng, Q., Anilkumar, A. V. & Wang, T. G. 2009 The phenomenon of bubble entrapment during capsule formation. J. Colloid Interface Sci. 333 (2), 523–532.CrossRefGoogle ScholarPubMed
Do-Quang, M. & Amberg, G. 2009 The splash of a solid sphere impacting on a liquid surface: numerical simulation of the influence of wetting. Phys. Fluids 21 (2), 022102.Google Scholar
Driscoll, M. M. & Nagel, S. R. 2011 Ultrafast interference imaging of air in splashing dynamics. Phys. Rev. Lett. 107, 154502.CrossRefGoogle ScholarPubMed
Hicks, P. D., Ermanyuk, E. V., Gavrilov, N. V. & Purvis, R. 2012 Air trapping at impact of a rigid sphere onto a liquid. J. Fluid Mech. 695, 310–320.Google Scholar
Hicks, P. D. & Purvis, R. 2011 Air cushioning in droplet impacts with liquid layers and other droplets. Phys. Fluids 23, 062104.Google Scholar
Howison, S. D., Ockendon, J. R. & Wilson, S. K. 1991 Incompressible water-entry problems at small deadrise angles. J. Fluid Mech. 222, 215–230.Google Scholar
Klaseboer, E., Chevaillier, J. P., Gourdon, C. & Masbernat, O. 2000 Film drainage between colliding drops at constant approach velocity: experiments and modeling. J. Colloid Interface Sci. 229 (1), 274–285.Google Scholar
Klaseboer, E., Manica, R. & Chan, D. Y. C. 2014 Universal behavior of the initial stage of drop impact. Phys. Rev. Lett. 113, 194501.Google Scholar
Korobkin, A. A., Ellis, A. S. & Smith, F. T. 2008 Trapping of air in impact between a body and shallow water. J. Fluid Mech. 611, 365–394.Google Scholar
Korobkin, A. A. & Pukhnachov, V. V. 1988 Initial stage of water impact. Annu. Rev. Fluid Mech. 20, 159–185.Google Scholar
Lamb, H. 1957 Hydrodynamics, 6th edn. Cambridge University Press.Google Scholar
Mandre, S. & Brenner, M. P. 2012 The mechanism of a splash on a dry solid surface. J. Fluid Mech. 690, 148–172.Google Scholar
Marston, J. O., Vakarelski, I. U. & Thoroddsen, S. T. 2011 Bubble entrapment during sphere impact onto quiescent liquid surfaces. J. Fluid Mech. 680, 660–670.CrossRefGoogle Scholar
Moore, M. R. & Oliver, J. M. 2014 On air cushioning in axisymmetric impacts. IMA J. Appl. Maths 79, 661–680.Google Scholar
Oguz, H. N. & Prosperetti, A. 1993 Dynamics of bubble growth and detachment from a needle. J. Fluid Mech. 257, 111–145.CrossRefGoogle Scholar
Pozrikidis, C. 1997 Introduction to Theoretical and Computational Fluid Dynamics, 1st edn. Oxford University Press.Google Scholar
Reynolds, O. 1886 On the theory of lubrication and its application to Mr. Beauchamp Tower’s experiments, including an experimental determination of the viscosity of olive oil. Phil. Trans. R. Soc. Lond. 177, 157–234.Google Scholar
Smith, F. T., Li, L. & Wu, G. X. 2003 Air cushioning with a lubrication/inviscid balance. J. Fluid Mech. 482, 291–318.Google Scholar
Snoeijer, J. H., Brunet, P. & Eggers, J. 2009 Maximum size of drops levitated by an air cushion. Phys. Rev. E 79, 036307.CrossRefGoogle ScholarPubMed
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. J. Fluid Mech. 545, 203–212.Google Scholar
Thoroddsen, S. T., Thoraval, M. J., Takehara, K. & Etoh, T. G. 2012 Micro-bubble morphologies following drop impacts onto a pool surface. J.  Fluid Mech. 708, 469–479.Google Scholar
Tran, T., de Maleprade, H., Sun, C. & Lohse, D. 2013 Air entrainment during impact of droplets on liquid surfaces. J. Fluid Mech. 726, R3.Google Scholar
Wilson, S. K. 1991 A mathematical model for the initial stages of fluid impact in the presence of a cushioning fluid layer. J. Engng Maths 25, 265–285.CrossRefGoogle Scholar
Yiantsios, S. G. & Davis, R. H. 1990 On the buoyancy-driven motion of a drop towards a rigid surface or a deformable interface. J. Fluid Mech. 217, 547–573.Google Scholar
Yoon, Y., Borrell, M., Park, C. C. & Leal, L. G. 2005 Viscosity ratio effects on the coalescence of two equal-sized drops in a two-dimensional linear flow. J. Fluid Mech. 525, 355–379.CrossRefGoogle Scholar
Figure 0

Figure 1.(a) A solid sphere (radius $R$) approaches a liquid surface with velocity $U$. The gap height between the bottom of the sphere and the undisturbed water level ($z=0$) is $h(r,t)$, where $r$ and $t$ are the radial coordinate and time respectively, with $h(0,t)=h_{0}(t)$. (b) While the sphere moves downwards, the pool deflects by a small amount ${\it\delta}(r,t)$, as a result of the local pressure buildup in the air that is squeezed out. In the limit where ${\it\delta}\ll h_{0}$, which typically is valid up to very close to the impact time, the profiles are computed analytically. (c) This mechanism will finally result in air bubble entrapment.

Figure 1

Figure 2.Deflection of the pool interface for Stokes gas flow; $R=1~\text{mm}$, $U=5~\text{m}~\text{s}^{-1}$, starting height of the (bottom of the) sphere $h_{s}=0.5~\text{mm}$, current height $h_{0}=0.1~\text{mm}$. (a) Global view of the sphere and pool contours, (b) the pool deflection ${\it\delta}$ as a function of $r$ and (c) $\partial ^{2}{\it\delta}/\partial t^{2}$ as a function of $r$. The solid red lines result from the BI simulation. The theoretical result from (3.3) has been superimposed in (c) (blue dots). Note the difference in scales on the vertical axes of (a) and (b). The BI results agree perfectly with the theoretical predictions, as long as $|{\it\delta}|\ll h$.

Figure 2

Figure 3.Deflection of the pool interface on the axis, ${\it\delta}_{r=0}$, plotted against $h_{0}(t)$, for Stokes gas flow; $R=1~\text{mm}$, $U=5~\text{m}~\text{s}^{-1}$, starting height $h_{s}=0.5~\text{mm}$. The solid red line is the result from the BI simulation. The theoretical result from (3.6) has been superimposed. After a start-up regime for large $h_{0}$, the deflection ${\it\delta}|_{r=0}$ converges towards a $-1/2$ power law. The BI results agree perfectly with the theoretical predictions, until ${\it\delta}$ and $h_{0}$ become of comparable magnitude, pointed out by the crossing with the solid grey line ${\it\delta}|_{r=0}=h_{0}$. At that moment ${\it\delta}_{r=0}$ saturates to a constant value, which is the ‘dimple height’ $H_{d}$ of Bouwhuis et al. (2012).

Figure 3

Figure 4.Deflection of the pool interface for potential gas flow in the limit $h_{0}\gg R$; $R=1~\text{mm}$, $U=5~\text{m}~\text{s}^{-1}$, $h_{0}=h_{s}=10~\text{mm}$ (thus, $h_{0}/R=10$). (a) Global view plot of the sphere and pool contours, (b) ${\it\delta}$ against $r$ and (c) $\partial ^{2}{\it\delta}/\partial t^{2}$ against $r$. The solid red lines result from the BI simulation. The theoretical result from (3.14) has been superimposed in (c) (blue dots). Note the difference in scales on the vertical axes of (a) and (b). The BI results are in good agreement with the theoretical predictions, until $h_{0}/R$ becomes of order 1.

Figure 4

Figure 5.Inertial gas flow in the thin-gap limit. Theoretical prediction (red dashed line) and BI gas flow simulation result (blue solid line) of the velocity (a) and pressure (b) profiles within the gas; $R=1~\text{mm}$, $U=5~\text{m}~\text{s}^{-1}$, $h_{s}=h_{0}=100~\text{nm}$ (thus, $h_{0}/R=10^{-4}$). We find very good agreement between the theoretical predictions and the BI results.

Figure 5

Figure 6.Behaviour of the gas pressure on the axis, $P_{r=0}$, plotted against the gap height, $h_{0}(t)$, for potential gas flow; $R=1~\text{mm}$, $U=5~\text{m}~\text{s}^{-1}$, $h_{s}=10~\text{mm}$. The red dashed line is the theoretical prediction in the regime $h_{0}\gg R$; the green dashed line is the theoretical prediction in the regime $h_{0}\ll R$. The black dashed line points out the crossover $h_{0}=R$. The BI gas flow simulation result (blue solid line) indeed follows these predicted behaviours in the corresponding regimes, with a crossover at $h_{0}\sim R$.

Figure 6

Figure 7.Deflection of the pool interface on the axis, ${\it\delta}_{r=0}$, plotted against the gap height, $h_{0}(t)$, for potential gas flow; $R=1~\text{mm}$, $U=5~\text{m}~\text{s}^{-1}$, $h_{s}=10~\text{mm}$. The solid red line is the result from the BI solid sphere on liquid pool simulations (§ 3.2.3). The theoretical result from (3.15) for the large-gap regime has been superimposed (blue dashed line); ${\it\delta}|_{r=0}$ saturates to a constant. The BI results agree perfectly with the large-gap predictions in the regime $h_{0}\gg R$. In the regime $h_{0}\ll R$, ${\it\delta}_{r=0}$ deviates from this prediction, but the difference is relatively small. The dashed grey line points out the crossover $h_{0}=R$; the solid grey line points out ${\it\delta}|_{r=0}=h_{0}$.