3.1. Global overview

We discover that the mechanism behind air crater formation is (figure 1c) intricately related to the penetration dynamics of the impinging droplet through the liquid pool and depends on the impact Weber number. Figure 1(c) depicts the schematic representation of the surface crater formation on the air–water interface of the impacting droplet. Figure 2 shows the high-speed snapshots of a water droplet impacting a silicone liquid pool at various impact Weber numbers ($We=4,16,65,145$). The timestamps are in milliseconds, and the scale bars represent 1 mm. The interaction between the silicone pool and water droplet occurs through the interplay of viscous, buoyancy, droplet weight, air–water and water–oil surface tension forces. Beyond a critical impact Weber number ($We_c\sim 10$), air craters of significant depths (crater depth comparable to droplet size) are formed. No significant air craters are observed during the penetration phase for impact Weber number smaller than the critical Weber number ($We< We_c$). The penetration time scale is largely increased due to the lubrication effect of the entrapped air layer formed under the droplet. The air layer rupture occurs at a time scale of the order of $t_r\sim 2\times 10^{-1}\ {\rm s}$. The transient dynamics of the droplet (surface oscillation) occurs early above the pool for $t< t_r$. For $t>t_r$, the droplet penetrates smoothly without any surface oscillation. The deceleration of the droplet is delayed due to the entrapped air and is non-impulsive. The retarding force during the penetration phase is of the order of $F\sim {O}(4\times 10^{-4})~{\rm N}$ for $We=4$. The air craters formed for $We>We_c$ occur due to the rapid deceleration experienced by the impacting droplet over a short period (impulsive deceleration). For impact Weber number larger than the critical Weber number ($We>We_c$), the air layer rupture time scale is of the order of $t_r\sim {O}(5\times 10^{-4})\ {\rm s}$, approximately two orders of magnitude faster than for low Weber number ($We\sim 1$). We have shown using scaling analysis that the force responsible for the sudden deceleration is the viscous drag. The viscous drag force is of the order of $F_v\sim {O}(1.5\times 10^{-3}\ {\rm N})$ for $We=16$ which is comparatively larger than the buoyancy/weight of the droplet (${O}(5\times 10^{-5}\ {\rm N})$). Further, we deduce a criterion for significant air crater depths. We show that the ratio of Stokes viscous drag pressure to capillary pressure (${\varGamma }=\Delta p_s/\Delta p_c$) is the parameter characterizing the formation and the size of air craters; $\Delta p_s/\Delta p_c\sim 1$ corresponds to low Weber number states ($We=4$). The ratio ${\varGamma }$ increases as we increase the impact Weber number. The shape of the air crater central jet depends on the response of the droplet subjected to impulsive retardation. We show that the shape of the central air jet can be approximated with a fourth-order Legendre polynomial. We observe that the air crater formation/response time scales for all Weber numbers greater than the critical lie within the capillary time scale of $t_{c}\sim \sqrt {m/{\sigma }}\sim {{O}(9\times 10^{-3}\ {\rm s})}$. However, a complete air crater retraction time scale ($T$) is a monotonic function of impact Weber number ($T\sim We^{1/2}$). We further demonstrate that the immiscible impact problem could be mapped to a constrained Rayleigh droplet model by incorporating a dynamic contact line relative to the impacting droplet. The crater retraction time scale increases with an increase in impact Weber number. When the air crater dynamics subsides, the droplet penetrates slowly under the influence of weight and buoyancy and attains an almost steady-state submerged configuration.

Figure 2.

Regime diagram showing the effect of impact Weber number on the evolution of droplet and air crater dynamics. The timestamps mentioned are in milliseconds and the scale bar denotes 1 mm.

3.2. Droplet penetration dynamics (scales of various forces acting on the droplet)

Figure 3(a) depicts a schematic representation of an impacting droplet on a silicone-oil pool. The oil pool is shown in light yellow. The schematic is drawn at a time instant $t$ ($t>t_r$) greater than the air rupture time scale, and the droplet is at a partially submerged state. Here, $CM$ represents the centre of mass of the impinging droplet. The contact line of three different fluids (air, water and oil) is marked as a white dotted line. The submerged portion/volume of the droplet immersed below the oil surface is quantitatively characterized and measured by the penetration depth $h(t)$ and penetration width $w(t)$ (figure 3a). We assume that the droplet deforms into an ellipsoid during the penetration process. The ellipsoid is characterized by the semi-major axis $a(t)$ and the semi-minor axis $b(t)$, respectively (figure 3a). The forces acting on the droplet are the weight $F_{mg}$ downwards, buoyancy force $F_b$ upwards and the viscous drag force $F_v$ upwards. The coordinate axis $z$ is chosen to be positive downwards. From mass conservation of the impacting droplet we have

(3.1)\begin{equation} {\rho}_w\tfrac{4}{3}{\rm \pi} R_0^3={\rho}_w\tfrac{4}{3}{\rm \pi}a^2(t)b(t)=m,\end{equation}

where ${\rho }_w$ is the density of the impacting droplet (water here), $R_0$ is the initial radius of the droplet and $m$ is the mass of the impacting droplet. For droplet radius of $R_0=1.1\ {\rm mm}$, the mass of the impinging water droplet is $m\sim {O}({5.58\times 10^{-6}}\ {\rm kg})$. Simplifying (3.1) we have

(3.2)\begin{equation} a^2(t)b(t)=R_0^3.\end{equation}

