2.1. Experimental set-up

An experimental set-up was developed to investigate the heat transfer mechanisms during droplet impact. The configuration is schematically illustrated in figure 1. Droplets of double-distilled water are generated within a temperature-controlled cooling chamber (marked 4 in figure 1), ensuring consistent and well-defined initial and ambient conditions. The temperature of the droplets is measured using a type-T thermocouple, marked 5, with an accuracy of $\pm 0.5\,^\circ$ C, positioned inside the droplet-generation needle. The cooling system, marked 3, ensures precise control of the droplet’s initial temperature. The moving substrate consists of a graphite-coated sapphire disc, marked 8, that serves as the target surface. This substrate is rotated using a geared motor, marked 9, with precise speed control, allowing for accurate setting of the angular velocity. A magnetic encoder, marked 10, is used to provide feedback on the disc’s rotation, which is managed by a motor driver controller, marked 11. To synchronise the measurements, a data acquisition system, marked 16, is used, while environmental parameters, such as humidity and temperature, are monitored using a dedicated sensor, marked 13, with data retrieved by a microcontroller, marked 14. The Sensirion SHT85 sensor measures temperature with $\pm 0.1\,^\circ$ C accuracy and relative humidity with $\pm$ 1.5 % error. The personal computer, marked 15, is used for data logging and processing. The impact radius, $ R$ , is fixed at $ 73 \pm 1$ mm, such that no edge effects are present in the boundary layer (Pier Reference Pier2013). A high-speed camera, marked 2, positioned laterally to the substrate, records the impact dynamics, enabling measurement of the droplet’s normal velocity relative to the substrate, $ U_{{d,n}}$ , and its diameter, $ D$ . The Phantom V2012 high-speed camera operates at 98 900 Hz, with spatial resolution of either 30.9 or 26.7 $\unicode{x03BC} \textrm{m}\,\textrm{pixel}^{-1}$ , depending on the experiment. A high-speed infrared camera, marked 1, is positioned orthogonally to the substrate’s axis of rotation to capture the thermal interactions during droplet impact. The Telops FAST M3k infrared camera covers a spectral range of 1.5–5.4 $\unicode{x03BC}$ m, operates at 2300–8600 Hz, and has a spatial resolution of 73 or 200 $\unicode{x03BC} \textrm{m}\,\textrm{pixel}^{-1}$ , depending on the experiment.

Figure 1. Illustration of the experimental set-up developed to investigate heat transfer during droplet impact. The set-up comprises a droplet-generation system, a rotating sapphire disc as the substrate, high-speed cameras for visualising droplet dynamics, an infrared camera for capturing temperature fields and environmental monitoring equipment.

Figure 2. Representative images of droplet impact on a moving substrate. Panels (a) and (b) show the droplet impact at a substrate velocity of 0.38 m s⁻¹, with a droplet diameter $ D = 2.45 \, \textrm{mm}$ , normal impact velocity $ U_{{d,n}} = 3.13 \, {\textrm{m s}^{- 1}}$ , initial droplet temperature $ T_{{d,0}} = 0.5\,^\circ \textrm{C}$ and initial substrate temperature $ T_{{w,0}} = 26.1\,^\circ \textrm{C}$ . Panels (c) and (d) depict the droplet impact at a substrate velocity of 6.12 m s⁻¹ under otherwise identical conditions. An asymmetry in the droplet spreading is visible, arising from the increased substrate motion.

Exemplary side-view images of the spreading drop are shown in figure 2. These were captured using the high-speed video system, marked 2. Thermocouple measurements are used to extract the corresponding temperature-dependent properties of water from VDI (Reference VDI2013). In the laboratory reference frame, the droplet possesses no tangential velocity. Consequently, its tangential velocity relative to the substrate equals the substrate speed, $ \omega R = U_{{d,t}}$ , where $ \omega$ is the angular velocity of the substrate, derived from the motor’s rotational speed, $ n$ .

