Insights into air cushion dynamics during drop impact on heated substrate at low impact energy

TL;DR

This study investigates the air cushion behavior beneath a droplet impacting a heated surface at low energy, revealing distinct flow regimes and asymmetric wetting phenomena through high-speed imaging and theoretical analysis.

Contribution

It introduces a combined continuum and non-continuum framework to understand the air layer dynamics during droplet impact on heated substrates, highlighting the role of gas kinetic effects.

Findings

Air film subdivided into central dimple and peripheral disc regions.

Dimple shape approximated by a Gaussian profile, weak dependence on temperature.

Air layer rupture time and radius increase with substrate temperature.

Abstract

We study the air layer dynamics beneath a drop impinging a heated surface at low impact energy using high-speed reflection interferometry imaging and theoretical analysis. The air film has been subdivided into two distinct disjoint regions, the central dimple and the peripheral disc. We decipher that a gaussian profile can approximate the dynamic shape evolution of the central air dimple. We further observe that the dimple geometry is a function of impact energy and its dependence on surface temperature is relatively weak. The air layer rupture time and rupture radius increases with increase in substrate temperature. We characterize the air layer profile as a 2D Knudsen field and show that a unified treatment, including continuum and non-continuum mechanics, is required to comprehend the air layer dynamics coherently. The airflow dynamics in the central dimple region falls within the purview of continuum stokes regime. In contrast, the peripheral air disc falls within the non-continuum (gas kinetic effects) slip flow and transition regime characterized by a high Knudsen number. However, the initial average air disc expansion dynamics could be understood in terms of stokes approximation. In non-continuum regimes of the peripheral air disc, we discover intriguing asymmetric interface perturbations. The asymmetric wetting of the substrate initiates at the edge of the peripheral disc region.These perturbative structures cause asymmetric wetting/contact between the droplet and the substrate. Due to the asymptotic effects of capillary and van der Waals interaction in the disc region, the sub-micron spatial structures can exist at short time scales.

Equations

$$I(x,y)=I_1(x,y)+I_2(x,y)+2\sqrt{I_1(x,y)I_2(x,y)}cos(\varphi (x,y))$$

$$I(x,y)=B(x,y)+A(x,y)cos(\varphi (x,y))$$

$$F(x,y)=\frac{I(x,y)-B(x,y)}{A(x,y)}=cos(\varphi (x,y))$$

$${\varphi }_0(x,y)=co{s}^{-1}(F(x,y))$$

$$C(x,y)=\sum _{({\delta }_1,{\delta }_2)\in N_{xy}}[F({\delta }_1,{\delta }_2)-cos(\varphi (x,y;{\delta }_1,{\delta }_2)){]}^2$$

$$\varphi (x,y;{\delta }_1,{\delta }_2)={\varphi }_0(x,y)+{\omega }_x(x,y)({\delta }_1-x)+{\omega }_y(x,y)({\delta }_2-y)$$

$$\nabla F(x,y)\cdot \nabla \varphi (x,y)=-sin(\varphi (x,y))|\nabla \varphi (x,y){|}^2$$

$$\nabla {\varphi }_0(x,y)=\omega (x,y)e^{i\theta (x,y)}$$

$$h(x,y)=\frac{{\lambda }_s}{2\pi n}\varphi (x,y)$$

$$sin\varphi =\frac{r}{R_0};cos\varphi =\frac{\sqrt{R_0^2-r^2}}{R_0}$$

$$z(r)=h_0+R_0\left(1-\sqrt{1-{\left(\frac{r}{R_0}\right)}^2}\right)$$

$$z(r)=R_0ϵ+\frac{1}{2}\frac{r^2}{R_0}$$

$$\bar{z}=\frac{z}{R_0ϵ}=1+\frac{r^2}{2R_0^2ϵ}$$

$$\bar{z}=h(\bar{r})=\frac{z}{R_0ϵ}=1+\frac{1}{2}{\bar{r}}^2$$

$$\frac{1}{r}\frac{\partial }{\partial r}\left(r{u}_r\right)+\frac{\partial u_z}{\partial z}=0$$

$$\frac{{ϵ}^{1/2}}{V_0}\frac{1}{\bar{r}}\frac{\partial }{\partial \bar{r}}\left(\bar{r}{u}_r\right)+\frac{\partial {\bar{u}}_z}{\partial \bar{z}}=0$$

$$\frac{1}{\bar{r}}\frac{\partial }{\partial \bar{r}}\left(\bar{r}{\bar{u}}_r\right)+\frac{\partial {\bar{u}}_z}{\partial \bar{z}}=0$$

$${\rho }_a\left(\frac{\partial u_r}{\partial t}+u_r\frac{\partial u_r}{\partial r}+u_z\frac{\partial u_r}{\partial z}\right)=-\frac{\partial p}{\partial r}+{\mu }_a\left(\frac{{\partial }^2{u}_r}{\partial z^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_r}{\partial r}\right)\right)$$

$$R{e}_{ϵ}\left(\frac{\partial {\bar{u}}_r}{\partial \bar{t}}+{\bar{u}}_r\frac{\partial {\bar{u}}_r}{\partial \bar{r}}+{\bar{u}}_z\frac{\partial {\bar{u}}_r}{\partial \bar{z}}\right)=-\frac{\partial \bar{p}}{\partial \bar{r}}+\left(\frac{{\partial }^2{\bar{u}}_r}{\partial {\bar{z}}^2}+\frac{ϵ}{\bar{r}}\frac{\partial }{\partial \bar{r}}\left(\bar{r}\frac{\partial {\bar{u}}_r}{\partial \bar{r}}\right)\right)$$

$$\frac{{\partial }^2{\bar{u}}_r}{\partial {\bar{z}}^2}=\frac{\partial \bar{p}}{\partial \bar{r}}$$

$${\rho }_a\left(\frac{\partial u_z}{\partial t}+u_r\frac{\partial u_z}{\partial r}+u_z\frac{\partial u_z}{\partial z}\right)=-\frac{\partial p}{\partial z}+{\mu }_a\left(\frac{{\partial }^2{u}_z}{\partial z^2}+\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_z}{\partial r}\right)\right)$$

$$R{e}_{ϵ}ϵ\left(\frac{\partial {\bar{u}}_z}{\partial \bar{t}}+{\bar{u}}_r\frac{\partial {\bar{u}}_z}{\partial \bar{r}}+{\bar{u}}_z\frac{\partial {\bar{u}}_z}{\partial \bar{z}}\right)=-\frac{\partial \bar{p}}{\partial \bar{z}}+ϵ\left(\frac{{\partial }^2{\bar{u}}_z}{\partial {\bar{z}}^2}+\frac{ϵ}{\bar{r}}\frac{\partial }{\partial \bar{r}}\left(\bar{r}\frac{\partial {\bar{u}}_z}{\partial \bar{r}}\right)\right)$$

$$\frac{\partial \bar{p}}{\partial \bar{z}}=0$$

$${\bar{u}}_r(\bar{r},\bar{z})=\frac{\partial \bar{p}}{\partial \bar{r}}{\bar{z}}^2+C_1(\bar{r})\bar{z}+C_2(\bar{r})$$

$${\bar{u}}_r={\bar{u}}_z=0;at\bar{z}=0$$

$${\bar{u}}_r=0;{\bar{u}}_z=-1;at\bar{z}=h(\bar{r})$$

$${\bar{u}}_r=\frac{\partial \bar{p}}{\partial \bar{r}}\left(\frac{{\bar{z}}^2}{2}-\frac{\bar{z}h(\bar{r})}{2}\right)$$

$${\int }_0^{h(\bar{r})}\frac{1}{\bar{r}}\frac{\partial \bar{r}{\bar{u}}_r}{\partial \bar{r}}d\bar{z}+{\int }_0^{h(\bar{r})}\frac{\partial {\bar{u}}_z}{\partial \bar{z}}d\bar{z}=0$$

$${\int }_0^{h(\bar{r})}\frac{1}{\bar{r}}\frac{\partial \bar{r}{\bar{u}}_r}{\partial \bar{r}}d\bar{z}+{\bar{u}}_z{|}_{\bar{z}=h(\bar{r})}-{\bar{u}}_z{|}_{\bar{z}=0}=0$$

$${\int }_0^{h(\bar{r})}\frac{1}{\bar{r}}\frac{\partial \bar{r}{\bar{u}}_r}{\partial \bar{r}}d\bar{z}+(-1-0)=0$$

Full text

1 Introduction

Drop impact phenomena is quite ubiquitous and spectacular. Appreciating the elegant structures that appear during drop impact dynamics on solids and liquids does not require scientific training. This phenomenon thus appears in a variety of advertisements/artistic contexts in the media industries. The physics underlying the aesthetic beauty of drop impact phenomena is equally beautiful and attracts the scientific community to study it in various contexts. The original study of drop impact physics began with the seminal work of Worthington (Worthington, 1877; Worthington et al., 1877; Worthington, 1883) in the late nineteenth century. Using ingenious mechanical arrangements to synchronize the fall of the drop with that of an electric discharge, Worthington was able to study drop impact phenomena on solids and liquids. Initial observations made by Worthington were manual and hand drawn. Later on, Worthington was able to record the impact phenomena on photographic plates (Worthington & Cole, 1897, 1900).

