2.1. Dimensional analysis and numerical method

Since the deformation of the interfaces responsible for air entrapment exhibits axial symmetry, we make use of an axisymmetric framework throughout the paper. The Reynolds number of the liquid drops we model, which is defined as $Re_{l}\equiv {\it\rho}_{l}RU/{\it\eta}_{l}$ , is assumed to be large, $Re_{l}\gg 1$ . Here, ${\it\rho}_{l}$ and ${\it\eta}_{l}$ are respectively the density and the dynamic viscosity of the liquid, $U$ is the impact velocity and $R$ is the radius of the drop. The flow can thus be regarded as irrotational, that is $\boldsymbol{{\rm\nabla}}\times \boldsymbol{u}=\mathbf{0}$ . Under the additional constraint of incompressible flow inside the drop this allows the liquid dynamics to be modelled with a harmonic function ${\it\phi}$ , to which the velocity field $\boldsymbol{u}$ is related through

(2.1) $$\begin{eqnarray}\boldsymbol{u}=\boldsymbol{{\rm\nabla}}{\it\phi}.\end{eqnarray}$$

The fact that the velocity potential ${\it\phi}$ obeys the Laplace equation $\boldsymbol{{\rm\nabla}}^{2}{\it\phi}=0$ is used to efficiently solve the potential problem, and thus the dynamics of the liquid, using the BIM. We use a BIM based on codes that are described in detail in Oguz & Prosperetti (Reference Oguz and Prosperetti1993) and Gekle & Gordillo (Reference Gekle and Gordillo2011). While the Reynolds number of the drop is large, the Reynolds number of the thin gaseous air layer, $Re_{g}\equiv {\it\rho}_{g}H_{d}U/{\it\eta}_{g}$ , is typically small. Here, ${\it\rho}_{g}$ is the gas density and $H_{d}$ is the air film thickness in the centre of the film which is referred to as the dimple height. When inserting typical parameters, ${\it\rho}_{g}$ is of order $1$ , $H_{d}$ is of order $10^{-6}$ , $U$ is of order $1$ and ${\it\eta}_{g}$ is of order $10^{-5}$ . This results in $Re_{g}$ being of the order of $10^{-1}$ , which justifies the assumption. The length scale characterizing the air layer in the lateral extension of the air film is denoted by $L$ , see figure 2(a). As shown in Bouwhuis et al. (Reference Bouwhuis, Van Der Veen, Tran, Keij, Winkels, Peters, Van Der Meer, Sun, Snoeijer and Lohse2012), $H_{d}\ll L$ , which in combination with the low Reynolds number of the gas allows the film to be described with viscous lubrication theory, see for example Leal (Reference Leal1992). For droplets impacting with a higher speed well outside the parameter range currently considered, lubrication theory should be extended to include the effect of inertia of the gas. The dimensionless group reflecting the presence of air is the Stokes number $\mathit{St}\equiv {\it\rho}_{l}RU/{\it\eta}_{g}$ , which compares the viscous force of the air layer with the inertial force in the drop. This number is relevant for describing dimple formation, since, for high enough impact velocity $U$ , this process is determined by two competing forces: the force of the viscous air layer trying to deform the drop in the centre and the opposing inertial force of the drop, which must be slowed down locally in order to form a dimple. Additional dimensionless numbers incorporating the surface tension ${\it\gamma}$ are the Weber number $\mathit{We}$ and the capillary number $\mathit{Ca}$ based on the gas properties. Summarizing, we thus have the following dimensionless parameters:

(2.2a−e ) $$\begin{eqnarray}Re_{l}\equiv \frac{{\it\rho}_{l}RU}{{\it\eta}_{l}},\quad Re_{g}\equiv \frac{{\it\rho}_{g}H_{d}U}{{\it\eta}_{g}},\quad St\equiv \frac{{\it\rho}_{l}RU}{{\it\eta}_{g}},\quad We\equiv \frac{{\it\rho}_{l}RU^{2}}{{\it\gamma}},\quad Ca\equiv \frac{We}{St}.\end{eqnarray}$$

Figure 2.

(a) Schematic of drop impact onto a pool. The methods used are indicated in the figure: both the liquid domains are modelled with potential flow, while the air layer is described with Stokes flow. The grey arrows indicate that the flow of the air film is coupled to the dynamics of the liquid domains, and vice versa. (b) Definition of the ( $n,s$ )-coordinate system, where $s$ is aligned along the drop curve and $n$ is the unit normal with respect to the drop.