The Reynolds number for the airflow induced by the moving substrate is denoted by ${\textit{Re}}_{{s}} = \omega R^2/\nu _{{a}}$ , where $ \nu _{{a}}$ is the kinematic viscosity of air. In these experiments, $ {Re}_{{s}}$ ranges from 2000 to 30 000. According to Kobayashi (Reference Kobayashi1994), the critical Reynolds number for the onset of instabilities in the boundary layer, $ {\textit{Re}}_{{cr}}$ , is $ 4.5 \times 10^4$ . As the maximum $ {\textit{Re}}_{{s}}$ investigated remains below this threshold, the boundary layer surrounding the droplet remains laminar throughout the experimental conditions. Understanding the airflow is crucial, as a strongly turbulent boundary layer, as highlighted by Stumpf et al. (Reference Stumpf, Qezeljeh, Kamal, Dezitter, Martuffo, Roisman and Hussong2024), can lead to droplet deformation and subsequent aerodynamic rebound. The disc is composed of an infrared-transmissive material (aluminium oxide, $\alpha$ - $\text{Al}_2\text{O}_3$ , also known as sapphire). To ensure accurate infrared measurements, the impact surface of the disc is coated with a highly emissive graphite-based paint. The influence of this coating on the measurements is considered negligible, an assumption widely supported in the literature (Chaze et al. Reference Chaze, Caballina, Castanet, Pierson, Lemoine and Maillet2019; Gholijani et al. Reference Gholijani, Schlawitschek, Gambaryan-Roisman and Stephan2020; Schmidt et al. Reference Schmidt, Breitenbach, Roisman and Tropea2020). The estimated time required for the thermal boundary layer within the coating to develop fully is given by $t \propto \epsilon ^2/\alpha _{{c}}$ , where $\epsilon$ denotes the thickness of the coating, and $\alpha _{{c}}$ is its thermal diffusivity. In these experiments, the average thickness of the coating is $ 8 \, \unicode{x03BC} \textrm{m}$ , and the thermal diffusivity is approximately $ 2 \times 10^{-6} \, \textrm{m}^2\,\textrm{s}^{- 1}$ . This yields an estimated development time of $t \approx 0.032$ ms, which is significantly shorter than the temporal resolution of $\varDelta t = 0.15$ ms given by the infrared camera. The infrared camera is positioned beneath the disc to capture transient temperature measurements during droplet impact. In the experiments, droplet temperature, diameter and substrate speed are systematically varied. The substrate temperature, monitored with temperature sensors, varied with ambient room conditions and ranged from $ 24.1$ – $ 33.0 \,^\circ \textrm{C}$ across all experiments.

Table 1 summarises the experimental conditions tested in this study. A full factorial design was used to investigate the interactions between these parameters, with the droplet’s normal impact velocity held constant.

Table 1. Summary of the experimental conditions tested in this study. A full factorial design was used to systematically vary droplet diameter ( $D$ ), substrate speed ( $U_{{d,t}}$ ) and droplet temperature ( $T_{{d,0}}$ ). The droplet’s normal impact velocity ( $U_{{d,n}}$ ) was held constant, while key non-dimensional parameters, including the normal and tangential Reynolds ( ${\textit{Re}}_{{d,n}}$ , ${\textit{Re}}_{{d,t}}$ ) and Weber numbers ( ${\textit{We}}_{{d,n}}$ , ${\textit{We}}_{{d,t}}$ ), were computed over the range of experimental conditions.

2.2. Calibration of the infrared camera

The infrared camera is calibrated in situ using the set-up illustrated in figure 3. A copper plate, cooled by a chiller, extends beyond the camera’s entire field of view (FOV), ensuring uniform calibration across the measurement area.

Figure 3. Schematic of the calibration set-up for the infrared camera. The set-up includes a temperature-controlled copper plate thermally coupled to the sapphire disc via a thin layer of thermal paste, ensuring uniform heat transfer. A thermocouple is embedded at the surface to monitor and maintain the substrate temperature accurately. The sapphire disc, coated with a high-emissivity material, serves as the substrate for infrared measurements. The infrared camera captures the temperature distribution at the coating’s bottom surface within its FOV.

Temperature images are not directly obtained from the infrared camera; instead, the digital level ( $DL$ ) recorded by the camera is converted into temperature values. The calibration process involves incrementally increasing the temperature of a reference surface from $16$ – $35\,^\circ \textrm{C}$ , with calibration images captured at $1\,^\circ \textrm{C}$ intervals. The variation in $DL$ is recorded as a function of the surface temperature ( $T$ ), establishing a calibration curve. This relationship between $DL$ and $T$ is governed by Planck’s law, which describes the spectral radiant emittance ( $W_{b}$ ) of an object as a function of its temperature,

(2.1) \begin{align} W_{{b}}(\lambda , T_{\textit{obj}}) = \frac {2 \pi h c^2}{\lambda ^5} \frac {1}{\left (e^{\frac {h c}{\lambda k T_{\textit{obj}}}} - 1\right )}, \end{align}

where $h$ is Planck’s constant, $k$ is Boltzmann’s constant, $c$ is the speed of light, $\lambda$ is the wavelength and $T_{\textit{obj}}$ is the temperature of the object.

Numerous studies, including those by Martiny et al. (Reference Martiny, Schiele, Gritsch, Schulz and Wittig1996), Schulz (Reference Schulz2000), Ochs et al. (Reference Ochs, Horbach, Schulz, Koch and Bauer2009), Elfner et al. (Reference Elfner, Schulz, Bauer and Lehmann2017), Falsetti, Sisti & Beard (Reference Falsetti, Sisti and Beard2021) and Sisti et al. (Reference Sisti, Falsetti, Beard and Chana2021), have based their infrared camera calibration on the above equation. By rearranging the terms and isolating $T$ , the following semiempirical relation is obtained:

(2.2) \begin{align} T = \frac {r}{\ln \left (\frac {b}{I_{\textit{obj}}} + f\right )}, \end{align}

where $I_{\textit{obj}}$ is the digital level recorded by the camera, and $r$ , $b$ and $f$ are calibration parameters. These parameters were estimated in situ using a thermocouple as the reference and a nonlinear least-squares fitting procedure applied to the calibration images obtained with the set-up illustrated in figure 3.