With the advent of high-speed imaging, experimental drop impact physics was revolutionized as scientists and engineers could probe smaller and finer spatio-temporal scales (Thoroddsen et al., 2008; Yarin et al., 2017). Experimental drop impact research was majorly driven by various technologies and industrial processes like inkjet printing (Wijshoff, 2018; Lohse, 2022), food processing industries (Andrade et al., 2013), spray technologies (Breitenbach et al., 2017, 2018) to mention a few. The landscape of outcomes from any drop impact scenario are enormous. The variety of impact outcomes depends on the impact geometry, substrate type (solid vs. liquid), ambient pressure, and surface temperature (Yarin, 2006; Xu et al., 2005; Roy et al., 2019, 2023). A general drop impact phenomenon can be characterized experimentally by a set of non-dimensional numbers like Weber number, Reynolds number, Ohnesorge number, Bond number, Capillary number, and Mach number, to name a few. Each of these non-dimensional numbers represents the balance of dominant scales related to a particular phenomenon observed during impact events. Using appropriate non-dimensional quantities various kinds of impact events like solid impact, drop impact, jet impact (Roy & Basu, 2022) to name a few could be understood in a unified framework (Yarin et al., 2017). Such unified description can explain a wide range of impact phenomena ranging from meteor impacts to granular matter impacts to jet impacts on solids and thin films like bio-mechanical systems (Dowling & Dowling, 2013; Roy et al., 2021; Gilet & Bourouiba, 2015; Wang & Bourouiba, 2018; Bourouiba, 2021). Some phenomena related to drop impacts are spreading, splashing, bouncing, corona formation, and various kinds of Worthington jets. These phenomena are very well studied experimentally Yarin et al. (2017); Thoroddsen et al. (2012) especially spreading and splashing due to their wide application in industries. Recently, scientists are also studying impact phenomena from infectious disease and pathogen transmission perspective (Shetty et al., 2020; Roy et al., 2021) (Hariharan et al., 2022). Drop impact generally is multi-scale, ranging from nano-scale (Li et al., 2015; De Ruiter et al., 2015) to millimeter scales. Despite immense research in drop impact studies for the last one and a half-century, several aspects of drop impact physics still remain elusive and need further investigation.

One such area is the air layer dynamics beneath an impacting drop. When a drop approaches an impacting interface, the bottom part (leading edge) of the drop experiences a rise in pressure due to lubrication pressure buildup. The lubrication pressure, in general, is inversely proportional to the distance between the leading edge and the impacting surface. The increased pressure at the leading edge of the drop causes a deformation of leading edge of the drop, entrapping air sandwiched between the the drop and the impinging interface/substrate. As the air between the drop and substrate thins out various interactions like capillary and van der Waals becomes dominant. The stability of the air layer between the drop and the impacting surface depends on the interaction between capillary and van der Waals interaction. For high speed impact, inertial and compressible effects will also be important (Lesser & Field, 1983) along with the capillary and van der Waals interaction effecting the air layer dynamics. Recent experimental works have shown that drops in the size range of millimeters are supported by micro and nano air films (Kolinski et al., 2012; Langley et al., 2017; Langley & Thoroddsen, 2019). The destabilization of the air layer is inherently related to bubble formation in various drop impact systems. Bubbles in industrial and various engineering impact systems like printing, painting, and cooling are often detrimental to the optimum working efficiency of the system. However, bubbles could also be helpful, as is seen in various natural systems like raindrops impacting liquid bodies. In these natural systems, bubbles help in gaseous exchange between the atmosphere and the liquid bodies supporting aquatic/marine life forms (Woolf et al., 2007).

On approaching the impacting surface, the drop experiences a pressure rise at the leading edge of the drop due to air drainage (outward radial flow) in between the drop and the substrate. Based on the relative air layer thickness, the entrapped air layer can be subdivided into two distinct disjoint regions; the central dimple and the peripheral air disc region (Roy et al., 2022). The dimple forms beneath the impacting drop when the lubrication pressure exceeds the capillary pressure. This is in contrast to solid particle impact at low to moderate velocities, where deformation of the impacting solid object is negligible since the particle behaves as a rigid body. The actual contact between the impacting objects occur, when the dimensions of the entrapped draining air film reaches and becomes smaller than the mean free path of non continuum molecular regimes. Non-continuum analysis in lubrication drainage flows was incorporated by Hocking (Hocking, 1973). Hocking modified the standard lubrication analysis by using the Maxwell slip boundary condition instead of the no-slip boundary condition. Using Hockings analysis as a base, rigid particle impacts like the squeeze flow problem (between rigid spheres or a sphere and a flat plate) typically found in gas-solid suspension, were studied in extensive detail (Davis, 1984). The monodisperse gas-solid suspension was studied mathematically using classical and non-continuum lubrication flows by many scientists (Koch, 1990; Sundararajakumar & Koch, 1996; How et al., 2021). In this work we use some of the asymptotic methods and concepts developed in the context of entrapped air layer dynamics during solid particle impact. The central dimple that forms in drop impact systems has been explored by scientists in the last decade experimentally and computationally (Hicks & Purvis, 2010, 2011; de Ruiter et al., 2012; De Ruiter et al., 2015; de Goede et al., 2019; Driscoll & Nagel, 2011). However, surrounding the central air dimple exists a region that we refer to as the peripheral air disc (Roy et al., 2022). The air layer dynamics in the peripheral region remain relatively unexplored, and scientific literature is sparse. Some of the important works in the context of drop impact exploring the central dimple experimentally and computationally are van der Veen et al. (Van Der Veen et al., 2012), Bouwhuis et al. (Bouwhuis et al., 2012), Mandre et al (Mandre et al., 2009; Mandre & Brenner, 2012). The work by van der Veen et al. (Van Der Veen et al., 2012), is an experimental study that measures the air layer layer profile under impacting drop using high speed color interferometry. The major focus was to measure the time evolution of the air layer profile before and during the contact process with the substrate. Further, from the air layer profile, the velocity of the air exiting the gap between the drop and the substrate during the wetting and bubble formation process was calculated. Bouwhuis et al. (Bouwhuis et al., 2012), explored the competing mechanisms that contols and minimizes the entrained air bubble size during drop impacts on solid surfaces. For large impact size and velocity, the bubble is flattened due to the inertia of the droplet. For smaller size and velocity, capillary effects minimizes the droplet size. Bouwhuis et al. (Bouwhuis et al., 2012) therefore showed that the entrained bubble volume can have a maximum size under the effect of inertial and capillary action. Mandre et al. Mandre et al. (2009) studied the splashing dynamics during drop impact on a solid surface and shows that the initial spreading of the droplet occurs on top of a thin air layer that prevents the initial wetting of the droplet and the substrate. During the initial spreading phase, the drop-air bottom interface develops regions of very high curvature which further generates capillary waves. Another study by Mandre et al. (Mandre & Brenner, 2012) provides an ab-initio description of sheet formation during drop impacts on dry solid surfaces. Further, the major part of the study by Mandre et al. was about the liquid phase and the prediction of sheet formation as a function of experimental parameters like drop size and velocity. The sheet formation dynamics also provides a semi-quantitative method to study various processes related to the sheet formation dynamics and understand various parameter dependence like viscosity, gas pressure, surface roughness, and surface tension. Josserand et al. Josserand & Thoroddsen (2016) wrote an excellent review article that provides a state of the art research for drop impacts on solid surfaces. The article also discusses some of the work related to the measurement of the central air dimple profile. It is to be noted that the above mentioned works does not discuss anything related to non-isothermal impacts, Knudsen field, structures in the peripheral disc, which are some of the major contribution of the current work. In this work, we study the entire air layer profile beneath the impacting droplet on heated substrate experimentally and theoretically under a unified framework using high speed reflection interferometry imaging, asymptotic methods and scaling analysis. Further, the present study also highlights the need of non-continuum hydrodynamics to properly understand the complete picture. We show experimentally how the air layer profile consists of regions that encompasses the entire regime from continuum to non-continuum hydrodynamics beneath the impacting droplet impinging on heated substrate. Recently Chubynsky et al. (Chubynsky et al., 2020) shed some light on the effect of van der Waals force related to gas kinetic effects (GKE) in drop impacts for low Weber number computationally.

In our previous work (Roy et al., 2022), we have explored the peripheral air disc in the context of drop impact on immiscible liquid pools. We discovered that new structure could exist in the peripheral disc due to balance of capillary and van der Walls (vdW) interaction. This work focuses on the air layer dynamics during drop impact on heated substrates at low impact energy/Weber number below the Leidenfrost temperature. We study the dynamic evolution of the central air dimple and the peripheral air disc for different substrate temperatures ranging from $300$K to $473$K. The impact Weber number ($We={\rho}_{l}{V_{0}}^{2}R_{0}/{\sigma}_{aw}$) based on drop radius $R_{0}$, liquid density ${\rho}_{l}$, impact velocity $V_{0}$, and air-water surface tension ${\sigma}_{aw}$ was kept low ($We{\sim}1$) so that the effect of inertia and compressibility could be neglected. At low Weber number regime, lubrication approximation holds good. We have detected non-isotropic structures in the peripheral air disc using high-speed reflection interferometric imaging. The structures in the peripheral air disc region cause asymmetric wetting between the drop and the substrate. The film mode and kink mode are the two fundamental contact modes designated by some recent works (de Ruiter et al., 2015b, a; Li et al., 2015; Langley & Thoroddsen, 2019). The film mode, or the first kink mode, occurs just outside the dimple, whereas the kink mode, or the second kink mode, occurs at the maximum extension of the drop (Chubynsky et al., 2020). The film mode occurs for $We>>1$ whereas the kink mode occurs for $We{\sim}1$. Based on our previous (Roy et al., 2022) and current work, we propose that the underlying cause of both the film and kink mode are the non-isotropic structures in the peripheral air disc region.