The impact of a liquid drop onto a pool of the same liquid and the impact of a rigid sphere onto a liquid pool can be described with the same dimensionless numbers. As the initial geometries of the problems are identical, the difference lies in the deformability of the object, which is zero in the case of the solid. The two effective control parameters that we will use here in our theoretical framework are $\mathit{St}$ and $\mathit{We}$ . In this work, the depth of the pool is considered to be infinite. In the case where the thickness of the pool is finite, the dynamics of the pool may be altered due to the presence of a solid boundary at the bottom of the pool. For the impact onto a liquid film with finite thickness we refer the reader to Korobkin, Ellis & Smith (Reference Korobkin, Ellis and Smith2008), in which the impact onto a wetted solid is discussed.

In figure 2(a), an illustration of the impact of a drop onto a pool, together with the method used, is shown. As is clear from this figure, the coupling between the dynamics of the air layer and the dynamics of the liquid is essential since the two liquid domains feel each other through the pressure build-up in the viscous air layer. The lubrication pressure $P_{g}$ acts on the liquid surface and appears in the unsteady Bernoulli equation which serves as a dynamic boundary condition in the BIM applied at the liquid surface, see also Bouwhuis et al. (Reference Bouwhuis, Van Der Veen, Tran, Keij, Winkels, Peters, Van Der Meer, Sun, Snoeijer and Lohse2012):

(2.3) $$\begin{eqnarray}\left(\frac{\partial {\it\phi}}{\partial t}+\frac{1}{2}\left|\boldsymbol{{\rm\nabla}}{\it\phi}\right|^{2}\right)=-\frac{{\it\gamma}}{{\it\rho}_{l}}{\it\kappa}(s,t)-\frac{P_{g}(s,t)}{{\it\rho}_{l}}.\end{eqnarray}$$

Here, $P_{g}$ is the pressure in excess of the ambient pressure, due to lubrication. The curvature of the interface is represented by ${\it\kappa}(s,t)$ , which is a function of the curvilinear coordinate $s$ , which follows the liquid surface and time $t$ . It should be noted that unlike Bouwhuis et al. (Reference Bouwhuis, Van Der Veen, Tran, Keij, Winkels, Peters, Van Der Meer, Sun, Snoeijer and Lohse2012) we do not include gravitation, to make sure that the impact speed of the impacting droplet stays constant. The small deformation of the pool justifies the assumption of neglecting the pressure due to hydrostatic gradients. As we have two liquid domains, two separate boundary integral equations are solved. We take the width of the pool large enough to approach the dynamics of an infinite liquid pool. In this case, a width of 4.5 times the drop radius is found to be sufficient. We focus on quantifying the amount of entrapped air by integrating the enclosed air pocket up to the moment the air layer reaches a physical minimum thickness of $0.4~{\rm\mu}\text{m}$ . At this point the volume of the enclosed air has converged and a subsequent rupture of the air film will prevent further drainage, which results in an entrapped air bubble (Bouwhuis et al. Reference Bouwhuis, Van Der Veen, Tran, Keij, Winkels, Peters, Van Der Meer, Sun, Snoeijer and Lohse2012). As we focus on the dynamics just prior to rupture we can make use of an axisymmetric framework. In Hicks et al. (Reference Hicks, Ermanyuk, Gavrilov and Purvis2012), a similar approach was used to predict the radius of the entrapped bubble that occurs when a solid sphere approaches a liquid free surface. We restrict ourselves to the inertial regime (Bouwhuis et al. Reference Bouwhuis, Van Der Veen, Tran, Keij, Winkels, Peters, Van Der Meer, Sun, Snoeijer and Lohse2012) for which experimental results (Tran et al. Reference Tran, De Maleprade, Sun and Lohse2013) are available for a direct comparison. Since the air layer continually deforms and translates during the impact, lubrication equations have been developed in a moving coordinate system that is aligned with the interface of the drop. These equations will be derived in the next section.

2.2. Lubrication in a moving and tilted coordinate system

In this section we develop an expression for the pressure $P_{g}$ in the air film based on lubrication theory in a moving ( $n$ , $s$ )-coordinate system that is aligned along the drop surface; see the sketch in figure 2. The reason for doing this (rather than just using the standard ( $r,z$ )-coordinate system) is that especially for the drop onto pool impact, the moving ( $n,s$ )-coordinate system is not necessarily oriented like the ( $r,z$ )-coordinate system and therefore only the first guarantees an accurate description of the draining air film. In appendix A a case is described that shows the difference between lubrication calculated in the two coordinate systems. The drop surface is taken as a reference, and the curvilinear coordinate $s$ is defined along the drop, starting at the axis of symmetry. At some large radial coordinate $s_{\infty }$ we assume atmospheric pressure. The coordinate perpendicular to $s$ is defined to be $n$ . The gap thickness $h(r,t)$ is defined as the length of the perpendicular line from the drop projected onto the liquid pool. The two surfaces in the impact zone are assumed to be nearly parallel ( $\left|\partial _{s}h\right|\ll 1$ ), so we can apply lubrication theory.