From the calibration process, two types of errors are quantified: the temporal noise in the temperature readings and the corresponding noise-equivalent heat flux. The temporal noise, $\sigma _{{T}}$ , is defined as

(2.3) \begin{equation} \sigma _{{T}} = \frac {1}{N_{\textit{pixels}}} \sum _{j=1}^{N_{\textit{pixels}}} \sigma _{{T,j}}, \end{equation}

where $\sigma _{{T,j}}$ represents the standard deviation of temperature fluctuations for pixel $j$ across multiple frames $i$ , given by

(2.4) \begin{equation} \sigma _{{T,j}} = \sqrt {\frac {1}{N_{\textit{frames}}} \sum _{i=1}^{N_{\textit{frames}}} \left (T_{{i,j}} - \bar {T}_{{j}}\right )^2}. \end{equation}

Here, $N_{\textit{pixels}}$ is the total number of pixels in the region of interest, and $N_{\textit{frames}}$ is the number of frames captured during the calibration measurement. Here $T_{{i,j}}$ denotes the temperature of pixel $j$ in frame $i$ , and $\bar {T}_{{j}}$ is the average temperature of pixel $j$ across all frames. The overall temperature noise, $\sigma _{{T}}$ , is computed as the average of $\sigma _{{T,j}}$ over all pixels within the region of interest, assuming a spatially uniform temperature distribution.

The calculated temporal noise, $\sigma _{{T}}$ , is $0.21\,^\circ \textrm{C}$ . This value quantifies the variation in temperature readings over time. The corresponding noise-equivalent heat flux is determined separately through numerical simulation of the transient Fourier heat equation applied to the substrate. By analysing the artificial heat flux generated by temperature fluctuations between consecutive frames, the heat flux noise is estimated to be $10^4 \, \textrm{W m}^{- 2}$ . Further details on the heat flux calculations are provided in § 2.4.

2.3. Temperature measurements

The experimental contact temperature values were directly obtained from the infrared camera measurements. A time-lapse sequence of the thermal evolution during droplet impact is shown in figure 4. The colder region represents the area wetted by the droplet, delineated by a sharp temperature gradient in the vicinity of the three-phase contact line. The contact area between the droplet and the substrate increases over time. Figure 4(a–d) correspond to four successive time steps: 0.57 ms, 1.03 ms, 1.61 ms and 2.07 ms after impact. As time progresses, the expansion of the contact area results from both droplet spreading and entrainment due to the relative motion with the substrate. To determine a single representative value for the contact temperature, $T_{{c}}$ , the average temperature over the droplet–surface contact area was used. This approach is justified by the uniformity of the temperature across the contact region, as evident in figure 4. A slight deviation, specifically a lower temperature difference of approximately $0.2 {-} 0.3 \,^\circ \textrm{C}$ , is observed in the initial contact area of the droplet, as evidenced in figure 4(d). This phenomenon can be attributed to the influence of the droplet’s initial normal velocity, which reduces the contact resistance (Aziz & Chandra Reference Aziz and Chandra2000) and consequently results in a slightly lower contact temperature in this region.

The contact temperature evolutions derived from experiments are given as a function of time for four different substrate velocities in figure 5(a–d). In addition to the experimental values and their associated errors, we compare our model with the models of Seki et al. (Reference Seki, Kawamura and Sanokawa1978) and Roisman (Reference Roisman2010), (1.2) and (1.3), respectively. These models predict a constant contact temperature, marked in figure 5 by dashed lines. At the early stages of impact, experimental results show a different behaviour, with the contact temperature decreasing monotonically before reaching the predicted values. We hypothesise that the thermal resistance at the substrate surface causes the temperature evolution. In this study, an improved model has been proposed, which accounts for both convective effects and thermal resistance at the interface. Exemplary results of this model are plotted as blue curves in figure 5(a–d). The theoretical framework for the current model is detailed in § 3.

Figure 4. Time-lapse sequence showing the thermal evolution during the impact of a droplet on a moving substrate, recorded using an infrared camera. The droplet, with a diameter of 1.95 mm and an initial temperature of 0.5 $\,^\circ$ C, impacts a substrate rotating at 800 r.p.m. (6.12 m s⁻¹). Snapshots are taken at different times following impact: 0.57 ms, 1.03 ms, 1.61 ms and 2.07 ms. The temperature distribution within the droplet-substrate contact region remains mostly uniform, with slight deviations observed at the leading edge during initial contact.

Figure 5. Contact temperature evolution ( $T_{{c}}$ ) over time at different rotational speeds ( $n$ ). The panels correspond to (a) $n = 200\,\text{r.p.m.}$ , (b) $n = 400 \, \text{r.p.m.}$ , (c) $n = 600 \, \text{r.p.m.}$ and (d) $n = 800 \, \text{r.p.m.}$ , for a droplet diameter of $D = 2.45 \, \text{mm}$ and an initial droplet temperature of $T_{d,0} = 0.5 \, \,^\circ \text{C}$ . The experimental data (red curve) is compared with the theoretical predictions of the current model (blue curve), with the shaded area indicating the measurement error of the experimental data.

2.4. Heat flux evaluation