As the drop approaches the heated impacting surface (refer to Fig. 1 for the experimental impact configuration), i.e., the gap between the bottom most point of the drop and the surface decreases (Fig. 2(a)), the lubrication pressure increases (since the lubrication pressure is inversely proportional to the gap between the drop and the surface). The lubrication pressure when exceeds the local capillary pressure at the bottom most point of the droplet results in the formation of a dimple. The air-water interface at the bottom of the drop therefore is essentially a curved surface entrapping a thin layer of air having thickness ranging from $h{\sim}5-10{\mu}$m at the center to $h{\sim}70-100$nm at the edge. The air layer profile beneath the drop hence, in general can be thought of as consisting of two disjoint regions based on the relative air layer thickness $h$, the central air dimple, $h{\sim}5-10{\mu}$m(refer to ’*’ in Fig. 2 (b) depicting the central air dimple) and the surrounding region around the central dimple that the authors refer to as the peripheral air disc, $h{\sim}70-100$nm (refer to ’#’ in Fig. 2 (b) and 2(c) for a schematic of the peripheral air disc and a 3D non-dimensional height profile respectively). The central air dimple has been generally studied for iso-thermal drop impact conditions both experimentally and theoretically. It is important to note that non-isothermal drop impact has also been studied extensively by the drop impact community. The major focus though were generally the experimental characterization of impact processes in various regimes like spreading, recceding, leidenfrost, bouncing, and atomization for various non dimensional parameters like Weber number, Capillary number, Reynolds number, Ohnesorge number, Stokes number, viscosity ratio and density ratio between the drop and the surrounding gas. However, there are relatively few studies that explores the central air dimple under non iso-thermal drop impact conditions using high speed interferometry (Qi et al., 2020). It is important to note that the central air dimple thickness is at least one order of magnitude larger than the peripheral air disc thickness. The peripheral air disc region thus is very sparsely studied and still remains elusive. The peripheral disc region is important as the first point of contact generally occurs in this region and is essential for an in-depth quantitative understanding of the bubble formation process. The first point of contact in drop impact scenarios were previously reported to occur majorly in two modes, the first kink mode or film mode for large impact Weber numbers and second kink mode or kink mode for Weber number of the order of unity. In this work we generalize based on the current and our previous impact experiments on liquid and solids that the first and second kink mode are a special case of various structures that are present in the peripheral air disc region as a result of the asymptotic balance between capillary and van der Waals forces. These structure acts as nucleation point of air film destabilization and form rupture points that subsequently leads to air film dewetting and bubble formations. Fig. 2(d) shows a typical Knudsen field beneath an impacting droplet. Two distinct regions are clearly visible, the central dimple having relatively small value of Knudsen number ($Kn<0.01$) can be understood using continuum hydrodynamics. On the contrary the peripheral air disc region is not smooth and has structures of various length scales. The Knudsen number in the disc region is relatively large ($0.01<Kn<0.1$) and hence falls within the purview of slip flow and non-continuum hydrodynamics regime where the first point of contact occurs.

The present study explores the air layer dynamics at low-impact Weber number regime ($We{\sim}1$) on heated glass substrate below Leidenfrost temperature using high speed reflection interferometry imaging and theoretical analysis. The main theme of this work is to understand the mechanics of the central air dimple and its immediate neighbourhood (the peripheral air disc). The detailed characterization for various experimental parameter ranges is outside the scope of the present work.

2 Materials and Methods

We study liquid water drop impact on heated solid glass substrate below the Leidenfrost temperature. The drop on approaching the impacting surface causes lubrication pressure rise in the thin entrapped air layer between the drop and the substrate. The lubrication pressure on exceeding the capillary pressure across the drop interface results in the formation of the central air dimple and the peripheral air disc. We study evolution of the central air dimple and the peripheral air disc using high speed imaging and theoretical analysis.

2.1 Experimental Set up

Using high-speed reflection interferometric technique, we study air layer dynamics beneath the impacting drop on a glass substrate at various surface temperatures. The schematic and assembly of the actual experimental setup is shown in Fig. 1(a) and 1(b), respectively. De-ionized water drop of radius $R_{0}=1.1\pm 0.1{\>}$mm generated from a syringe pump (New Era Pump Systems, NE-1010) impinges a glass substrate of dimensions $75\times 25\times 1$mm3 (procured from Blue Star) with an impact velocity of $V_{0}=0.25$m/s ($We{\sim}1$). The corresponding impact energy of the drop is $0.3{\mu}$J. The drop was allowed to fall freely under gravity from a designated height to achieve the desired impact velocity. We used cleaned glass slides for the experiments. The glass slides were sonicated in a bath of isopropanol for 10 minutes. The glass slide was then placed on the annular-shaped heater (copper surface) connected to a PID controller for maintaining the precise temperature conditions on the heated surface. The surface temperature of the glass plate at the impact region (placed on the annular heater) was measured using a K-type thermocouple within the error range of $\pm 5^{\circ}$C=$\pm 5$K. We capture the drop impact phenomena using two simultaneous high speed imaging methods (shadowgraphy and interferometry) at identical frame rates to understand and characterize the impact and air-layer dynamics for substrate temperatures $300$K, $353$K, $423$K and $473$K. Shadowgraphy is performed from the side view, and reflection interference imaging from the bottom view as shown in Fig. 1(a) and Fig. 1(b). For the shadowgraphy we use a LED (Light Emitting Diode) light source for backlight and a high speed camera connected with a zoom lens (Navitar). The shadowgraphy and the interferometry imaging were both performed at 50000 FPS (frames per second). Sample shadowgraphy and interference image sequence are shown in the top and bottom image sequence panel of Fig. 1(c) respectively for substrate temperature $T_{s}=300$K. The drop impacts the glass with a velocity $V_{0}=0.25$m/s. The dotted rectangle in the side view shows the field of view for the bottom view interferograms. The timestamps are in microseconds for both the shadowgraphy and interference image sequence. The scale bar shown in red for the side and interferometric views represent $1.3$mm and $0.296$mm respectively. The reflection interference imaging setup was assembled on the bottom side of the glass plate to record the dynamic interference patterns during drop impact. The interference system consists of a pulse-diode laser source (Cavitar Cavilux smart UHS, 400 W power, ${\lambda}_{s}=640$nm) for illumination purpose; a beam-splitter; a 4x microscope objective; a dichoric mirror; and a high speed camera as shown in Fig. 1(a). The reflected light beam from the beam-splitter passes through the microscope objective onto the glass surface. The reflected light from the top glass surface and the bottom surface of the air-water interface of the impacting drop superimposes to produce an interference pattern. The interference pattern was focused using a zoom lens (Navitar) with the assistance of in-line dichoric mirror, as shown in the schematic of Fig. 1(a) and the actual experimental setup Fig. 1(b). The dynamic fringe patterns were recorded using a high-speed camera (Photron SA5) with similar frame rates to the side-view imaging (i.e., at 50000 FPS). A pixel resolution of $1.5{\mu}$m per pixel was used to resolve the constructive and destructive interference patterns.

2.2 Interferometry post-processing

A monochromatic beam from the diode laser source passes from the bottom to the top through the microscope objective and gets transmitted through the glass substrate (Fig. 1(a)). Some part of the transmitted beam suffers reflection at the air-glass interface (top surface of the glass substrate) and travels back to the CCD of the high-speed camera with an intensity $I_{1}(x,y)$ where $(x,y)$ represents the image coordinate. The other part of the transmitted beam coming out of the glass substrate travels through the air layer beneath the impacting drop and undergoes reflection at the air-water drop interface (bottom surface of the impacting drop) (refer to Fig. 2(a) and 2(b)). The reflected beam from the air-water drop interface travels back to the CCD of the high-speed camera with an intensity of $I_{2}(x,y)$. Beams with intensity $I_{1}(x,y)$ and $I_{2}(x,y)$ superimpose to produce interference pattern recorded by the high-speed camera. The interference occurs, because the beams $I_{1}(x,y)$ and $I_{2}(x,y)$ are coherent. Further, the beam with intensity $I_{2}(x,y)$ is phase shifted to the beam with intensity $I_{1}(x,y)$ due to the additional path traversed through the air layer. The path difference between beams $I_{1}(x,y)$ and $I_{2}(x,y)$ causes a phase shift producing a resultant interference pattern of intensity $I(x,y)$ given by (Daniel et al., 2017; Limozin & Sengupta, 2009; Sugiyama et al., 2006)

[TABLE]

where ${\phi}(x,y)$ is the phase difference between beams having intensity $I_{1}(x,y)$ and $I_{2}$(x,y). $I(x,y)$ is the resultant intensity of the interference pattern recorded by the camera at an image coordinate $(x,y)$. Sample dynamic side views image sequence and interferograms are shown in Fig. 1(c) during drop impact on $300$K glass substrate at $We{\sim}1$. The timestamps are in microseconds. The main feature of equation (2.1) is to represent the resultant intensity distribution on the image coordinate plane due to interference of beams having intensity $I_{1}(x,y)$, $I_{2}(x,y)$ and the cosine-modulated phase field given by $cos{\phi}(x,y)$. If the phase gradient is high, high-density fringes are observed. Similarly, a smaller phase gradient corresponds to low fringe density. The dark fringes correspond to destructive interference, whereas the bright fringes correspond to constructive interference. From equation (2.1), we can observe that a general resultant intensity field could be represented mathematically as

[TABLE]

where $B(x,y)=I_{1}(x,y)+I_{2}(x,y)$ and $A(x,y)=I_{1}(x,y)I_{2}(x,y)$. Physically $A(x,y)$ and $B(x,y)$ represents the fringe amplitude and background intensity fields respectively. On proper background subtraction and fringe amplitude normalization, the resultant intensity field $F(x,y)$ could be represented as

[TABLE]

Therefore from an input experimental image $I(x,y)$, a background subtracted normalized image $F(x,y)$ could be computed using equation (2.3). The unknown phase field ${\phi}(x,y)$ could be computed from the inverse cosine function, i.e.,

[TABLE]

The inverse cosine function being multi-valued, the true phase ${\phi}(x,y)$ could be determined using a proper signature of ${\phi}_{0}(x,y)$ if represented in a wrapped form ($\phi(x,y)=\pm{\phi}_{0}(x,y)$). Given that $-1{\leq}F(x,y){\leq}1$, ${\forall}$ $(x,y)\in[0,x_{max}]{\times}[0,y_{max}]$, i.e. for all $(x,y)$ coordinates in the image plane, $\phi_{0}(x,y)\in[0,\pi]$ for single valuedness of the inverse cosine function. Here the pair $(x_{max},y_{max})$ represents the region of interest boundary coordinates diagonally opposite to $(0,0)$ which is chosen as the origin in the image plane. For the entire image captured from the CCD sensor, the pair $(x_{max},y_{max})$ could be mapped to the image resolution in appropriate length units. In general, frequency guided sequential demodulation (FSD) methods (Kemao & Soon, 2007; Wang & Kemao, 2009), the actual phase field is calculated using linear and quadratic approximation methods, where the gradient phase fields (also known as local frequencies) are calculated by optimization an appropriate cost function using least square methods. A typical cost function that can be used for the optimization process and computing the gradient fields $({\omega}_{x}(x,y),{\omega}_{y}(x,y))$ is given by

[TABLE]

where $N_{xy}$ denotes a neighborhood around the point $(x,y)$. Assuming linear approximation to compute ${\phi}(x,y;{\delta}_{1},{\delta}_{2})$ in the neighborhood $N_{xy}$ from ${\phi}_{0}(x,y)$, we have using the linear Taylor’s series approximation in two dimensions

[TABLE]

The proper sign of the phase field is determined by allowing the gradients to be continuous. The optimization is generally done using exhaustive search algorithms and more recently using Levenberg-Marquardt optimization methods to improve the computational time (Wang & Kemao, 2009). In this work, the dynamic fringe patterns during drop impingement have been processed using the fast-frequency guided sequential demodulation (FFSD) method (Kai & Kemao, 2010; Wang & Kemao, 2009; Roy et al., 2022). Fast-frequency guided algorithms are suitable for transient phenomena (Roy et al., 2022) and consequent comprehensive extraction of the two-dimensional phase field of the complex fringe patterns. The FFSD algorithm has been used to calculate the signature of the unknown phase. Using equation (2.3) we can construct a scalar field that will be used for signature determination as

[TABLE]

$\bm{{\nabla}}F(x,y)$ can be computed from the field $F(x,y)$ and $\bm{{\nabla}}{\phi}(x,y)$ can be computed using ${\phi}_{0}(x,y)$ from equation (2.4) in the neighborhood $N_{xy}$ around the point (x,y). From equation (2.7) we can see that the sign(${\phi}(x,y)$)=-sign($\bm{{\nabla}}F(x,y){\cdot}\bm{{\nabla}}{\phi}(x,y))$ since $|\bm{{\nabla}}{\phi}(x,y)|^{2}{\geq}0$ and sign($sin({\phi}(x,y))$)=sign(${\phi}(x,y)$). Therefore, to determine the signature of the $\phi(x,y)$, we have to determine the signature of $\bm{{\nabla}}{\phi}(x,y)$. The gradient of phase distribution for the obtained fringe patterns could be represented in terms of complex numbers for efficient computing as

[TABLE]

Here, ${\omega}(x,y)$ signifies the total amplitude of the phase gradient field vector and ${\theta}(x,y)$ represents the direction of the vector field at any point $(x,y)$. The complex field from equation (2.8) is smoothed by averaging the square of the gradient phase field values over some neighborhood region $N_{xy}$ to obtain ${\omega}_{s}(x,y)$ and ${\theta}_{s}(x,y)$. The suffix ’s’ denotes the smoothed field. Using smoothed version of the phase gradient field, the phase gradient field could be computed without any sign ambiguity by requiring spatial continuity on the gradient of the phase field, as $\bm{\nabla}{\phi}(x,y)=\bm{\nabla}{\phi}_{0s}(x,y)$ if $\bm{\nabla}{\phi}(x,y){\cdot}\bm{\nabla}{\phi}_{0s}(x_{*},y_{*}){\geq}0$ where $(x_{*},y_{*})$ denotes the adjacent pixels of $(x,y)$. Similarly, $\bm{\nabla}{\phi}(x,y)=-\bm{\nabla}{\phi}_{0s}(x,y)$ if $\bm{\nabla}{\phi}(x,y){\cdot}\bm{\nabla}{\phi}_{0s}(x_{*},y_{*}){<}0$. Using the proper sign of the gradient of the phase field in equation (2.7), the phase field signature can be determined uniquely without any ambiguity. However, the phase field can be ambiguous to a factor of $2k{\pi}$ in general (the period of cosine modulated function), where $k{\in}Z_{I}$, $Z_{I}$ represents set of all integers. To circumnavigate the phase ambiguity issue, the obtained phase are constructed to be continuous and unwrapped. Using the unwrapped phase field ${\phi}(x,y)$, the two dimensional thickness profile $h(x,y)$ could be computed. The phase ambiguity can be handled using various interferometry techniques like color interferometry, dual/multi wavelength interferometry (Josserand & Thoroddsen, 2016). Further, the thickness profile can also be measured without a $2{\pi}$ phase ambiguity using monochromatic interferometry. The fringes obtained from monochromatic and coherent source of illumination is known as Haidinger fringes Hecht (2012) in general. Some of the commonly used techniques to determine the thickness without a phase ambiguity for the monochromatic illumination are the peak and valley methods Gillen & Guha (2005); Shukla et al. (2006), reference thickness methods Bouwhuis et al. (2012), Haidinger quadrature Kim et al. (2017); Park et al. (2019) to name a few. We use the reference thickness methods in our measurements to determine the absolute thickness profile, since the absolute thickness are known at certain positions of the interference images. One such reference point is the hole nucleation spots (air layer rupture locations) where the absolute thickness monotonically approaches zero in time (refer to ’’ in Fig. 1 of the supplementary document) and hence any arbitrary integer multiple of wavelength will not be able to match the actual thickness profile evolution close to the rupture (nucleation) locations. Therefore temporal causality for the evolution of the thickness profiles at the hole nucleation point can lead to the determination of absolute thickness from a monochromatic illumination. The evolution of air layer thickness in general is a causal process in time and hence allow the computation of actual air layer thickness based on the rupture locations, where the thickness profile goes to zero provided the phase field is continuous. Further, based on the observations in the current experiment, that most rupture points are very close to the outer expanding rim, the outer rim location are used to set up a zero thickness reference for very early time data when rupture locations have not formed (refer to ’’ in Fig. 1 of the supplementary document). It is important to observe that the average air layer thickness decreases with radial coordinate upto the outer expanding rim of the air disc. The average thickness at the outer expanding rim therefore will be smallest in the the entire air layer region well beyond the measurement range and hence can be considered as reference points of zero thickness.

The height of the air-layer profiles can then be estimated for the preferred wavelength of the light source ${\lambda}_{s}$ and refractive index $n$ of the fluid medium as (Kitagawa, 2013):

[TABLE]

The images acquired with the shadowgraphy technique (i.e., side-view imaging) were processed using a combination of open-source tools like ImageJ (Schneider et al., 2012) and python (Van Rossum & Drake, 2009) for one to one comparison for variation of drop shape during impact and simultaneous change in interference pattern for varying surface temperatures. Sufficient number of experimental runs were conducted to ensure statistically significant datasets.

3 Results and Discussions

3.1 Global overview

Fig. 2(a) depicts the schematic and the geometry of the impact condition under study. The air layer dynamics beneath an impacting drop on a heated substrate at various substrate temperatures ($T_{s}=300$K, $353$K, $423$K and $473$K) has been studied using high-speed reflection interferometry imaging at low impact kinetic energy ($We{\sim}1$). The air layer beneath the drop is subdivided into two distinct disjoint regions that we refer to as the central air dimple, surrounded by the peripheral air disc (Fig. 2(b)). Fig. 2(b) depicts the air layer profile schematically. The * denotes the central air dimple, and # represents the peripheral air disc. The thickness of the central dimple, in general, is larger than the peripheral air disc region by one order of magnitude (i.e., $h_{*}/h_{\#}{\sim}\mathcal{O}(10)$). We observe using high-speed interferometry imaging that the geometrical features features of the dimple (i.e., central dimple maximum thickness and central dimple diameter) are weakly dependent on the substrate temperature $T_{s}$ and depends on the impact energy/Weber number ($We$) in general. The average maximum height of the central air dimple is ${\sim}\mathcal{O}$($6{\mu}$m). We also observe that the air layer rupture time scale and rupture radius increases with increase in the substrate temperature. The rupture time scale lies within the range of ${\sim}\mathcal{O}$($100{\mu}$s) to ${\sim}\mathcal{O}$($400-500{\mu}$s) for a substrate temperature variation between $300$K to $473$K respectively. Similarly, the rupture radius variation from $300$K to $473$K substrate temperature lies within the range ${\sim}\mathcal{O}$($200{\mu}$m) to ${\sim}\mathcal{O}$($600{\mu}$m) respectively. Further, we show that a Gaussian profile can approximate the shape of the central air dimple during its temporal evolution. A standard continuum Stokes flow regime governs the central dimple evolution. In contrast, the dynamics in the peripheral disc region falls in the non-continuum regime characterized by a relatively high value of Knudsen number (Fig. 2(c) and 2(d)). Fig. 2(c) shows a 3d wireframe plot of the non-dimensional height profile of air layer at a surface temperature of $T_{s}=423$K normalized with respect to the mean free path of air. The non-dimensional height profile is a one-to-one map of the inverse Knudsen field map. The dotted red rectangle in Fig. 2(b) and Fig. 2(c) depicts the interfacial perturbations in the peripheral air disc. Fig. 2(d) depicts the Knudsen contour field corresponding to Fig. 2(c). We observe intriguing 3d structures along the peripheral air disc region that acts as a seed for asymmetric wetting by the drop on the substrate. Such structural length scales exist at short time scales due to balance of capillary and molecular interactions like van der Waals. The structures in the peripheral air disc region lead to asymmetric wetting between the drop and substrate. The structures in the peripheral disc are the underlying cause of the kink and film mode of first contact between the drop and substrate. For the impact Weber number in this study ($We{\sim}1$), we observe kink mode of contact in contrast to the film mode of contact observed in our previous study on liquid pools ($We{\sim}10$) (Roy et al., 2022). This is in accordance with the numerical findings of Chubynsky et al. (Chubynsky et al., 2020).

The remaining article is arranged as follows. From section 3.2 onwards, we analyze the air layer dynamics using experimental and theoretical methods. We subdivide the analysis into various sub-sections. Section 3.2 introduces the coordinate system and the initial condition before dimple formation. Section 3.3 focuses on the asymptotic analysis for the velocity and pressure field in the dimple region. Section 3.4 provides an estimate of the dynamic dimple radius. Section 3.5 deals with the radial expansion of the peripheral air disc. We discuss and explore the air layer dynamics and perturbative structures in the peripheral air disc in section 3.6.

We begin the analysis by considering the geometry and initial condition for understanding air layer dynamics underneath the impacting drop.

3.2 Coordinate system and Air profile just before dimple formation

Fig. 2(a) depicts a schematic representation of the impact geometry just before the central air dimple begins to form. We use a cylindrical coordinate system to pose the mathematical problem.

From the geometry of the drop, using elementary trigonometric correlations (Fig. 2(a)) we have

[TABLE]

where ${\phi}$ is an angle parametrizing the bottom surface of the drop, $r$ is the radial coordinate from the axis of the cylindrical coordinate system (vertical $z$ axis), and $R_{0}$ is the impacting drop radius. The build of lubrication pressure in the air layer just before the dimple formation is shown schematically as a red halo in Fig. 2(a). The pressure in the liquid drop and the air layer is $p_{l}$ and $p_{a}$ respectively. The closest approach distance between the drop and the substrate is denoted by $h_{0}={\epsilon}R_{0}$, where ${\epsilon}$ is a small asymptotic parameter that we will be used later in computing the pressure and velocity field respectively. The bottom part of the impacting drop just before the dimple formation can be approximated as part of a sphere (Fig. 2(a)). The initial air layer thickness profile therefore can be written as

[TABLE]

where $z(r)$ is the air thickness at any radius $r$ between the drop and the substrate. $h_{0}$ is the minimum air layer thickness between the drop and the substrate before dimple formation begins to occur at $r=0$. Being part of a sphere, equation (3.1) is symmetric about the vertical $z$ axis in a cylindrical coordinate system. Fig. 3(a) visualizes the drop profile along a sectional mid-plane passing through the drop centre. The white solid 1D curve in Fig. 3(a) represents the drop shape in the non-dimensional cartesian coordinate space.

3.3 Asymptotic analysis

The asymptotic analysis in this section is based on the concepts developed in the context of solid particle impacts in gas flows (Koch, 1990; Sundararajakumar & Koch, 1996; How et al., 2021). We follow the analysis of How et al. (How et al., 2021) closely. We begin with the asymptotic analysis of the air layer beneath the impacting drop. The asymptotic analysis analyzes the pressure and velocity field of the air layer just before the dimple starts to form. Equation (3.1) could be rewritten for $r<<R_{0}$ using binomial expansion (region of interest is very close to the bottom part of the drop) as

[TABLE]

where ${\epsilon}=h_{0}/R_{0}$ is the asymptotic parameter characterizing the initial condition. Equation (3.2) depicts a parabolic (quadratic) profile approximating the bottom surface of the impacting drop. Normalizing the air layer thickness $z(r)$ with respect to minimum thickness of the air layer just before dimple formation $h_{0}={\epsilon}R_{0}$, equation (3.2) could be written in terms of a non dimensional thickness profile $\bar{z}$

[TABLE]

where $\bar{z}=h(\bar{r})$ represents the non-dimensional height profile of the drop surface just before the dimple formation process. Defining a non-dimensional radial coordinate $\bar{r}=r/(R_{0}{\epsilon}^{1/2})$ equation (3.3) could be written as

[TABLE]

The white solid curve in Fig. 3(a) depicts equation (3.4) graphically along the sectional mid plane of the drop in cartesian coordinates $X-Z$ plane. The continuity equation in cylindrical axi-symmetric coordinates can be written as

[TABLE]

Using $z=\bar{z}R_{0}{\epsilon}$, ${r}=\bar{r}R_{0}{\epsilon}^{1/2}$, $u_{z}=\bar{u}_{z}V_{0}$ in equation (3.5) we have

[TABLE]

where $V_{0}$ is the impact velocity of the drop. Defining $\bar{u}_{r}=u_{r}/(V_{0}/{\epsilon}^{1/2})$ the dimensionless continuity equation becomes

[TABLE]

The radial momentum equation in cylindrical coordinates can be written as

[TABLE]

Using the scalings for $z=\bar{z}R_{0}{\epsilon}$, ${r}=\bar{r}R_{0}{\epsilon}^{1/2}$, $u_{z}=\bar{u}_{z}V_{0}$, $u_{r}=\bar{u}_{r}V_{0}{\epsilon}^{-1/2}$, $t=\bar{t}R_{0}{\epsilon}/V_{0}$, $p=\bar{p}{\mu}_{a}V_{0}/(R_{0}{\epsilon}^{2})$ in equation (3.8) the radial momentum equation becomes

[TABLE]

where $Re_{\epsilon}={\rho}R_{0}V_{0}{\epsilon}/{\mu}_{a}=Re{\>}{\epsilon}$. $Re$ is the standard Reynolds number defined with respect to drop radius $R_{0}$, whereas $Re_{\epsilon}$ is the modified Reynolds number with respect to $R_{0}$. In the limit of ${\epsilon}{\rightarrow}0$ equation (3.9) reduces to

[TABLE]

Equation (3.10) is the governing flow field equation near the leading edge of the drop. Equation (3.10) essentially describe continuum Stokesian flow in a quasi-steady framework (Sundararajakumar & Koch, 1996; How et al., 2021). The $z$ momentum equation in cylindrical coordinates can be written as

[TABLE]

Using the scalings for $z=\bar{z}R_{0}{\epsilon}$, ${r}=\bar{r}R_{0}{\epsilon}^{1/2}$, $u_{z}=\bar{u}_{z}V_{0}$, $u_{r}=\bar{u}_{r}V_{0}{\epsilon}^{-1/2}$, $t=\bar{t}R_{0}{\epsilon}/V_{0}$, $p=\bar{p}{\mu}_{a}V_{0}/(R_{0}{\epsilon}^{2})$ in equation (3.11) the $z$ momentum equation becomes

[TABLE]

In the limit of ${\epsilon}{\rightarrow}0$ equation (3.12) reduces to

[TABLE]

Integrating equation (3.10) with respect to $\bar{z}$ twice we have

[TABLE]

The quantities $C_{1}$ and $C_{2}$ are evaluated from the respective boundary conditions. The continuum boundary conditions for the velocity field components are

[TABLE]

[TABLE]

Evaluating $C_{1}$ and $C_{2}$ from equation (3.15) and equation (3.16); equation (3.14) for the radial velocity component reduces to

[TABLE]

The pressure gradient term ${{\partial}\bar{p}}/{{\partial}\bar{r}}$ is still unknown at this stage and is computed from the continuity equation (3.7). Integrating equation (3.7) with respect $\bar{z}$ we have

[TABLE]

Using the limits in the second integral we have

[TABLE]

Applying the boundary conditions from equation (3.15) and (3.16) in equation (3.19), we have

[TABLE]

[TABLE]

Substituting equation (3.17) in equation (3.21) we have

[TABLE]

[TABLE]

At $\bar{r}=0$, ${{\partial}\bar{p}}/{{\partial}\bar{r}}=0$ (symmetry boundary condition). Therefore $D_{1}=0$

[TABLE]

Substituting the radial pressure gradient from equation (3.24) in equation (3.17), the radial velocity could be written as

[TABLE]

Using the undeformed asymptotic drop shape profile $h(\bar{r})=1+{\bar{r}}^{2}/2$ and integrating equation (3.24) we have

[TABLE]

Further, substituting the dimple profile in equation (3.25) the radial velocity becomes

[TABLE]

Fig. 3 depicts the asymptotic normalized pressure and velocity fields in the air layer beneath the drop just prior to dimple formation. Fig. 3(a) represents the axisymmetric pressure field in the $\bar{X}-\bar{Z}$ cartesian plane according to equation (3.26). The respective cartesian coordinates $x$, $y$, $z$ are expressed as non-dimensional cartesian coordinates using the scales ${x}=\bar{X}R_{0}{\epsilon}^{1/2}$, ${y}=\bar{Y}R_{0}{\epsilon}^{1/2}$, $z=\bar{Z}R_{0}{\epsilon}$. We observe from equation (3.26) that the non-dimensional pressure is maximum at the origin of coordinate system ($\bar{p}(\bar{r}=0)=3$). The white curve shows the asymptotic air-water drop interface prior to dimple formation begins according to equation (3.4). The actual experimental value of ${\epsilon}$ observed was ${\epsilon}=1.146{\times}10^{-3}$. Fig. 3(b) represents the asymptotic velocity field data in the $\bar{X}$-$\bar{Z}$ plane using equation (3.27) and the linear approximation of the z-component of the velocity field $\bar{u}_{z}$ satisfying the boundary condition given by equation (3.15), (3.16) owing to the very small value of ${\epsilon}$ prior to dimple formation. The black dotted line represents $\bar{Z}=0.5$ horizontal plane parallel to $\bar{X}-\bar{Y}$ plane. Fig. 3(c) represents the axisymmetric pressure field in the $\bar{X}-\bar{Y}$ plane for $0<\bar{Z}<1$. Fig. 3(d) represents the axisymmetric velocity field in a plane given by $\bar{Z}=0.5$ parallel to the $\bar{X}-\bar{Y}$ plane according to equation (3.27). The velocity field denotes radially outward air flow from the center representing the air draining out due to the impacting drop. The dimple starts to form when the air pressure at the center of the dimple becomes greater than the capillary pressure.

[TABLE]

The non dimensional pressure at the dimple centre is $\bar{p}(\bar{r}=0)=3$ from equation (3.26). Using the scale transformation between $p$ and $\bar{p}$ i.e., $p=\bar{p}{\mu}_{a}V_{0}/(R_{0}{\epsilon}^{2})$ we have

[TABLE]

On simplifing further and solving the inequality for ${\epsilon}$ we have

[TABLE]

Therefore the critical value of ${\epsilon}$ before central dimple formation begins is

[TABLE]

where $Ca={\mu}_{a}V_{0}/{\sigma}_{aw}$ is the capillary number.

The actual central dimple thickness profile evolution was measured from high speed interferometric imaging as shown in Fig.4. Fig. 4(a), 4(b), 4(c) and 4(d) represents the central dimple evolution with time for substrate temperatures of $300$K, $353$K, $423$K and $473$K respectively. The time coordinate has been renormalized appropriately for non-dimensional purposes as $\bar{t}=tV_{0}/(R_{0}{\epsilon})$. $\bar{t}=0$ is referenced with respect to the closest approach of the drop to the substrate prior to dimple formation, i.e. the time instant when ${\epsilon}=1.146{\times}10^{-3}$. The shape of the drop close to the origin of coordinates is concave up shape represented by equation (3.4) locally (quadratic profile). Since at $\bar{t}=0$, ${\epsilon}<{\epsilon}_{crit}$, the pressure below the drop is greater than the capillary pressure (refer to equation 3.28-3.31). The excess pressure causes radial air flow outwards as depicted in Fig. 3(b) and 3(d). In addition to the flow as a response to the increase in air pressure, the air-water interface deforms and its curvature changes as depicted in Fig. 4 at different time instants labelled by $\bar{t}$. We discover that the central air dimple evolution profiles can be approximated by a Gaussian curve(refer to black curves in Fig.4) of the form

[TABLE]

where $a$ characterizes the dimple thickness at the origin and ${\sigma}$ characterizes the spread of the profile (the standard deviation in terms of Gaussian), $\bar{r}$ is the radial coordinate.

The black curves in Fig.4 depicts the Gaussian approximation to the central dimple profile evolution at different time instants. The Gaussian approximation is a good estimation of the the central air dimple height profile as can been seen from Fig. 5 for all substrate temperature. Fig. 5 depicts the error percentage $\bar{\epsilon}$ plotted as a function of non dimensional x-coordinate $\bar{X}$ with non dimensional time $\bar{t}$ as a parameter for various surface temperature $T_{s}=300,{\>}353,{\>}423,{\>}473{\>}$K in Fig. 5(a), 5(b), 5(c), 5(d) respectively. $\bar{\epsilon}$ is calculated as $\bar{\epsilon}(\bar{X})=(|h_{exp}(\bar{X})-h_{G}(\bar{X})|{\times}100)/h_{exp}(\bar{X})$, where $h_{exp}(\bar{X})$ and $h_{G}(\bar{X})$ is the air layer thickness and the Gaussian approximation of the air layer height profile at a particular $\bar{X}$ coordinate respectively. In Fig.5, $\bar{X}=-5,5$ is depicted as vertical dashed red lines. The horizontal dashed black line represents a $10\%$ error. It could be observed that the average error percentage (averaged over $\bar{X}$) in approximating the central air dimple is less than $5\%$ for $\bar{X}$ in the range $[-5,5]$ for all substrate temperature except $T_{s}=473$K. The error percentage $\bar{\epsilon}$ increases for $|\bar{X}|>5$, i.e. in the peripheral air disc region. Further, the author also wants to highlight that the error between the Gaussian fit and the experimental value increases with increase in temperature as can be seen for $T_{s}=473$K. However, the average error for $T_{s}=473$K shown in Fig. 5(d) is still smaller than $10\%$. This higher percentage error is due to the flattening nature of the dimple profile at higher temperatures caused due to combined effect of the weight of the droplet and the relatively stable air cushioning effect. Since, the air film evolution time scale increases with temperature, gravity and hence weight of the droplet along with the relatively thicker and stable air cushion can therefore effect the shape of the dimple.

Substituting equation (3.32) in equation (3.24) we have

[TABLE]

Integrating with respect to $\bar{r}$ we have

[TABLE]

[TABLE]

Further on substituting equation (3.32) and (3.33) in equation (3.17) the radial velocity field becomes

[TABLE]

The radial velocity field variation with time is plotted as a contour field in Fig. 6, Fig. 7, and Fig. 8. Fig. 6 depicts the radial velocity field contours (Fig. 6(a)) and velocity profiles (Fig. 6(b)) in the $\bar{Z}$-$\bar{X}$ plane for $T_{s}=300$K Figs. 7(a), 7(b) and 7(c) represents the dimple evolution and its associated radial velocity contour fields at three different time instants $\bar{t}=3.36$, $\bar{t}=6.71$ and $\bar{t}=10.07$ for $T_{s}=300$K. The white dotted curve depicts the air-dimple interface evolution as a function of time according to equation (3.32) where $a$ and ${\sigma}$ are function of time. The time $\bar{t}=10.07$ corresponds to the instant where the dimple is at its maximum height and has attained a steady profile. It is interesting to note that the maximum radial velocity scale decreases as time progresses as can be observed from Fig. 7. This is in accordance with the continuity equation that as time progresses, the velocity profiles becomes elongated in the vertical direction due to increasing dimple height and narrower in the horizontal direction denoting a decrease in velocity scale at a given fixed radial coordinate. Fig. 8 depicts the dynamic dimple evolution and the contour maps for impact on substrate at temperature $T_{s}=353K$. Fig. 8(a), 8(b), 8(c) and 8(d) represents the various non-dimensional time instants represented by $\bar{t}=6.71$, $\bar{t}=10.07$, $\bar{t}=13.42$ and $\bar{t}=20.13$ respectively. Similar radial velocity contours are plotted for other substrate temperatures in Fig. 2 ($423$K) and Fig. 3 ($473$K) of the supplementary document. We can observe that in all substrate temperature the non-dimensional radial velocity scales are of similar order of magnitude indicating the dynamics in the central dimple region is inherently independent of substrate temperature in the range considered in this work. The important idea that these plots reveal is the nature of temporal evolution of the radial velocity field. It should be noted that the asymptotic radial velocity and the its corresponding pressure field in the central dimple evolves as a quasi-steady field with the central dimple air-water interface (shown white dotted lines in the contour field plots) as the moving boundary. The maximum central dimple thickness therefore is weakly dependent on the substrate temperature. The region around the central dimple consists of the peripheral air disc where the air thickness is of the order of $100$nm. Small perturbations develops at the air-water interface in the peripheral air disc region due to capillary and van der Waals interaction (discussed in detail in a later section). This results in asymmetric air layer dewetting/rupture (or it could be called wetting of the drop with the substrate).

As the air between the drop and the substrate drains and thins out, structures appear in the peripheral air disc. The air thickness in the disc region reaches close to mean free path before first point of contact between the drop and the substrate. The first point of contact is designated as the air layer rupture point. The time it takes to form the rupture point is known as the rupture time scale and the radial distance between the rupture location and the origin of coordinates is known as the rupture radius. The rupture time scale increases with increase in substrate temperature as can be observed in the box-plot depicted in Fig. 9(a). We also observe correspondingly in Fig. 9(b) that the rupture radius increases with the substrate temperature $T_{s}$. The green horizontal line in the box plot denotes the median value of the rupture time and rupture radius. The rupture time and rupture radius are measured in ${\mu}$s and ${\mu}$m respectively. Corresponding to ${\epsilon}_{crit}$ from equation (3.31), the critical air layer thickness prior to dimple formation is given as

[TABLE]

where $Ca={\mu}_{a}V_{0}/{\sigma}_{aw}$ is the capillary number and $R_{0}$ is the impacting drop radius. The capillary number depends on the temperature due to surface tension and viscosity dependence of temperature as ${\sigma}_{aw}=a-bT_{s}$ and ${\mu}_{a}={\mu}_{ref}(T_{s}/T_{ref})^{\omega}$ (Valentini et al., 2023) respectively. Here, ${\mu}_{ref}$ is the air viscosity at the reference temperature $T_{ref}$ which is $300K$ for the current study. Incorporating the surface tension and viscosity dependence on temperature, the critical thickness from equation (3.37) becomes

[TABLE]

In general the maximum central dimple thickness is proportional to the critical thickness given by equation (3.38) as the pressure at the center of the dimple that causes deformation of the drop interface is correlated to ${\epsilon}_{crit}$ and hence to $h_{crit}$. Therefore, the maximum dimple thickness at $\bar{r}=0$ has a temperature dependence similar to equation (3.38) i.e. $h(\bar{r}=0)|_{max}{\sim}\sqrt{T^{\omega}_{s}}/\sqrt{a-bT_{s}}$. Fig. 9(c) compares the experimental central dimple maximum thickness with the theoretical scaling of $h(\bar{r}=0)|_{max}$. The experimental values agrees with the theoretical scale within the experimental uncertainty. The lubrication force prior to dimple formation is given by

[TABLE]

Using the value of non-dimensional pressure from equation (3.26) in equation (3.39) we have

[TABLE]

The lubrication force varies inversely as ${\epsilon}$ and hence blows up in the limit ${\epsilon}{\rightarrow}0$. Physically for impact of solid objects the singularity is real and the singularity cutoff occurs through non-continuum effects. For the present case of impact of liquid drop, the singularity at the center of coordinates is avoided due to deforming nature of the air water interface. i.e., for drop impacts on solids or liquids, the singularity is never reached at the center of coordinates due to dimple formation mechanisms discussed above. For impact of drops on solids and liquids at moderate to low impact energies the system approaches the mathematical singularity in the peripheral air disc region. However, physically the singularity is never reached and is cut off by non-continuum 3d structures that occurs due to balance of capillary and van der Waal interactions. Using ${\epsilon}=h_{0}/R_{0}$ in (3.40), the lubrication pressure prior to dimple formation becomes

[TABLE]

Applying Newtons second law provides us the dynamical equation of motion for the impacting drop,

[TABLE]

where $F_{mg}$ is the weight of the drop, $F_{b}$ is the buoyancy force, $F_{lub}$ is the lubrication force, $F_{v}$ is the viscous drag, $m$ is the mass of the drop, $V_{CM}$ centre of mass velocity of the drop and $t$ is time. As the drop approaches the impacting surface just prior to dimple formation, the center of mass velocity of the drop attains approximately a constant value for the time duration under consideration, and hence $dV_{CM}/dt=0$. The non dimensional centre of mass $z$ coordinate $z_{CM}/R_{0}$ is plotted as a function of non dimensional time $t/{\tau}$ in Fig. 4 of the supplementary document. From the figure, we can observe that $V_{CM}$ is constant as ${z}_{CM}$ decreases linearly with time. Therefore, the centre of mass acceleration is zero ($dV_{CM}/dt=0$) just prior to dimple formation. The Reynolds number is given as $Re={\rho}_{a}V_{0}R_{0}/{\mu}_{a}$, where ${\rho}_{a}$ is the air density, $V_{0}$ the impact velocity, $R_{0}$ is the drop radius, and ${\mu}_{a}$ is the air viscosity. For, the current experimental condition using ${\rho}_{a}=1.293$kg/m3, $V_{0}=0.25$m/s, $R_{0}=1.1$mm, ${\mu}_{a}=1.846{\times}10^{-5}$Pa.s, we have $Re=19.26{\sim}20$. In general, the drag on a spherical drop scales as $F_{v}{\sim}{\rho}_{a}V_{0}^{2}{\pi}R_{0}^{2}C_{d}$, where $C_{d}$ is the drag coefficient. For flow past spheres, $C_{d}{\sim}Re^{-1}$ for $Re<300$. Therefore the drag force on the drop scales as $F_{v}{\sim}{\pi}{\mu}_{a}R_{0}V_{0}$ (i.e. a linear drag law). Therefore, we use a linear drag model, the Stokes drag law. Further the deviation of the actual drag force from the Stokes drag is given by $C_{d}Re/24=F_{v}/6{\pi}{\mu}_{a}R_{0}V_{0}$. For $Re<200$, the factor $(C_{d}Re/24){\sim}\mathcal{O}(1)$ Langmuir & Blodgett (1946). Therefore the drag force is given by $F_{v}{\sim}6{\pi}{\mu}_{a}R_{0}V_{0}$, i.e. the Stokes drag law. Using the Stokes drag as the scale of the viscous force, equation (3.42) scales as

[TABLE]

On simplifying and solving for the scale of $h_{0}$ we have

[TABLE]

Writing equation (3.44) in a non-dimensional manner we have

[TABLE]

where $Ca={\mu}_{a}V_{0}/{\sigma}_{aw}$, and $Bo=({\rho}_{w}-{\rho}_{a})R_{0}^{2}g/{\sigma}_{aw}$ is the Bond number based on the drop radius. The theoretical value of ${\epsilon}$ calculated using equation (3.45) is $1.25{\times}10^{-3}$ and agrees with the experimental value of $1.15{\times}10^{-3}$ within the experimental uncertainty.

3.4 Estimating a scale for the dimple radius during dimple formation process

The bottom surface of the drop and hence the air layer profile prior to dimple formation is given by (refer to equation (3.1))

[TABLE]

The condition for the air profile height becoming equal to $h_{crit}$ at any particular radial coordinate $r_{dimp}$ is given by

[TABLE]

The following inequalities are being satisfied in general impact scenarios. The condition for dimple to form in terms of air layer thickness is given by

[TABLE]

Condition given by equation (3.48) holds due to the following inequality that states, in a general impact scenario where the kinematic time scale of the interface is smaller than the capillary deformation time scale for the dimple formation, i.e.

[TABLE]

Rearranging equation (3.47) we have

[TABLE]

Simplifying equation (3.50) we have

[TABLE]

The dimple radius therefore becomes

[TABLE]

where

[TABLE]

Note equation (3.52) predicts the approximate radius of dimple during its formation stage. Therefore equation (3.52) should provide a dimple scale that should be greater than the radial length scale $R_{0}{\epsilon}^{1/2}$ and smaller than the maximum dimple radius, i.e.

[TABLE]

Using the drop radius $R_{0}$ and the value of ${\epsilon}$ and ${\epsilon}_{crit}$ from equations (3.45) and (3.27) respectively, the inequality provided in (3.54) is satisfied (the numerical translation of inequality (3.54) is $46{\mu}m<142{\mu}m<509{\mu}m$).

3.5 The expansion of the peripheral air disc

Based on the Stokes approximation in the air layer region, the radial pressure gradient balances the viscous stresses and hence the radial momentum equation becomes

[TABLE]

The pressure in the air layer is a superposition of various effects, however capillary is the dominant scale once the dimple has formed. Therefore using the scales ($r{\sim}r_{d}$, $z{\sim}h_{0}$, $p{\sim}2{\sigma}_{aw}/R_{0}$, $u_{r}{\sim}<V>$ where $<V>=dr_{d}/dt$ is the average expansion velocity of the air disc) in equation (3.55)

[TABLE]

[TABLE]

Integrating with respect to $t$ we have

[TABLE]

The radial disc expansion as a function of time becomes

[TABLE]

Fig. 9(d) depicts the air disc expansion both theoretically and experimentally. The experimental and theoretical scales agree within the experimental uncertainities. Notice that from equation (3.59) we have $r_{d}{\sim}t^{1/2}$ (Fig. 9(d)) and the coefficient of $t^{1/2}$ predicted using equation (3.59) agrees with the experiments.

3.6 Spatial length scales in the peripheral air disc

The central air dimple is surrounded by a thin region that we refer to as the peripheral air disc. The peripheral air disc thickness is one order smaller than the central dimple thickness and approaches the mean free path of air. The mean free path (${\lambda}$) of the air molecules is given by

[TABLE]

where ${\mu}_{a}$ is the air viscosity, and $p$, $T$, $R_{s}$ denote air pressure, temperature and specific gas constant respectively. At standard room temperature of $300$K and pressure of $10^{5}{\>}$Pa the mean free path is approximately $67$nm. The ratio of the mean free path of the air to local air layer thickness profile represents a 2D Knudsen field that we use to characterize the air layer beneath the impacting drop. The Knudsen field is defined as

[TABLE]

where $h$ represents the two dimensional air layer thickness profile. Note that the reciprocal of the Knudsen field represents a non-dimensional height profile measure. Figs. 10, 11, 12, 14 depicts the Knudsen field plotted for substrate temperatures of $300$K, $353$K, $423$K, $473$K respectively at various non dimensional time instants. Fig. 5 and Fig. 6 of the supplementary document shows the Knudsen field variation for a larger temporal interval for $T_{s}=423$K and $T_{s}=473$K respectively. Comparing the Knudsen field, we can observe that the structures in the peripheral disc regions are distinct. A unique circular pattern is observed for all the temperature cases in the peripheral disc region where $Kn{\sim}{\mathcal{O}(10^{-1})}$ (refer to Figs. 10(c), 11(e), 12(a), 14(d)). Subjected to very high values of Knudsen number in the peripheral air disc region, the dynamics is governed by an interaction of continuum and non-continuum effects. The flow field in the peripheral disc region lies in the slip flow and transition flow regime in contrast to the central air dimple. Notice the the evolving nature of the structures in the peripheral air disc denoting a time varying unsteady field (Fig. 13). Fig. 13 shows the 3D visualization of a typical air layer thickness profile as a function of non dimensional time for substrate temperature $T_{s}=423$K. The evolving perturbative structures (air layer height fluctuations) in the peripheral air disc region are clearly visible. This is in contrast to the central dimple which is essentially quasi-steady and relatively smooth. The structure in the peripheral air disc and the relatively smooth profile could be understood in terms of capillary flows as discussed below. As the air beneath the impacting drop drains out in the radial direction, the air-water interface at the bottom of the drop is subjected to capillary waves. As the air drains out below the impacting drop, the rise in capillary pressure causes the formation of the dimple. The draining air along with high curvatures caused during during dimple formation produces capillary waves at the bottom air-water interface of the drop (Mandre et al., 2009). Capillary waves could be sustained if the capillary wavelength given by ${\lambda}_{c}$ is smaller than equal to the drop length scale $R_{0}$, i.e.,

[TABLE]

where ${\rho}_{l}$ is the density of the liquid, ${\sigma}_{aw}$ is the air-water surface tension, and $V$ is the air velocity between the drop and the substrate. Equation (3.62) is a necessary but not a sufficient condition for capillary waves to form. For the impact condition in this work, the inequality is satisfied and hence the interface is subjected to capillary perturbations. From Figs. 10-14, we observe that the Knudsen number before the first point of contact in certain locations of the peripheral air disc is one order larger than the mean free path ($Kn{\sim}0.1$). Owing to the very small gap between the substrate and drop, the pressure in the gap region is subject to disjoining pressure that originates due to molecular interactions between surfaces. Derjaguin (Derjaguin & Churaev, 1978) showed that the pressure field in such small gaps becomes anisotropic in comparison to the large scale isotropic pressure. The balance between disjoining and capillary pressure, establishes the existence of new length scales for short time scales.

The molecular (disjoining pressure term) and capillary interaction balance becomes (Roy et al., 2022)

[TABLE]

where $A$ is the Hamacker’s constant ($A{\sim}\mathcal{O}$($4.76{\times}10^{-20}$J)), $h_{pd}$ is the thickness in the peripheral air disc region. Solving for $h_{pd}$ from equation (3.63) becomes

[TABLE]

The wavelength corresponding to thin film instability is given by (Vrij, 1966; Roy et al., 2022)

[TABLE]

where $h$ represents the air layer thickness in the peripheral air disc region. Substituting $h_{pd}$ from equation (3.64) in equation (3.65) we have

[TABLE]

Simplifying equation (3.66) we get

[TABLE]

Fig. 7 in the supplementary document shows the length scale measurement process in the peripheral air disc. Fig. 7(a) in the supplementary document shows a typical Knudsen field for $T_{s}=300$K. The perturbative structures caused due to the combined effect of capillary effects, thin film instability and intermolecular forces can be observed by a distinct cyan color in the disc region. A typical black dotted line shown depicts the measurement pathway of the length scale. Fig. 7(b) in the supplementary document depicts the intensity profile along the black dotted line. The length scale is determined from the width of the peak near the maximum intensity value. The experimental length scale is approximately $12{\mu}$m (${\sim}\mathcal{O}$($10{\mu}$m)) with a standard deviation of $4{\mu}$m. The theoretical capillary wavelength ${\lambda}_{c}$ and the thin film instability length scale ${\lambda}_{TF}$ predicted from equation (3.63) and (3.67) respectively corroborates with the experimental measured scales in the peripheral disc region. The theoretical scale comes out to be of the order of ${\lambda}_{c}{\sim}8{\mu}$m and ${{\lambda}_{TF}}{\sim}7{\mu}$m. The experimental length scale of the order of $10{\mu}$m corroborates with the theoretical values of the capillary wavelength and the thin film wavelength within the experimental uncertainty. The central air dimple interface is smoother than the peripheral air disc interface as could be seen in the various Knudsen field contour plots. The relative smoothness of the interface could be understood based on the capillary wave characteristics given by equations (3.62). We can observe that the capillary wavelength scale ${\lambda}_{c}$ is inversely proportional to the air velocity field between the drop and the substrate. It could be observed from equation (3.36), that the average radial air velocity increases with radial coordinate. Therefore, for higher velocities observed in the peripheral air disc corresponds to lower wavelength capillary waves in comparison to smaller velocities and larger wavelength capillary waves near the central dimple region. The central air dimple is therefore subjected long wave capillary perturbations in contrast to short wave capillary perturbations in the peripheral air disc region. The short wave capillary perturbations along with the van der Waals interactions results in the formation of hole nucleation points that result in air layer dewetting.

For the impact Weber number of unity considered in this work ($We{\sim}1$), the first point of contact occurs through kink mode of contact as shown numerically by Chubynsky et al. For larger Weber number regime that we have observed in our previous work (for impact on liquid pools) we observe the first point of contact occurs through film mode of contact. Therefore, the structures we discover in the peripheral air disc for a variety of imapct events are essentially the experimental realizations of the kink and film mode of contact discussed numerically in some recent articles (de Ruiter et al., 2015b; Chubynsky et al., 2020).

4 Conclusion

In conclusion, we study the air layer dynamics beneath an impacting drop on a heated substrate. We show that a Gaussian profile could represent the evolving central air dimple profile. The geometry of the dimple depends on the substrate temperature very weakly and is a function of the impact Weber number majorly. Further, we observe that the air layer rupture time scale and rupture radius increase with substrate temperature slowly. We characterize the entire air layer profile beneath the drop by a time-varying Knudsen field. We show that a unified treatment of continuum and non-continuum mechanics is required to comprehend the flow field and the local anisotropic structures observed in the peripheral air disc. The airflow dynamics in the central dimple region lies within the continuum Stokes regime, whereas the peripheral air disc falls within the non-continuum slip flow and transition regime owing to the large Knudsen number. However, the average air disc expansion length scale on an average could be explained based on the continuum Stokes approximation. We further show that the structures observed in the peripheral air disc provide the mechanism for the film and kink contact modes between the drop and the substrate. The structures in the peripheral disc region initiate asymmetric contact between the drop and the substrate. The structures in the peripheral air disc occur due to the asymptotic effects of van der Waals and capillary interactions. We hope that the current work will lead to future investigations at the intersection of classical statistical and continuum mechanics to probe the mechanisms and the associated spatiotemporal length scales in the peripheral air disc region in further detail.

Supplementary movie captions

Movie 1: High speed interferometry visualization for drop impact on substrate temperature $T_{s}=300K$

Movie 2: High speed interferometry visualization for drop impact on substrate temperature $T_{s}=353K$

Movie 3: High speed interferometry visualization for drop impact on substrate temperature $T_{s}=423K$

Movie 4: High speed interferometry visualization for drop impact on substrate temperature $T_{s}=473K$

Declaration of Interests

The authors declare no conflict of interest.