It can be shown (see appendix B) that the continuity equation in this new ( $n$ , $s$ )-coordinate system reads

(2.4) $$\begin{eqnarray}\frac{u_{r}}{r}+\partial _{s}u_{s}+\partial _{n}u_{n}=0.\end{eqnarray}$$

At the interface of the liquid pool ( $n=h$ ) we know that the fluid particles have to move with the interface. This is mathematically described with the kinematic boundary condition

(2.5) $$\begin{eqnarray}\partial _{t}h+(u_{s}\partial _{s}h)|_{n=h}=u_{n\mid n=h}-u_{n\mid n=0}.\end{eqnarray}$$

Here, $\partial _{t}h$ is the time derivative of $h$ . We now integrate (2.4) along the gap thickness $h$ and obtain

(2.6) $$\begin{eqnarray}\int _{0}^{h}\frac{u_{r}}{r}\,\text{d}n+\int _{0}^{h}\partial _{s}u_{s}\,\text{d}n=-\int _{0}^{h}\partial _{n}u_{n}\,\text{d}n=u_{n\mid n=0}-u_{n\mid n=h}.\end{eqnarray}$$

Using the Leibniz integral rule for the second integral on the left-hand side we find

(2.7) $$\begin{eqnarray}\int _{0}^{h}\frac{u_{r}}{r}\,\text{d}n+\partial _{s}\int _{0}^{h}u_{s}\,\text{d}n-\left(u_{s}\left.\partial _{s}h\right)\right|_{n=h}=u_{n\mid n=0}-u_{n\mid n=h}.\end{eqnarray}$$

We now use the kinematic boundary condition formulated in (2.5) for the third term on the left-hand side to obtain

(2.8) $$\begin{eqnarray}\int _{0}^{h}\frac{u_{r}}{r}\,\text{d}n+\partial _{s}\int _{0}^{h}u_{s}\,\text{d}n+\partial _{t}h-u_{n\mid n=h}+u_{n\mid n=0}=u_{n\mid n=0}-u_{n\mid n=h}.\end{eqnarray}$$

Cancelling the terms $u_{n\mid n=h}$ and $u_{n\mid n=0}$ on both sides gives

(2.9) $$\begin{eqnarray}\int _{0}^{h}\frac{u_{r}}{r}\,\text{d}n+\partial _{s}\int _{0}^{h}u_{s}\,\text{d}n+\partial _{t}h=0.\end{eqnarray}$$

We still have to describe $u_{r}$ within the new ( $s$ , $n$ )-coordinate system. Therefore, we substitute $u_{r}=u_{n}\cos {\it\theta}-u_{s}\sin {\it\theta}$ in the equation above to obtain

(2.10) $$\begin{eqnarray}\int _{0}^{h}\frac{1}{r}u_{n}\cos {\it\theta}\,\text{d}n-\int _{0}^{h}\frac{1}{r}u_{s}\sin {\it\theta}\,\text{d}n+\partial _{s}\int _{0}^{h}u_{s}\,\text{d}n+\partial _{t}h=0.\end{eqnarray}$$

We assume that the main flow of the air that is squeezed out from the gap is along the $s$ coordinate, which implies that $u_{n}$ is relatively small, so we neglect the first term. The second term is an integral with respect to $n$ containing the variable $r$ . This radial coordinate $r$ across $h$ is a function of $n$ : $r=n\cos {\it\theta}+c(s)$ . Here, $c(s)$ is the value of $r$ at the drop surface ( $n=0$ ) for some coordinate $s$ . We thus substitute this expression for $r$ into (2.10) and neglect the first term to find

(2.11) $$\begin{eqnarray}-\int _{0}^{h}\frac{\sin {\it\theta}}{n\cos {\it\theta}+c(s)}u_{s}\,\text{d}n+\partial _{s}\int _{0}^{h}u_{s}\,\text{d}n+\partial _{t}h=0.\end{eqnarray}$$

2.2.1. Flow profile within the air film