Heat flux at the wall cannot be directly retrieved from the infrared camera; therefore, a numerical simulation is performed based on the contact temperature measured by the camera. The Fourier heat diffusion equation is solved within the substrate as follows:

(2.5) \begin{equation} \frac {\partial T}{\partial t} + \boldsymbol{w} \boldsymbol{\cdot }\boldsymbol{\nabla }T = \alpha _{{w}} {\nabla} ^2 T, \end{equation}

where $\alpha _{{w}}$ is the thermal diffusivity of the substrate, and the velocity vector $\boldsymbol{w}$ represents the substrate’s motion, as the measurements are taken in the laboratory reference frame. Thermal properties are assumed to remain constant throughout the analysis.

The computational mesh is designed such that the top layer emulates the camera pixels in a one-to-one mapping, with each pixel assigned a velocity vector derived from calibration images acquired during the motion-calibration procedure. The mesh spacing through the substrate thickness and the numerical time step are $3.8\,\unicode{x03BC} \text{m}$ and $30.16\,\unicode{x03BC} \text{s}$ , respectively. Grid- and time-independence studies are reported in Appendix A. At the interface, a Dirichlet boundary condition is applied, where the temperature is specified using data from the infrared camera. Along the sides of the domain, zero-gradient boundary conditions are imposed. At the bottom, a zero-gradient condition is also applied. This choice is based on the consideration that the substrate is assumed to remain at room temperature, and the heat flux at the bottom is expected to be negligible.

The validity of the assumption is assessed by analysing the effect of the entrained air within the boundary layer, which can lead to heating and a slight increase in temperature. To estimate this effect, the problem is treated analytically by solving a system of ordinary differential equations (ODEs) as described by Millsaps & Pohlhausen (Reference Millsaps and Pohlhausen1952), with the boundary condition modified to impose an adiabatic wall instead of a constant temperature. The details of the scaled temperature solution are provided in Appendix B.

Under these conditions, the maximum expected temperature increase for the highest substrate speed is

(2.6) \begin{equation} \Delta T_{\textit{max}} = \frac {(R\omega )^2}{c_{{p,a}}S(0)} + \frac {\nu _{{a}} \omega }{c_{{p,a}}}Q(0) \approx 0.2\, ^\circ \textrm{C}, \end{equation}

where $S(0)$ and $Q(0)$ are reduced parameters evaluated at the wall based on the ODE solution.

Aside from the theoretical description, this statement is further supported by the experimental findings of Cardone, Astarita & Carlomagno (Reference Cardone, Astarita and Carlomagno1994). While the Nusselt number $\textit{Nu}$ reported in Cardone et al. (Reference Cardone, Astarita and Carlomagno1994) could, in principle, be used to calculate the heat transfer coefficient $h_{\textit{tc}}$ for implementing a Robin-type boundary condition, the combination of a very small temperature difference and the uncertainty associated with the coefficient led us to adopt a zero-gradient boundary condition instead. Moreover, it is well established that for elliptic and parabolic equations, boundary effects diminish rapidly with increasing distance from the boundary (Evans Reference Evans2010).

Finally, preliminary simulations with both bottom boundary conditions (adiabatic and Robin) show that the heat flux at the impact-side interface ( $z=0$ ) is unchanged. Any differences occur only near the bottom of the substrate, where the magnitude of $\varDelta q = q_{{Robin}}-q_{{Neumann}}$ is one to two orders of magnitude below the experimental noise, and thus irrelevant for the reconstructed interfacial flux.

3.1. Problem formulation

The flow and temperature fields in the spreading droplet must satisfy the continuity, momentum and energy equations (Bird, Stewart & Lightfoot Reference Bird, Stewart and Lightfoot1966). We assume incompressible flow, as the time scale and spatial region affected by compressibility are negligible in our experiments. This assumption has been thoroughly discussed and validated by Weiss & Yarin (Reference Weiss and Yarin1999). Under these conditions, the governing mass balance, momentum balance and energy balance equations for the Newtonian flow $\boldsymbol{\upsilon }$ , pressure $p$ and temperature $T$ are written as

(3.1a) \begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{\upsilon } &= 0, \end{align}