The dynamical equation (Newton's second law of motion) for the impinging liquid droplet is given by

(3.3)\begin{equation} F_{mg} - F_b(t) - F_v(t) = m\frac{{\rm d}V_{CM}(t)}{{\rm d}t} = ma_{CM},\end{equation}

with the initial condition for $V_{CM}$ as

(3.4)\begin{equation} V_{CM}(t=0)=V_0, \end{equation}

where $V_{CM}$ is the centre of mass velocity, $V_0$ is the impact velocity. The parameter $t=0$ defines the time instant just before the penetration of the droplet, $F_{mg}=mg$ is the weight of the impacting droplet where $g=9.81\ {\rm m}\ {\rm s}^{-2}$ is the acceleration due to gravity, $F_b$ is the buoyant force on the droplet and $F_v$ is the viscous drag force on the droplet due to the interaction with the silicone-oil pool. For $R_0=1.1\ {\rm mm}$, the weight of the droplet is of the order of $F_{mg}\sim {O}(5.47\times 10^{-5}\ {\rm N})$. The submerged droplet volume in the silicone-oil pool can be approximated by an ellipsoidal cap approximation (refer to figure 3a)

(3.5)\begin{equation} v(t) = \frac{{\rm \pi} a^2(t)h^2(t)}{3b^2(t)}(3b(t)-h(t)), \end{equation}

where $h(t)$ is the depth of penetration of the droplet. The buoyancy force is the weight of the displaced silicone oil, and therefore $F_b(t)$ becomes

(3.6)\begin{equation} F_b(t)=v(t){\rho}_sg=\frac{{\rm \pi} a^2h^2{\rho}_sg}{3b^2}(3b-h). \end{equation}

Using (3.2) in (3.6) the buoyancy force becomes

(3.7)\begin{equation} F_b(t) = \frac{{\rho}_sg{\rm \pi} R_0^3h^2}{3b^3}(3b-h). \end{equation}

At $h=b$ (i.e. when the droplet (ellipsoid) is half-submerged) the buoyancy force becomes

(3.8)\begin{equation} F_b = \tfrac{2}{3}{\rho}_s{\rm \pi}R_0^3g\sim {O}(2.65\times 10^{{-}5}\ {\rm N}). \end{equation}

Notice that the order of magnitude of the buoyancy force is the same as the weight of the droplet, and the buoyancy force becomes important at long time scales to determine the maximum steady-state submerged depth of the droplet. During the early stage of impact ($t\sim 2\ {\rm ms}$) the semi-major axis is approximately equal to the semi-minor axis ($a\simeq b$) and hence the orders of magnitudes of $a$ and $b$ are ${O}(R_0)$. Equation (3.6) can therefore be approximated for computing the buoyancy force at early impact time as

(3.9)\begin{equation} F_b(t)=v(t){\rho}_sg\sim\frac{{\rm \pi} h^2{\rho}_sg}{3}(3R_0-h). \end{equation}

The viscous force on the droplet $F_v$ can be approximated by the Stokes drag law (Moghisi & Squire Reference Moghisi and Squire1981) as given by

(3.10)\begin{equation} F_v(t)\sim 3{\rm \pi} \mu_sw(t)\frac{{\rm d}h}{{\rm d}t},\end{equation}

where ${\mu }_s$ is the viscosity of the silicone-oil pool. Figure 3(b) represents the scales of viscous force, buoyancy forces compared with the weight of the droplet plotted in a logarithmic scale against logarithmic time normalized with respect to the capillary time scale $t_c$ for impact Weber numbers of $4$, $16$ and $145$. Stokes’ approximation used in (3.10) becomes valid since the Reynolds number inside the silicone-oil pool is very small ($Re=V_0R_0/{\nu }_s\sim 1$). The maximum viscous force for $We\sim 16$ scales as $F_v\sim {O}(1.5\times 10^{-3}\ {\rm N})$ in the initial phase of penetration. Note that the viscous drag initially is a hundred times larger than the buoyant force and weight of the droplet. Therefore, the initial dynamics of droplet penetration and droplet response is governed by the Stokes viscous drag.

Figure 3.

(a) Schematic representation depicting various geometrical scales and forces acting on the impacting droplet while it penetrates the silicone-oil pool; (b) log–log plot depicting the strength of various forces in comparison with weight of the droplet as a function of time normalized with respect to capillary time scale $t_c=\sqrt {m/{\sigma }_{aw}}$, where $m$ is the mass of the droplet and ${\sigma }_{aw}$ is the air–water surface tension with impact Weber number as a parameter; (c,d) $f(t)/R_0^2$ as a function of time in milliseconds. The function $f(t)$ is a measure of a spherical drop shape.

3.3. Characterizing droplet deformation from spherical shape

In general, as the droplet penetrates the liquid pool, the droplet undergoes deformation and evolves as an ellipsoid (figure 3a). We have measured and characterized the amount of deviation from a spherical shape by calculating the function $f(t)$ defined as

(3.11)\begin{equation} f(t) = h^2(t) + \frac{w^2(t)}{4} - 2R_0h(t).\end{equation}

Here, $f(t)=0$ represents a spherical undeformed droplet while the penetration occurs. Larger non-zero values of $f(t)$ depict the deviation from spherical symmetry. Also, $f(t)$ has dimensions of length squared and therefore can be represented in a dimensionless way by normalizing with respect to $R_0^2$. Equation (3.11) therefore becomes

(3.12)\begin{equation} \frac{f(t)}{R_0^2}=\frac{h^2(t)}{R_0^2}+\frac{w^2(t)}{4R_0^2} - 2\frac{h(t)}{R_0}.\end{equation}

Figures 3(c) and 3(d) show the evolution of $f(t)/R_0^2$ plotted as a function of time using (3.12) for various impact Weber numbers ($We=4,16,145$). It can be inferred from the variation of $f(t)$ that, for low Weber numbers ($We\sim 4$), the droplet shape during penetration is very close to spherical. However, on increasing the impact Weber number ($We\sim 16,145$) the droplet deformation is substantial as can be observed by the non-zero values of $f(t)/R_0^2>1$. Regardless of the impact Weber number, beyond the capillary time $t_c$, the droplet penetrates, maintaining a spherical shape. It could be inferred that the deviation from the spherical shape due to penetration is a short time scale phenomenon ($t< t_c$). The dominant force experienced by the droplet on penetration is the viscous force (refer to figure 3b). We observe the formation of significant air craters due to sudden viscous drag acting on the droplet beyond critical impact Weber number of $10$. Figure 4(a) depicts the geometrical characteristics of the deforming droplet/air craters formed during penetration for three different impact Weber numbers equal to $4$, $16$ and $145$. We observe that significant air craters do not form for low Weber number ($We\sim 4$ shown here) due to delayed air layer rupture ($t_r\sim {O}(2\times 10^{-1}\ {\rm s})$). The impulsive nature of the drag force does not manifest during the penetration process due to delayed air layer rupture, and the viscous force becomes gradual for $We< We_c$. Air craters of significant shape and size are observed for larger Weber numbers ($We\sim 16$ and $We\sim 145$ shown here).

Figure 4.

(a) Snapshots depicting the penetration depth ($h(t)$) and width ($w(t)$) for impact Weber numbers of $4$, $16$ and $145$; (b) ${\varGamma }(t)=\Delta p_s/\Delta p_c$ as a function of time ($t_r< t< t_c$) for impact Weber numbers of $4$, $16$ and $145$.

3.4. Air crater formation criterion

The air craters observed for $We>We_c$ (refer to figure 4a) occur due to impulsive viscous drag force acting on the droplet across the immersed surface/volume in the silicone oil bounded by the contact line (refer to the schematic shown in figure 3(a). We introduce a non-dimensional parameter ${\varGamma }(t)$ as the ratio of the stokes drag pressure to capillary pressure across the air–water interface in the air crater formation time scale ($t_r< t< t_c$). The parameter can be represented as

(3.13)\begin{equation} \varGamma(t)=\frac{\Delta p_s}{\Delta p_c} \sim \frac{{\mu}_sw(t)V_0R_0}{{\sigma}_{aw}w^2(t)}, \end{equation}

where $\Delta p_s\sim \mu _sw(t)V_0/w^2(t)$ is the pressure force due to the stokes drag and $\Delta p_c\sim \sigma _{aw}/R_0$ is the capillary pressure across the air–water interface. Here, ${\varGamma }(t)\sim 1$ signifies that the impulsive drag is almost balanced by the capillary force acting across the air–water interface of the droplet.

Detectable air craters during the penetration process can only be observed if the viscous force becomes comparable to and overcomes the capillary force across the droplet interface, which corresponds to ${\varGamma }(t)>1$. We also show that, for $\varGamma (t)\leqslant 1$ for $t< t_c$, no significant air craters are observed, as is the case for low impact Weber number. Increasing Weber number increases ${\varGamma }(t)$ and hence the ratio of the Stokes viscous force to capillary force across the air–water interface of the droplet. Equation (3.13) can be rewritten in terms of the capillary number $Ca={\mu }_sV_0/{\sigma }_{aw}$ as

(3.14)\begin{equation} {\varGamma}(t)=\frac{\Delta p_s}{\Delta p_c} \sim Ca\frac{R_0}{w(t)}. \end{equation}

Figure 4(b) represents ${\varGamma }(t)=\Delta p_s/\Delta p_c$ as a function of time for $t_r< t< t_c$ for Weber numbers $4$, $16$ and $145$. It can be inferred from figure 4(b) that significant air craters are formed for ${\varGamma }>1$ where the impulsive viscous forces become larger than the capillary force across the air–water interface. For small Weber number, $We< We_c=10$, $\varGamma \lesssim 1$ signifying that the viscous forces are smaller than the capillary force across the air–water interface due to surface tension.

3.5. Evaluating the centre of mass velocity ($V_{CM}$)

We observe that the order of magnitude of the air crater jet velocity is comparable to the centre of mass velocity of the droplet ($V_{CM}$). Using (3.9) and (3.10) for the buoyancy and viscous terms, respectively, in (3.3), the centre of mass acceleration of the droplet can be expressed as

(3.15)\begin{equation} \frac{{\rm d}V_{CM}(t)}{{\rm d}t}=\frac{1}{m}\left(mg - \frac{{\rm \pi}h^2{\rho}_sg}{3}(3R_0-h) - 3{\rm \pi}{\mu}_sw(t)\frac{{\rm d}h}{{\rm d}t}\right). \end{equation}

Integrating (3.15) we have

(3.16)\begin{equation} \int_{V_0}^{V_{CM}(t)}{\frac{{\rm d}V_{CM}(t)}{{\rm d}t}}{\rm d}t= \int_{0}^{t}\frac{1}{m}\left(mg - \frac{{\rm \pi}h^2{\rho}_sg}{3}(3R_0-h) - 3{\rm \pi}{\mu}_sw(t)\frac{{\rm d}h}{{\rm d}t}\right){\rm d}t. \end{equation}

Simplifying (3.16) we have

(3.17)\begin{equation} V_{CM}(t) = V_0 + gt - \frac{{\rm \pi}{\rho}_sg}{3}\int_0^th^2(3R_0-h)\,{\rm d}t - 3{\rm \pi}{\mu}_s\int_0^tw(t)\frac{{\rm d}h}{{\rm d}t}{\rm d}t.\end{equation}

Notice that the acceleration due to gravity term $gt$ and the buoyancy term (third term on the right-hand side) are negligible compared with the viscous term (last term on the right-hand side) for $t_r< t< t_c$. The last term representing the viscous forces is the most dominant for $t_r< t< t_c$.

Figures 5, 6 and 7 show various kinematic and dynamic quantities for impact Weber numbers of $4$, $16$ and $145$, respectively. Figure 5(a) depicts the high-speed snapshots during the penetration process for an impact Weber number of $4$. The penetration dynamics for low Weber numbers is drastically different as the droplet and air layer interaction is significant. The air layer ruptures and the formation of a subsequent bubble can be observed in figure 5(a). Further, as the droplet penetrates through the silicone-oil pool, the droplet shape is close to spherical (refer to figure 3c,d). Figure 5(b) shows the temporal evolution of submerged droplet normalized width ($w(t)/R_0$) and normalized depth ($h(t)/2R_0$). Using the experimental values of $w(t)$ and $h(t)$, we have calculated the resultant force acting on the droplet. All the individual forces ($F_{mg}$, $F_b$ and $F_v$) along with the resultant force (in red $F_{mg}+F_b+F_v$) is plotted in figure 5(c). Notice that, during the early time of penetration for $t< t_c$, the dominant force is only the viscous drag. In general, we observe that the resultant force is almost entirely due to viscous drag with buoyancy and weight almost balancing each other as time increases.

Figure 5.

(a) High-speed snapshots depicting the evolution of air layer rupture forming a bubble during penetration through silicone oil for impact Weber number $We=4$. The timestamps are in milliseconds and the scale bar represents $500\ \mathrm {\mu }{\rm m}$. (b) Temporal evolution of kinematic parameters like penetration width and depth. (c) Temporal evolution of the various forces acting on the impacting droplet. The resultant of the various forces is also shown in red. Refer to supplementary movie 1 for $We=4$ available at https://doi.org/10.1017/jfm.2022.1071.

Figure 6.

(a) High-speed snapshots depicting the evolution of air layer rupture forming a bubble during penetration through silicone oil for impact Weber number $We=16$. The timestamps are in milliseconds and the scale bar represents $500\ \mathrm {\mu }{\rm m}$. (b) Temporal evolution of kinematic parameters like penetration width and depth. (c) Temporal evolution of the various forces acting on the impacting droplet. The resultant of the various forces is also shown in red. Refer to supplementary movie 2 for $We=16$ here.

Figure 7.

(a) High-speed snapshots depicting the evolution of air layer rupture forming a bubble during penetration through silicone oil for impact Weber number $We=145$. The timestamps are in milliseconds and the scale bar represents $500\ {\mathrm {\mu }}{\rm m}$. (b) Temporal evolution of kinematic parameters like penetration width and depth. (c) Temporal evolution of the various forces acting on the impacting droplet. The resultant of the various forces is also shown in red. Refer to supplementary movie 3 for $We=145$ here.

Using the resultant external forces acting on the droplet (figure 5c), the centre of mass velocity ($V_{CM}$) can be evaluated using (3.17). We evaluate centre of mass velocity ($V_{CM}$) for early penetration time ($t< t_c$) as the approximation under which (3.17) is derived becomes increasingly accurate for small time duration. The effect of the impact Weber number can be understood by analysing the scale of $F_v$ as a function of the Weber number. It can be inferred from figure 3(b) that, on increasing the impact Weber number, the viscous force $F_v$ changes significantly. The resultant force on the droplet is approximately equal to that of the viscous force. Therefore, on increasing the impact Weber number, the ratio of the Stokes viscous pressure to capillary pressure of the droplet ${\varGamma }$ increases, resulting in the formation of air crater/reversed jets of various shapes and sizes (refer to figure 4a,b).

Figure 6(a) shows a sequence of high-speed snapshots showing the evolution of the air crater for $We=16$. Figure 6(b) plots the evolution of normalized penetration depth and width of the impacting droplet for an impact Weber number of $16$. The individual and resultant forces acting on the droplet are shown in figure 6(c).

Figure 7(a) shows the air crater formation and evolution for impact Weber number of $145$. The crater formed is significantly different compared with the low Weber number impact. For a Weber number of $16$, a single air jet is formed due to the capillary waves interfering. Whereas for a Weber number of $145$, a cascade of capillary waves is formed, forming a bigger and more intricate air crater. The air crater normalized depth ($h(t)/2R_0$) and normalized width ($w(t)/R_0$) are plotted as a function of time in figure 7(b). Using the experimental air crater geometrical parameters, the individual and the resultant force acting on the impacting droplet were calculated and are shown in figure 7(c). The resultant force ${\varSigma }F(t)$ acting on the impinging droplet during the penetration process is plotted as a function of time for various impact Weber numbers in figure 8(a) for comparison purposes ($We=416\,145$). The centre of mass velocity ($V_{CM}$), the solution for (3.17), is plotted in figure 8(b) for various impact Weber numbers.

Figure 8.

(a) Comparison of the resultant force ${\varSigma }F(t)$ in mN acting on the droplet during the penetration process as a function of time $t$ in ms for impact Weber numbers ($We$) of $4$, $16$ and $145$ evaluated using the right-hand side of (3.15). (b) Comparison of the centre of mass velocity of the droplet ($V_{CM}$) in ${\rm m}\ {\rm s}^{-1}$ plotted as a function of time $t$ in ms for impact Weber numbers ($We$) of $4$, $16$ and $145$. (c) A high-speed snapshot depicting the crater depth $h_c$. (d) Crater depth ($h_c$) evolution as a function of time in ms.

The centre of mass velocity decreases as the droplet experiences retarding forces as it penetrates the liquid silicone-oil pool. We can further observe that the centre of mass velocity shows a damped oscillatory response for $t\sim t_c$. Figure 8(c) shows a typical air crater geometry and its characteristics scale (air crater depth $h_c$) formed on the surface of an impacting droplet. The evolution of air crater depth $h_c$ as a function of time $t$ is plotted in figure 8(d) for impact Weber numbers of $16$, $65$ and $145$. Note that the initial slope of the $h_c$ vs $t$ plot is ordered according to the centre of mass velocity for various Weber numbers; $We=145$ has the highest initial slope followed by $We=65$ and $We=16$.

3.6. Transient droplet response (local analysis)

The air crater/jets are formed due to the transient response of the droplet to impulsive forces acting on the droplet. The sudden deceleration of the droplet induces droplet oscillation that causes air craters of various shapes and sizes to form. Further, relative to the droplet the external force acting could also be formulated based on the relative effect of the surrounding air field. The transient response of the droplet surface based on sudden retardation can be characterized by surface displacement ${\eta }=r({\mu },t)/R_0$; where $r({\mu },t)$ is the radial location of the droplet interface at any angular coordinate ${\theta }$ at time $t$ (refer to figure 3a) and $R_0$ is the initial droplet radius. We closely follow this part of the analysis based on the work of Harper et al. (Reference Harper, Grube and Chang1972) and Simpkins & Bales (Reference Simpkins and Bales1972).

The perturbation result for the surface displacement based on a linear theory proposed by Harper et al. (Reference Harper, Grube and Chang1972) is given as

(3.18)\begin{equation} {\eta} = 1 + \epsilon\sum_{n=0}^{\infty}\frac{n(2n+1)}{4{\beta}_n^2}C_nP_n({\mu})[\cos({\beta}_nt)-1], \end{equation}

where $\epsilon =\rho _a/\rho _w$ is the density ratio of air to water characterizing the amplitude of the perturbation, $n$ is a positive integer and $C_n$ are the weighting coefficients given by

(3.19)\begin{equation} C_n = \int{P}_e^{(0)}P_n({\chi})\,{\rm d}\chi, \end{equation}

where $\chi =\cos \psi$. Here, $P_n(\mu )$ represents the $n$th-order Legendre polynomials where $\mu =cos\theta$ and

(3.20)\begin{equation} \beta_n={\left(\frac{(n-1)n(n+2)}{We_{R_0}}\right)^{1/2}}, \end{equation}

where ${\beta }_n$ denotes the frequency of $n$th Legendre mode and $We_{R_0}$ denotes the Weber number defined with respect to the initial radius $R_0$ ($We = 2 We_{R_0}$). Here, $P_e^{(0)}$ represents the leading term of the external pressure $P_e$ in the perturbation expansion with ${\epsilon }$ as the perturbation parameter acting on the droplet during the penetration process

(3.21)\begin{equation} \beta\sim We^{{-}1/2}. \end{equation}

Refer to figure 3(a) for the definition of various angular coordinates like ${\theta }$ and ${\psi }$. Note that, near the droplet top surface, that corresponds to ${\psi }=0$ and ${\theta }={{\rm \pi} }$, the droplet shape approximates that of a Legendre polynomial locally. Figures 9(a) and 9(b) show the comparison of the air jet profile at the centre with the fourth-order Legendre polynomial for $We=16$. Similar comparison is shown in figures 9(c) and 9(d) for $We=145$. In principle, the combined air crater deformation can be thought of as the superposition of various Legendre modes. The period of half-oscillation of the air jets scales as

(3.22)\begin{equation} T\sim \beta^{{-}1}\sim We^{1/2}.\end{equation}

Figure 9.

(a) Image depicting the profile of the reversed jet formed during air crater formation/collapse at impact Weber number of 16. (b) Comparison of the experimental profile of the central jet with fourth degree Legendre polynomial. (c) Image depicting the profile of the reversed jet formed during air crater formation/collapse at impact Weber number of 145. (d) Comparison of the experimental profile of the central jet with fourth degree Legendre polynomial. (e) Experimental and theoretical comparison of air crater time scale ($T$) as a function of impact Weber number ($We$) based on (3.22).

Therefore, from (3.22), we observe that the period of the air crater dynamics increases monotonically with the impact Weber number. As inferred from figure 9(e), we observe the corresponding monotonic behaviour experimentally.

3.7. Transient vibrational response of the droplet based on constrained Rayleigh drop model (global analysis)

In this section, we generalize the local analysis performed in § 3.6 with a global response model. When the entrapped air beneath the impacting droplet ruptures and forms a central air bubble, a three-fluid dynamic contact line is formed (refer to figures 1c, 3a, 10a and 10b). The contact line is a one-dimensional boundary separating the three immiscible fluids (air, water and silicone oil) and is one of the important aspects controlling the response of the droplet. The dynamic contact line separates the droplet interface into two disjoint curved surface areas. The dynamic contact line is characterized using an inclination angle ${\alpha }$ measured from the $-y$ axis. Figures 10(a) and 10(b) graphically represent the geometry of a penetrating droplet from the constrained Rayleigh drop perspective. Figure 10(a) represents the droplet penetration at any given instant for an acute value of ${\alpha }={{\rm \pi} /4}$. For ${\alpha }={{\rm \pi} /2}$, the schematic is shown in figure 10(b). The saffron coloured region represents the silicone-oil pool and $R=R_0$ represents the radius of the undeformed equilibrium shape of the droplet. The dynamic contact line is represented by the horizontal dotted blue line. The curved dotted blue line represents the bottom submerged portion of the droplet. The red dotted curved interface represents the equilibrium position of the air–water interface. The blue and green solid lines depict the oscillating air–water interface of the droplet that leads to the formation of air craters and $z({\theta })$ represents the amplitude of the air–water interface from the equilibrium position of the interface at a given polar angle ${\theta }$ and time $t$. In general, the dynamic contact line imposes the constraint of zero amplitude along its perimeter. However, owing to the large viscosity of the silicone pool, the entire curved surface of the droplet submerged in the silicone pool (water–oil interface) could be considered as an surface area of negligible (zero) amplitude. The amplitude response of the droplet from the mean equilibrium shape in general depends on the impact Weber number and increases monotonically with the same (refer figure 8d). For small impact Weber number the droplet interface amplitudes and the deviation from the average spherical equilibrium shape are relatively small (refer figure 3d). We therefore use the constrained Rayleigh drop model for such low energy impacts. We follow the analysis of Strani, Bostwick–Steen closely in this section (Strani & Sabetta Reference Strani and Sabetta1984; Bostwick & Steen Reference Bostwick and Steen2009). In a recent work, Kern et al. (Reference Kern, Bostwick and Steen2021) showed that, in drop impact systems on solids, contact angle hysteresis filters impact energy into modal interfacial vibrations.

Figure 10.

Graphical representation depicting the geometrical nomenclature and the air–water interfacial waves for the global deformation analysis. For (a) ${\alpha }={{\rm \pi} }/4$, (b) ${\alpha }={{\rm \pi} }/2$. (c) An experimental snapshot depicting the observed crater characteristics for impact Weber number $We=16$.

The continuity equation for incompressible flow inside the impacting droplet is given as

(3.23)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{V}=0. \end{equation}

The velocity field $\boldsymbol {V}$ inside the droplet could be extracted from a velocity potential ${\varPhi }$ provided that the flow inside the impacting droplet is inviscid. This could be done due to the fact that the impacting droplet is approximately $350$ times less viscous than the silicone-oil pool. The experimental response time scale of the droplet is within the capillary time scale ($t_c\sim {O}(9\times 10^{-3}\ {\rm s})$). The viscous effects inside the droplet are important at the scale of the viscous diffusion time scale (${\sim }{O}(4R_0^2/{\nu }_w)=5.43 {\rm {s}}$). We should note that the droplet response occurs at a scale of milliseconds, which is three orders of magnitude smaller than the viscous diffusion time scale inside the droplet. We therefore use a inviscid droplet model for computing the droplet response during the capillary time scale

(3.24)\begin{equation} \boldsymbol{V}=\boldsymbol{\nabla}{\varPhi}.\end{equation}

Substituting (3.24) in (3.23) we have

(3.25)\begin{equation} {\nabla}^2{\varPhi}=0, \end{equation}

where ${\nabla }^2$ is the Laplacian operator. We use spherical polar coordinates $(r,{\theta },{\varPsi })$, where $r$ represents the radial coordinate from the centre of the droplet, ${\theta }$ is the zenith/polar angle and ${\varPsi }$ is the azimuthal angle. We assume an axisymmetric geometry due to a cylindrical symmetry about the $y$-axis (figure 10a) in our experimental configuration. Note that the local analysis carried out in the previous section uses a different nomenclature for various angular coordinates ${\theta }$ and ${\psi }$ (refer to figure 3c). Equation (3.25) in spherical axisymmetric polar coordinates therefore becomes

(3.26)\begin{equation} \nabla^2\varPhi=\sin\theta(r^2\phi_r)_r+(\sin\theta\varPhi_{\theta})_{\theta}=0, \end{equation}

on

(3.27)$$\begin{gather} \mathscr{D}^i\equiv\{(r,{\theta})\,|\,0\leqslant r\leqslant R,0\leqslant \theta\leqslant {\rm \pi}\}, \end{gather}$$

(3.28)$$\begin{gather}\mathscr{D}^e\equiv\{(r,{\theta})\,|\,R\leqslant r<\infty,0\leqslant\theta\leqslant{\rm \pi}\}. \end{gather}$$

We use subscripts $r$, ${\theta }$ in (3.26) to represent partial derivatives with respect to radial coordinate ($r$) and the zenith angle (${\theta }$), respectively. This convention is used throughout the analysis in this section. $\mathscr {D}$ represents the two-dimensional axisymmetric equilibrium domain of the problem. The superscripts $i$, $e$ represents regions pertaining to the interior and exterior of the droplet equilibrium shape, respectively, and $\mathscr {D}^i$ in (3.27) represents the two-dimensional axisymmetric interior region of the equilibrium shape of the impacting droplet. Similarly, the exterior of the droplet in equilibrium conditions is given by $\mathscr {D}^e$ represented in (3.28). The pressure inside and outside of the droplet could be computed from the unsteady Bernoulli equation Kundu, Cohen & Dowling (Reference Kundu, Cohen and Dowling2015)

(3.29)\begin{equation} p(r,{\theta},t)=p_0+\rho\varPhi_t,\end{equation}

where $p_0$ is a reference pressure (the atmospheric pressure), and ${\rho }$ is the density of the fluid. The subscript $t$ refers to the partial derivative of ${\varPhi }$ with respect to time. The droplet interface is represented by $\partial \mathscr {D}$. The superscript $s$ will represent the submerged interface (water–silicone-oil interface) and the superscript $f$ represents the free surface of the droplet (air–water interface). On the submerged boundary portion of the droplet represented by

(3.30)\begin{equation} \partial\mathscr{D}^s\equiv\{(r,\theta)\,|\,r=R,{\rm \pi}-\alpha\leqslant\theta\leqslant{\rm \pi}\}, \end{equation}

the radial velocity ${\varPhi }_r$ and the amplitude $Z$ is zero. This approximation could be made due to the fact that the liquid pool is highly viscous and dissipates any interface oscillatory motion. The time scale of the viscous dissipation in the liquid pool is of the order of ${\sim }{O}(4R_0^2/{\nu }_s)\sim 10\times 10^{-3}s\sim t_c$

(3.31a,b)\begin{equation} \varPhi_r = 0;\quad Z=0. \end{equation}

On the free boundary portion of the droplet

(3.32)\begin{equation} \partial\mathscr{D}^f\equiv\{(r,{\theta})\,|\,r=R,0\leqslant \theta\leqslant {\rm \pi}-\alpha\},\end{equation}

the boundary condition at the interface preserving the continuity of normal velocity is given by

(3.33)\begin{equation} {\varPhi}_r ={-}Z_t. \end{equation}

The dynamic balance of momentum along the droplet radius gives

(3.34)\begin{equation} p^i - p^e = \sigma\left[\frac{2}{R}-\frac{1}{R^2}\left(\frac{1}{\sin\theta}\left(\sin\theta Z_{\theta}\right)_{\theta}+2Z \right)\right].\end{equation}

Recognize that the unknowns are the droplet deformation from equilibrium position given by $Z({\theta },t)$ and the velocity potential given by ${\varPhi }(r,{\theta },t)$

(3.35a,b)\begin{equation} Z = Z({\theta},t)=\,?;\quad \varPhi(r,{\theta},t)=\,?. \end{equation}

Conservation of mass for the impacting drop gives us

(3.36)\begin{equation} \int_0^{\rm \pi}Z\sin \theta \,{\rm d}\theta=0. \end{equation}

This is true if the droplet does not fragment into daughter droplets, which is true at low impact energies, as in the current experimental configuration. By considering periodic motion inside and on the free air–water interface, we can set

(3.37)$$\begin{gather} {\varPhi}(r,{\theta},t)={\phi}(r,{\theta})\exp({{\rm i}{\omega}t}), \end{gather}$$

(3.38)$$\begin{gather}Z({\theta},t) = z(\theta)\exp({{\rm i}({\omega}t+{\rm \pi}/2)}). \end{gather}$$

The phase difference of ${{\rm \pi} }/2$ is consistent with the boundary condition given by (3.33). Using the unsteady Bernoulli equation (3.29), the velocity potential given by (3.37) and the surface perturbation given by (3.38) in (3.34), the pressure difference term $(p^i-p^e)$ on the left-hand side of (3.34) is evaluated. The pressure difference calculation across the droplet interface and the simplification of (3.34) using the unsteady Bernoulli and the velocity potential is shown in Appendix A (refer to (A1)–(A15)). Further, on simplifying, (3.34) represents two equations from the real part and complex part of equation. The real part of (3.34) is given by (A14) and the complex part is given by (A15). On the two-dimensional axisymmetric domain $\mathscr {D}$ we have from (3.26)

(3.39)\begin{equation} {\nabla}^2{\phi}=0\end{equation}

and from (A3) we have

(3.40)\begin{equation} {\phi}_r=\omega z({\theta}). \end{equation}

On the free air–water interface of the droplet $\partial \mathscr {D}^f$ we have from (A15)

(3.41)\begin{equation} {\omega}[{\rho}^i{\phi}^i-{\rho}^e{\phi}^e]={-}\frac{\sigma}{R^2}\left[\frac{1}{\sin\theta} \left(\sin\theta Z_{\theta}\right)_{\theta}+2Z\right]. \end{equation}

On the submerged water–oil interface of the droplet $\partial \mathscr {D}^s$ we have from (3.31a,b), (3.33), (3.37) and (3.38)

(3.42a,b)\begin{equation} {\phi}_r=0,\quad z=0. \end{equation}

Using (3.38) in (3.36) we have

(3.43)\begin{equation} \int_0^{\rm \pi}z\sin\theta \,{\rm d}\theta=0. \end{equation}

Equations (3.39)–(3.43) define an eigenvalue problem in a continuous functional space. Further, (3.39) on $\mathscr {D}$ and (3.40) on $\partial \mathscr {D}$ define a Neumann problem. From the stagnant nature of the ambient at far field $r\rightarrow \infty$ we have

(3.44)\begin{equation} {\phi}^e{\rightarrow}0\quad {\rm as}\ r\rightarrow\infty. \end{equation}

Introducing the change of variable ${\mu }=\cos \theta$, the velocity potential solution for the interior and exterior of the droplet could be written as infinite series using the Legendre polynomials as the basis functions. The velocity potential therefore could be written as

(3.45)$$\begin{gather} {\phi}^i(r,{\mu})=\omega R\left({\phi}_0+\sum_{k=1}^{\infty}\frac{{\phi}_k}{k}\frac{r^k}{R^k}P_k({\mu})\right) \end{gather}$$

(3.46)$$\begin{gather}{\phi}^e(r,{\mu})={-}\omega R\sum_{k=1}^{\infty}\frac{{\phi}_k}{k+1}\frac{R^{k+1}}{r^{k+1}}P_k({\mu}), \end{gather}$$

where

(3.47)$$\begin{gather} {\phi}_k=\frac{\langle z,P_k\rangle }{\langle P_k,P_k\rangle }, \end{gather}$$

(3.48)$$\begin{gather}\langle\, f,g\rangle =\int_{{-}1}^{{+}1}fg\,{\rm d}{\mu}, \end{gather}$$

(3.49)$$\begin{gather}\langle P_k,P_k\rangle=\frac{2}{2k+1}. \end{gather}$$

Here, $P_k({\mu })$ denotes the $k$th Legendre polynomial and $\langle\, f,g\rangle$ represents the inner product of $f$ and $g$ and is square integrable in the domain ($L^2(-1,1)$). Equation (3.41) using the definition of ${\mu }$ could be written as

(3.50)\begin{equation} [(1-{\mu}^2)z{\mu}]_{\mu}+2z={-}\frac{{\rho}^i\omega R^2}{\sigma}\left[{\phi}^i(R,{\mu}) - \frac{{\rho}^e}{{\rho}^i}{\phi}^e(R,{\mu})\right]=f({\mu}), \end{equation}

for $a\leqslant \mu \leqslant 1$. Here, $a=-\cos ({\alpha })$ (refer to figure 10a). The solution to (3.50) is provided in Appendix B (refer to (B1a,b)–(B14)) in terms of a functional eigenvalue problem.

The continuous functional eigenvalue problem represented by (B13) could be reduced to an infinite set of algebraic equations by expanding $z({\mu })$ and $K({\mu },{\sigma })$ (refer to (C13) and (C14)) in the Legendre polynomial basis (refer to Appendix C). Equation (C8) describes an infinite-dimensional matrix-vector eigenvalue problem where $A_{lh}$ represents the matrix, $x_l$ or $x_h$ represents the eigenvector and ${\lambda }$ represents the eigenvalue corresponding to the eigenvector $x_h$. The elements of matrix $A$ are described in Appendix D. Equations (D21), (D22) and (D23) depict that $A$ is a symmetric real diagonal matrix with eigenvalues as the elements of the diagonal. The eigenvalues corresponding to matrix $A$ in limiting conditions are given as

(3.51a–c)\begin{gather} {\lambda}^1\rightarrow+\infty,\quad x_k^1\rightarrow \delta_{1k},\quad z_k^1\rightarrow \delta_{1k}, \end{gather}

(3.52a–c)\begin{gather} {\lambda}^n\rightarrow\frac{(n+1)+({\rho}^e/{\rho}^i)n}{n(n-1)(n+1)(n+2)};\quad x_k^n\rightarrow \delta_{nk},\quad z_k^n\rightarrow \delta_{nk}. \end{gather}

From (3.51a–c) and (3.52a–c) we observe that the first eigenmode corresponds to zero frequency which corresponds to a rigid body motion, whereas the various other eigenmodes tend towards pure Rayleigh modes for the case of negligible pinning at and below the dynamic contact line. The infinite-dimensional eigenvalue problem represented by (C8) is solved numerically using in-house program written in python using a ten-dimensional Legendre function space approximation. The elements of the matrix $A$ were represented using (D1). Note that the matrix $A$ is real and symmetric. Figure 11(a) depicts a visual representation of the matrix $A$. Note that the matrix $A$ is diagonal. Further, $A$ being real and symmetric has real eigenvalues and the eigenvalues are the ordered diagonal elements. Matrix $A$ represent a self-adjoint operator in an orthonormal basis of Legendre polynomials forming a real inner product space.

Figure 11.

(a) Visual representation of matrix $A$ depicting its elements. (b) Visual representation of the eigenvector matrix of $A$. (c) Bar plot depicting the strength (eigenvalues) of various eigenmodes of matrix $A$. (d) Experimental and theoretical comparison of the air–water interface crater profile for impact Weber number $We=16$.

Figure 10(b) depicts the eigenvector matrix for $A$. The columns of the matrix represent each eigenvector while the eigenvector matrix being diagonal represents the orthogonal nature of the eigenvectors in the function space. Figure 11(c) represents the eigenvalues of the matrix. The eigenvalues correspond to the various modes of the Legendre polynomials used in the approximate function space for calculating the amplitude $z({\theta })$, ${\omega }$ and ${\phi }(r,{\theta })$. Figures 11(a) and 11(c) depict that modes $3$, $4$ are the most dominant ones and can approximate the amplitude, frequency and velocity potential functions. Figure 11(d) shows the comparison of the theoretical curve obtained from the global analysis with the experimental data at the maximum penetration width for an impact Weber number of $16$. The experimental and theoretical curves conform within the experimental error bounds. Further, the global analysis agrees with the local analysis outlined in the previous § 3.6. For larger impact Weber number, extreme air–crater deformations occur (figure 7) and the above analysis needs to be generalized for high amplitudes ${z({\theta })}$. Further, in a general scenario the effects of viscosity have to be incorporated in the shape model and the pinned boundary conditions need to be applied on the dynamic one-dimensional contact line.