As has been previously described, the Reynolds number of the thin air film is small, $Re_{g}\ll 1$ , and the geometry of the problem, $H_{d}\ll L$ , allows us to use lubrication theory. In the ( $n,s$ )-coordinate system, the Stokes equations can then be reduced to

(2.12) $$\begin{eqnarray}\partial _{s}P_{g}={\it\eta}_{g}\partial _{n}^{2}u_{s}.\end{eqnarray}$$

We can integrate (2.12) twice with respect to $n$ , employing a no-slip boundary condition at the drop surface ( $u_{s}=U_{d}$ at $n=0$ ) as well as at the surface of the pool ( $u_{s}=U_{p}$ at $n=h$ ):

(2.13) $$\begin{eqnarray}u_{s}=\left((U_{p}-U_{d})\frac{n}{h}+U_{d}\right)+\frac{1}{2{\it\eta}_{g}}\partial _{s}P_{g}(n^{2}-nh).\end{eqnarray}$$

The first term of (2.13) can be associated with Couette flow, caused by the movement of the interfaces. The terms $U_{p}$ and $U_{d}$ can be extracted by differentiating the potential ${\it\phi}$ with respect to the tangential direction of the liquid surface. The second term can be associated with Poiseuille flow, which is driven by the pressure gradient, see also Klaseboer et al. (Reference Klaseboer, Chevaillier, Gourdon and Masbernat2000). Substituting this expression for $u_{s}$ in our equation for mass conservation, (2.11), we obtain

(2.14) $$\begin{eqnarray}\displaystyle & & \displaystyle -\int _{0}^{h}\frac{\sin {\it\theta}}{n\cos {\it\theta}+c}\left[\left((U_{p}-U_{d})\frac{n}{h}+U_{d}\right)+\frac{1}{2{\it\eta}_{g}}\partial _{s}(n^{2}-nh)\right]\text{d}n\nonumber\\ \displaystyle & & \displaystyle \quad +\,\partial _{s}\int _{0}^{h}\left[\left((U_{p}-U_{d})\frac{n}{h}+U_{d}\right)+\frac{1}{2{\it\eta}_{g}}\partial _{s}P_{g}(n^{2}-nh)\right]\text{d}n+\partial _{t}h=0.\end{eqnarray}$$

In the first integral we deal with a prefactor $\sin {\it\theta}/(n\cos {\it\theta}+c)$ . When taking the geometry of the problem into account we note that $n\cos {\it\theta}\ll c$ . We can thus write $\sin {\it\theta}/(n\cos {\it\theta}+c)\approx \sin {\it\theta}/c$ . Performing the integrals of (2.14) under this assumption yields

(2.15) $$\begin{eqnarray}\displaystyle -\frac{\sin {\it\theta}}{c}\left(\frac{h}{2}\left(U_{p}+U_{d}\right)-\frac{h^{3}}{12{\it\eta}_{g}}\partial _{s}P_{g}\right)+\partial _{s}\left(\frac{h}{2}\left(U_{p}+U_{d}\right)-\frac{h^{3}}{12{\it\eta}_{g}}\partial _{s}P_{g}\right)+\partial _{t}h=0.\qquad & & \displaystyle\end{eqnarray}$$

If we define $G(s)\equiv (h/2(U_{p}+U_{d})-(h^{3}/12{\it\eta}_{g})\partial _{s}P_{g})$ we can transform the above equation into a first-order inhomogeneous linear ordinary differential equation (ODE) for $G(s)$ :

(2.16) $$\begin{eqnarray}{\dot{G}}(s)-a(s)G(s)=f(s).\end{eqnarray}$$

Here, $a(s)$ and $f(s)$ are known functions of $s$ :

(2.17) $$\begin{eqnarray}a(s)=\frac{\sin {\it\theta}}{c(s)},\end{eqnarray}$$

(2.18) $$\begin{eqnarray}f(s)=-\partial _{t}h.\end{eqnarray}$$

2.2.2. Solving the first-order inhomogeneous ODE for $G(s)$

Equation (2.16) can be solved with the help of an integrating factor $I$ , defined as $I(s)\equiv \text{e}^{-\int a(s)\,\text{d}s}$ . Using the boundary condition $G(s)=0$ for $s=0$ , because we have zero pressure gradient in the centre of symmetry, and also zero tangential velocities, we can multiply (2.16) with $I(s)$ and solve for $G(s)$ :

(2.19) $$\begin{eqnarray}G(s)=\frac{1}{I(s)}\left(\int _{0}^{s}I(\tilde{s})f(\tilde{s})\,\text{d}\tilde{s}\right),\end{eqnarray}$$