(3.1b) \begin{align} \rho \frac {\partial \boldsymbol{\upsilon }}{\partial t} + \rho (\boldsymbol{\upsilon } \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol{\upsilon } &= -\boldsymbol{\nabla } p + \mu {\nabla} ^2 \boldsymbol{\upsilon }, \end{align}

(3.1c) \begin{align} \rho c_p \left ( \frac {\partial T}{\partial t} + \boldsymbol{\upsilon } \boldsymbol{\cdot }\boldsymbol{\nabla } T \right ) &= k {\nabla} ^2 T. \end{align}

The decoupling of the hydrodynamic problem from the thermal problem is a direct result of assuming constant thermophysical properties. This simplifies the analysis, allowing the temperature field to be solved independently of the velocity and pressure fields. This assumption is appropriate for the regime considered: temperature changes in the liquid remain moderate and are confined to a thin layer near the interface, variations of $\rho$ , $c_p$ and $k$ are much smaller than the change in viscosity, which, in turn influences the thermal problem only indirectly via the advective correction to the similarity velocity $w(z,t)$ . A more detailed validation and sensitivity analysis is provided in Appendix C. However, while valid in our parameter range, at larger temperature gradients (leading to more pronounced property variations), at higher shear rates where viscous heating is non-negligible, or when thermocapillary effects are sufficiently strong to modify the velocity field, a more thorough assessment would be required.

Figure 6. Illustration of droplet impact on a rough substrate, highlighting the interface between the droplet and the substrate. The zoomed-in view emphasises the roughness of the substrate and the liquid–substrate contact area, with the white regions representing entrapped air. The red arrows in the zoomed area represent the heat flux, illustrating local variations where regions with better contact exhibit higher heat flux, while areas with entrapped air show reduced heat flux.

Consider now the Cartesian coordinate system $\{x,y,z\}$ with the unit base vectors $\{\boldsymbol{e}_x, \boldsymbol{e}_y, \boldsymbol{e}_z\}$ with $\boldsymbol{e}_z$ normal to the substrate interface directed towards the drop region. The velocity field in the fluid, which also accounts for the viscous effects, is defined as $\boldsymbol{\upsilon } = u\boldsymbol{e}_x + \upsilon \boldsymbol{e}_y + w\boldsymbol{e}_z$ . The initial conditions and the boundary conditions far from the interface $z=0$ are expressed as

(3.2a) \begin{align} (\boldsymbol{\upsilon } - \boldsymbol{\upsilon }_0) \times \boldsymbol{e}_z = 0, & \quad T = T_{{d,0}} \quad\, \text{at } z \to \infty \; \forall \, t \gt 0 \;\,\,\,\, \text{ and at } t = 0 \; \forall \, z \gt 0, \end{align}

(3.2b) \begin{align} \boldsymbol{\upsilon } = 0,&\quad T = T_{{w,0}} \quad \text{at } z \to -\infty \; \forall \, t \gt 0 \; \text{ and at } t = 0 \; \forall \, z \lt 0, \end{align}

where $\boldsymbol{\upsilon }_0$ is the velocity in the outer region of the drop far from the viscous boundary layer. The expression for the velocity $\boldsymbol{\upsilon }_0$ is obtained as the solution for the inviscid flow in the drop (Yarin & Weiss Reference Yarin and Weiss1995; Roisman Reference Roisman2010). Here, $T_{{d,0}}$ and $T_{{w,0}}$ denote the initial temperatures of the droplet and the substrate, respectively.

The no-slip and zero-flux boundary conditions yield the following boundary conditions for the flow at the interface:

(3.3) \begin{equation} \boldsymbol{\upsilon }=0 \quad \textrm{at}\quad z =0. \end{equation}

In this study, the temperature discontinuity arising from the presence of contact resistance between the substrate and the spreading droplet is taken into account. The physical meaning of this thermal resistance is illustrated in figure 6. The thermal effects are influenced by the morphology of the substrate and air entrapment at the interface. The region in the substrate and in the drop, associated with the thermal resistance, is characterised by an effective thickness $h_{{R}}$ and effective heat conductivity $k_{{R}}$ . In this case the heat flux can be estimated in the form $q \approx k_{{R}} (T_{{d}}-T_{{w}})/h_{{R}}$ , where $h_{{R}}/k_{{R}} \equiv R_{{c}}$ is the thermal contact resistance as already proposed by Kapitza (Reference Kapitza1941), Swartz & Pohl (Reference Swartz and Pohl1989).

The thermal boundary conditions at the interface are expressed as

(3.4) \begin{equation} q = k_{{w}}\frac {\partial T}{\partial z} \bigg |_{\textit{wall}} = k_{{d}}\frac {\partial T}{\partial z} \bigg |_{\textit{drop}}, \quad T_{{d}} - T_{{w}} = qR_{{c}},\quad \textrm{at}\quad z =0. \end{equation}

The first equation in (3.4) represents the energy balance at the surface, neglecting the heat accumulated due to the heating of the microscopic contact region, and the second term represents the temperature jump at this contact region.

In this study, we treat $R_c$ as a single, casewise constant effective value, i.e. a lumped, area-averaged interfacial impedance over the inner contact zone. As microroughness and lamella pressure create heterogeneous microcontacts, the interfacial impedance is a function of time and space $(R_c = R_c(x,y,t))$ ; however, within our 1-D framework, only the area-averaged resistance enters the boundary condition. This simplification is supported by spatially and temporally smooth infrared temperature fields in the central patch of our wetted area (figure 4). We explicitly restrict model–data comparisons with the footprint interior, away from the lamella and rim where a constant $R_c$ assumption would be violated. The same constant $R_c$ closure is common in related droplet-impact and thermal-spray analyses to capture interfacial cooling processes (Dhiman & Chandra Reference Dhiman and Chandra2005; Nastic, Pershin & Mostaghimi Reference Nastic, Pershin and Mostaghimi2024). Using this boundary condition, rather than enforcing thermal continuity, reflects the physical reality of many practical systems, where contact resistance cannot be neglected, and highlights a key distinction from prior studies of drop-based cooling that assume perfect thermal continuity at the interface.

3.2. Viscous flow in the drop region

The similarity solution for flow within a lamella, expressed in Cartesian coordinates system $\{x,y,z\}$ , fixed at the substrate surface, has been obtained in Roisman (Reference Roisman2009, Reference Roisman2010),

(3.5) \begin{equation} \boldsymbol{\upsilon } = f(\xi ) \frac {(x + U_x t) \boldsymbol{e}_x + y \boldsymbol{e}_y}{t} - 2 \sqrt {\frac {\nu }{t}} g(\xi ) \boldsymbol{e}_z, \quad \xi \equiv \frac {z}{\sqrt {\nu t}}. \end{equation}

Here, $f(\xi )$ and $g(\xi )$ are, respectively, the dimensionless similarity profiles for the tangential and normal velocity components of the boundary-layer flow, where $\xi$ is the similarity variable. At large values of $\xi$ , this field approaches the outer inviscid flow in the spreading drop. The outer solution (Yarin & Weiss Reference Yarin and Weiss1995; Roisman Reference Roisman2010), which satisfies exactly the mass and momentum balance equations is

(3.6) \begin{equation} \boldsymbol{\upsilon }_0 = \frac {(x + U_x\tau )\boldsymbol{e}_x + y\boldsymbol{e}_y - 2z\boldsymbol{e}_z}{t}, \end{equation}

where $U_x$ is the tangential, $x$ -component of the impact velocity relative to the substrate. Here, $\tau$ is a time constant; for the present initial condition it takes the value $\tau \approx 0.25\,D_0/U_0$ (see Roisman et al.  Reference Roisman, Berberović and Tropea2009).

Substituting the expression (3.5) in the mass and momentum balance (3.1a ) and (3.1b ) yields the following system of ODEs:

(3.7a) \begin{align} f =& g', \end{align}

(3.7b) \begin{align} f'' + \frac {f'}{2} (\xi + 4g) + f =& f^2, \end{align}

which can be solved computationally, subject to the no-penetrability conditions and the no-slip conditions at the wall and continuity of the velocity in the outer region

(3.8a) \begin{align} f &= 1, \quad \text{as } \xi \to \infty , \end{align}

(3.8b) \begin{align} g' &= g = 0, \quad \text{as } \xi \to 0. \end{align}

For completeness, the functions $g(\xi )$ and $f(\xi )$ first reported in Roisman (Reference Roisman2010) are plotted in figure 7.

Figure 7. Numerical solutions for the scaled velocity components as functions of the similarity variable $\xi$ . Panel (a) shows the scaled vertical velocity $g(\xi )$ , while panel (b) presents the scaled in-plane velocity $f(\xi )$ , representing the velocity components in the $x$ – $y$ plane. These results are derived from the coupled ODEs described in (3.7a ), (a) The scaled vertical velocity component $g(\xi )$ as a function of $\xi$ . (b)The scaled in-plane velocity component $f(\xi )$ (in the $x$ – $y$ plane) as a function of $\xi$ .

3.3. Heat transfer in the drop and substrate at high Prandtl numbers ${\textit{Pr}}_{d}\gg 1$

Assuming contact resistance, the temperature fields in the spreading drop $z\geqslant 0$ and the solid substrate $z\leqslant 0$ are represented using the known analytical solution for heat conduction. For the solid region, the solution that exactly satisfies the energy equation together with the boundary and initial conditions is given in Carslaw & Jaeger (Reference Carslaw and Jaeger1959) as

(3.9a) \begin{align} T_{{w}} &= T_{{c0}} + (T_{{w0}}-T_{{c0}})\left [\exp \left (\frac {t}{t_R}-2 \zeta _w \sqrt {\frac {t}{t_{{R}}}} \right )\textrm{erfc}\left (\sqrt {\frac {t}{t_{{R}}}}-\zeta _w \right )-\textrm{erf}(\zeta _w) \right ]\!, \end{align}

(3.9b) \begin{align} \zeta _w &\equiv \frac {z}{2\sqrt {\alpha _{{w}} t}}, \end{align}

where $T_{c0}$ is a constant temperature equal to the asymptotic value of the contact temperature at long times, and $t_R$ is a characteristic time associated with the duration of the propagation of the thermal boundary layer inside the contact region of thickness $h_R$ . These constants will be determined from the boundary conditions. Here $\zeta _w$ is the similarity variable in the solid region. The energy equation in the liquid region has to account for the heat convection. Therefore, in contrast to the purely diffusive semi-infinite solution used for the solid (Carslaw & Jaeger Reference Carslaw and Jaeger1959), we adopt a similarity ansatz that extends that diffusion form to include advection. We thus write for the spreading drop

(3.9c) \begin{align} T_{{d}} &= T_{{c0}} + (T_{{d0}}-T_{{c0}})\left [c_{{d}} \exp \left (\frac {t}{t_{{R}}}+2 \zeta _d \sqrt {\frac {t}{t_{{R}}}} \right )\textrm{erfc}\left (\sqrt {\frac {t}{t_{{R}}}}+\zeta _d \right )+\theta (\zeta _d) \right ]\!, \end{align}

(3.9d) \begin{align} \zeta _d &\equiv \frac {z}{2\sqrt {\alpha _{{d}} t}}, \end{align}

where $\alpha _{{d}}$ is the thermal diffusivity of the drop, $c_d$ is a constant and the function $\theta (\zeta _d)$ is the self-similar component of the temperature field associated with convective transport in the drop. Here $\zeta _d$ is the similarity variable in the liquid region. In essence, the exponential factor in (3.9a ) encodes the effects of finite contact resistance, and the classical $\operatorname {erf}(\zeta )$ term describes the response of a semi-infinite conductor under pure thermal diffusion. In the droplet, this term is replaced by the self-similar convective kernel $\theta$ to account for thermal diffusion and advection in the flow. The kernel $\theta$ is obtained from early-time asymptotics of the coupled advection–diffusion problem, as described later.

Expression (3.9a ) satisfies the heat conduction equation in the wall exactly. However, the form (3.9c ) yields only an approximate solution for the temperature field in the drop if the convective terms are smaller than the conductive terms. In the case of a large Prandtl number,

(3.10) \begin{align} {\textit{Pr}}_{{d}} \equiv \sqrt {\frac {\nu _{{d}}}{\alpha _{{d}}}}, \end{align}

the thickness of the thermal boundary layer is much smaller than that of the viscous boundary layer. In this case, since both the value of the normal velocity and its gradient at the wall are equal to zero, the convection effects are indeed small. The validity of this assumption is assessed in Appendix D, where we compare the closed-form model with an independent advection–diffusion numerical reference. The asymptotic expression for the gradient of the drop temperature taking into account a small value of $t_R/t$ is obtained from (3.9c ) in the form

(3.11) \begin{equation} \frac {\partial T_{{d}}}{\partial z} \approx \frac {(T_{{d0}}-T_{{c0}})}{2 \sqrt {\alpha _{{d}} t}}\left [\exp \left(-\zeta _d^2\right) \theta '(0)+\theta '(\zeta ) \right ],\quad \textrm{if}\quad t\gg t_{{R}}. \end{equation}

This approximate expression can be used for the estimation of the small convective terms in the energy equation if the Prandtl number is large.

Figure 8. Theoretically predicted values of $\theta (\zeta _d)$ at various values of the Prandtl number.

Substitution of expressions (3.9) and (3.11) in the energy (3.1c ) yields the following similarity equation:

(3.12) \begin{equation} \theta ''(\zeta _d) + 2 \zeta _d \theta '(\zeta _d) + 4 {\textit{Pr}}_{{d}}^{1/2} g\left [\frac {\zeta _d }{{\textit{Pr}}_{{d}}^{1/2}}\right ] \left \{\theta '(\zeta _d) + \exp \left(-\zeta _d^2\right)\theta '(0)\right \}=0, \end{equation}

where the dimensionless function $g$ is determined from the similarity solution for the drop velocity field.

This ODE has to be solved subject to the boundary conditions

(3.13a) \begin{align} \theta &= 0, \quad \textrm{at} \quad \zeta _d =0 ,\end{align}

(3.13b) \begin{align} \theta &\rightarrow 1, \quad \textrm{as} \quad \zeta _d \rightarrow \infty . \end{align}

Figure 9. Theoretically predicted values of $\theta '(0)$ and $c_{{d}}$ as a function of the Prandtl number.

Figure 10. Comparison of experimental contact temperatures with theoretical predictions for various models: (a) Seki et al. (Reference Seki, Kawamura and Sanokawa1978); (b) Roisman (Reference Roisman2010); (c) the current model. Each plot includes experimental data (circles), a line of perfect agreement (black dashed) and experimental error bounds (red dashed).

Figure 11. Heat transfer coefficient, $h_{\textit{tc}}$ , as a function of time, $t$ , for different rotational speeds and droplet diameters. Experimental data (solid lines) are compared with theoretical predictions (dashed lines). Panels (a) and (b) correspond to substrate rotational speeds of 200 and 800 r.p.m., respectively. Shaded regions indicate measurement uncertainty.

In figure 8, the theoretical predictions for the dimensionless function $\theta (\zeta _d)$ are shown for various Prandtl numbers. The computed values of $\theta '(0)$ as a function of ${\textit{Pr}}_{{d}}$ are shown in figure 9. The value of the constants $T_{c0}$ and $t_R$ are determined from the boundary conditions (3.4) associated with the thermal contact resistance

(3.14) \begin{equation} T_{{c0}} = \frac { e_{{w}} T_{{w0}} + e_{{d}} c_{{d}} T_{{d0}}}{ e_{{w}} + e_{{d}} c_{{d}}}, \quad t_{{R}} = \left (\frac {e_{{d}} e_{{w}}}{e_{{d}}+e_{{w}}} R_{{c}}\right )^2, \quad c_{{d}} = \frac {\sqrt {\pi } \theta '(0)}{2}, \end{equation}

where $\theta (\zeta _d)$ is a dimensionless function given in (3.9c ), $e_{{d}}$ and $e_{{w}}$ are the thermal effusivities of the drop and substrate materials, respectively. Finally, the substrate temperature at the interface $z=0$ and the heat flux are obtained using (3.9) in the form

(3.15a) \begin{align} T_{{w},{\textit{interface}}} &= T_{{c0}} + (T_{{w0}} - T_{{c0}}) \text{exp}\left (\frac {t}{t_{ R}}\right ) \text{erfc}\left (\sqrt {\frac {t}{t_{ R}}}\right )\!, \end{align}

(3.15b) \begin{align} q_{{\textit{interface}}} &= \frac {T_{{d0}} - T_{{w0}}}{\sqrt {t_{R}}} \frac {e_{w}\,e_{ d}\,\sqrt {\pi }\,\theta '(0)} {2\,e_{ w} + e_{ d}\,\sqrt {\pi }\,\theta '(0)} \text{exp}\left (\frac {t}{t_{ R}}\right ) \text{erfc}\left (\sqrt {\frac {t}{t_{ R}}}\right ) .\end{align}

Equations (3.15a ) and (3.15b ) incorporate the influence of substrate morphology, resulting in a time-dependent interfacial temperature. Unlike (1.2) and (1.3), they preserve both the conductive mechanisms of Seki et al. (Reference Seki, Kawamura and Sanokawa1978) and the convective effects described by Roisman (Reference Roisman2010). The value of the thermal contact resistance, $R_{{c}}$ , required to determine the characteristic time scale $t_{{R}}$ was obtained by fitting the analytical substrate–temperature evolution (3.15a ) to each experimental curve and subsequently averaging the results. The mean value was found to be $2.39 \times 10^{-5}$ $\text{m} ^2{\rm K\, W^{- 1}}$ with a 95 % confidence interval of [ $1.73 \times 10^{-5}$ , $3.65 \times 10^{-5}$ ] $\text{m} ^2{\rm K\, W^{- 1}}$ .

A comparison between the experimental results and the three models can be seen in figure 10(a–c). For the models of Seki et al. (Reference Seki, Kawamura and Sanokawa1978) and Roisman (Reference Roisman2010), which predict a time-independent interfacial temperature, we perform the comparison at the instant when the measured interfacial temperature attains its local minimum. The model of Seki et al. (Reference Seki, Kawamura and Sanokawa1978) has a rather fair agreement as shown in figure 10(a), which can be explained by the fact that for our investigated ${\textit{Pr}}_{{d}}$ , the convective effects are small and the underestimation of the contact temperature is balanced by the missing contact resistance condition. On the other hand, the model of Roisman (Reference Roisman2010) consistently underestimates the contact temperature due to the missing contact resistance term. The present model incorporates both heat convection and contact resistance, resulting in a very good agreement with the experimental data.

The heat transfer coefficient at the interface is readily derived by scaling the heat flux to the initial temperature difference between the substrate and the spreading drop as

(3.16) \begin{equation} h_{\textit{tc}} = \frac {q}{\varDelta T} = \frac {e_{{w}} e_{{d}} \sqrt {\pi }\theta '(0)}{ \sqrt {t_{{R}}}(2e_{{w}}+e_{{d}} \sqrt {\pi } \theta '(0))} \exp \left (\frac {t}{t_{{R}}}\right )\textrm{erfc}\left (\sqrt {\frac {t}{t_{{R}}}}\right )\!, \end{equation}

where $\varDelta T = T_{{d0}} - T_{{w0}}$ is the initial temperature difference between the droplet and the substrate. This quantity is plotted in figures 11(a) and 11(b) for the two extreme substrate velocities tested ( $\approx$ 200 and 800 r.p.m.), each at varying droplet diameters. Intermediate velocities were also investigated and exhibited qualitatively similar trends; they are therefore omitted here for conciseness. As expected, neither the droplet diameter nor the substrate velocity significantly affects the heat flux. This is consistent with the theoretical model predictions. The results demonstrate good agreement between experimental data (solid lines) and theoretical predictions (dashed lines), as shown in figure 11, where shaded regions represent measurement uncertainty. The present model applies to early-time, conduction-dominated exchange with finite interfacial resistance and no phase change. Two practical bounds are worth noting: at very high substrate speeds, large shear rates can introduce viscous heating and strong advection, which are not included here; and for very small (micronscale) droplets, the semi-infinite-liquid and locally 1-D assumptions may fail and free-surface effects may become important. Conversely, for larger droplets, the assumptions remain valid. Furthermore, regimes involving phase change (boiling, evaporation, solidification) are outside the scope of this study. For completeness, we propagate the estimated uncertainty in the interfacial resistance $R_c$ into the predicted heat-transfer coefficient $h_{tc}(t)$ ; the procedure and the resulting 95 % confidence intervals are given in Appendix E. These intervals are consistent with the experimental ones. While the contact–temperature points in figure 10 differ only modestly through the models, the corresponding heat flux, and therefore $h_{tc}(t)$ , shows larger differences between the models. This arises from including the contact–resistance boundary condition, which yields a finite heat flux at $t=0$ , in contrast to the initial singularity predicted by the models of Seki et al. (Reference Seki, Kawamura and Sanokawa1978) and Roisman (Reference Roisman2010). This finite behaviour highlights the importance of incorporating realistic boundary conditions in heat transfer models to accurately capture interfacial dynamics during droplet impact.