with $I(s)=\text{e}^{-\int _{0}^{s}a(\tilde{s})\,\text{d}\tilde{s}}$ .

We can now substitute $G(s)\equiv (h/2(U_{p}+U_{d})-(h^{3}/12{\it\eta}_{g})\partial _{s}P_{g})$ back into (2.19) to obtain an equation for $\partial _{s}P_{g}$ :

(2.20) $$\begin{eqnarray}\partial _{s}P_{g}=-\frac{12{\it\eta}_{g}}{h^{3}}\left(\frac{1}{I(s)}\left(\int _{0}^{s}I(\tilde{s})f(\tilde{s})\,\text{d}\tilde{s}\right)-\frac{h}{2}\left(U_{p}+U_{d}\right)\right).\end{eqnarray}$$

We note that we have to evaluate two numerical integrals to calculate $\partial _{s}P_{g}$ . In order to find the pressure $P_{g}(s)$ , we integrate (2.20) using atmospheric pressure for some large value for $s_{\infty }$ well outside the thin air gap as a boundary value. As a check of our analysis we now orientate the ( $n$ , $s$ )-coordinate system in such a way that $s=r$ , to recover the lubrication equation in the conventional ( $r,z$ )-coordinate system. In that case, we have ${\it\theta}=-{\rm\pi}/2$ , and we can write for $a(s)$

(2.21) $$\begin{eqnarray}a(s=r)=\frac{\sin {\it\theta}}{r}=-\frac{1}{r}.\end{eqnarray}$$

The integrating factor $I$ now becomes

(2.22) $$\begin{eqnarray}I(s=r)=\text{e}^{-\int _{0}^{r}a(\tilde{r})\,\text{d}\tilde{r}}=\text{e}^{\ln r}=r.\end{eqnarray}$$

Substituting (2.22) into (2.20), using the proposition $s=r$ and setting $U_{b}=0$ and $U_{d}=0$ , we can now write (2.20) as

(2.23) $$\begin{eqnarray}\partial _{r}P_{g}=-\frac{12{\it\eta}_{g}}{h^{3}}\left(\frac{1}{I(r)}\left(\int _{0}^{r}I(\tilde{r})f(\tilde{r})\,\text{d}\tilde{r}\right)-\frac{h}{2}\left(U_{b}+U_{d}\right)\right)=\frac{12{\it\eta}_{g}}{h^{3}}\left(\frac{1}{r}\left(\int _{0}^{r}\tilde{r}\partial _{t}h\,\text{d}\tilde{r}\right)\right).\end{eqnarray}$$

Figure 3.

(a–h) Drop impact onto a liquid pool: (a,e)  $t=0$ , (b,f)  $t=0.13$ , (c,g)  $t=0.15$ , (d,h)  $t=0.17$ . Note the different length scales for the $r$ -axis and $z$ -axis in the shape plots. The impact speed is $U=0.42~\text{m}~\text{s}^{-1}$ and the drop radius is $R=0.95~\text{mm}$ . The density and surface tension of the liquid are respectively ${\it\rho}=916~\text{kg}~\text{m}^{-3}$ and ${\it\gamma}=0.020~\text{N}~\text{m}^{-1}$ . These impact parameters correspond to $\mathit{St}=2.0\times 10^{4}$ and $\mathit{We}=7.7$ . The simulation starts at time $t=0\;\text{ms}$ at a separation of $h_{r=0}=50~{\rm\mu}\text{m}$ . Due to the approach of the drop, the excess air pressure $P_{g}$ will increase and acts on both the drop and the liquid pool ( $t=0.13\;\text{ms}$ ). At the final stage ( $t=0.17\;\text{ms}$ ) the minimum separation of the interfaces reaches $0.4~{\rm\mu}\text{m}$ and the simulation is stopped. The bubble volume $V_{b}$ can thus be determined. (i–k) Part of the simulation domain with detailed snapshots of the air film at $t=0.17\;\text{ms}$ . In (k), the actual node distribution around the smallest separation point can be inspected. This is the most refined distribution of computational nodes that is used. For the region outside the gap a coarser node distribution is sufficient.

We check by inspection that this (2.23) is the equation for the radial pressure gradient for viscous lubrication theory in the conventional ( $r,z$ )-coordinate system (Bouwhuis et al. Reference Bouwhuis, Winkels, Peters, Brunet, Van Der Meer and Snoeijer2013), which gives a consistency check for our analysis. This was also numerically verified.