2.1 Quasi-steady state assumption

Leidenfrost drops levitate above a thin layer of vapour, of good insulating properties since $k_a \ll k$, where $k_a$ and $k$ are the air and water thermal conductivities, provided in Table 1 of the supplementary material. Evaporation, which mainly takes place at the drop base, is thus markedly reduced and we verify here that the drop is at quasi-static equilibrium. A water drop with initial radius $R_0$ initially close to $\approx 4$ mm is deposited on a plate brought to 300 $^{\circ }$C. We use a slightly curved substrate to immobilize the highly mobile liquid. The drop is observed using a top-view high speed camera, from which we extract $R$, the drop equatorial radius as evaporation proceeds. Figure 1 shows that $R$ decreases linearly with time $t$ at a rate $\mathrm {d}R/\mathrm {d}t\approx - 22\ \mathrm {\mu }{\rm m}\ {\rm s}^{-1}$, so that the drop survives approximately $\tau _0\approx 3$ minutes. Temporal variations in $R$ are best fitted by a linear law, plotted as the blue dotted line and with equation

(2.1)\begin{equation} R(t) = R_0(1-t/\tau_0),\end{equation}

where $R_0 = 3.7\pm 0.1$ mm and $\tau _0 = 176.5$ s. This time is much larger than the characteristic time of the internal motion $R/V$, since tracers inside the liquid and at the drop surface reveal flow velocities $V$ as high as a few centimetres per second. The separation of time scales $\tau _0 \gg R/V \approx 0.02$ s suggests that the evaporation-driven dynamics can be decoupled from the inner dynamics. As a result, a Leidenfrost drop can be considered in quasi-static equilibrium at any time, a key assumption to discuss the stability of Leidenfrost inner flows. Moreover, as will be discussed in § 4.3.2, the instability develops faster than $\tau _0$, which supports the quasi-static stability analysis of the Leidenfrost drop.

Figure 1. (a) Radius $R$ of a Leidenfrost drop levitating on a plate heated at 350 $^{\circ }$C as a function of time $R$ decreases as $R(t)=R_0(1-t/\tau _0)$ ((2.1)), plotted as the dotted line, denoting $R_0 = 3.7\pm 0.1$ mm as the initial radius and $\tau _0 = 176.5$ s as the lifetime. (b) Drop shape for some radii $R\in [0.9 ;4.5]$ mm readable in (a) (see coloured dots). Simulated shapes (full lines), obtained by numerically integrating (2.2), are compared with experimental ones (dotted lines), obtained for side-viewed drops pinned with a needle. The scale is indicated with the capillary length $\ell _c$, that is on the order of 2.5 mm for water. Surface flows, viewed from the top, successively self-organize into (c) six counter-rotating cells (mode $m=3$), (d) four counter-rotating cells (mode $m=2$) and eventually a unique rolling cell (mode $m=1$). Drop keeps on rolling but eventually stops. The consecutive snapshots are extracted from the supplementary movie SM1, available at https://doi.org/10.1017/flo.2022.5. ( f–h) The horizontal median cut views of the most unstable mode from the numerical stability analysis for radii corresponding to (c–e). The radii for the inner flow symmetry transitions from experiments and from stability analyses are plotted is (a) as black and orange lines, respectively.

2.2 Leidenfrost drop shapes

A consequence of the time scale separation is that a given volume of liquid adopts the static shape of a non-wetting drop (Reference Roman, Gay and ClanetRoman, Gay, & Clanet, 2001). If we denote $C(z)$ as the local curvature at a given height $z$ and $C_0$ as the curvature at the apex, the balance of hydrostatic and Laplace pressures can be written $C(z) = C_0 + z/\ell _c^2$. We introduce the curvilinear abscissa $s$, the horizontal radius at a given height $r(z)$ and the angle $\beta$ tangent to the interface. The previous equation can be recast into

(2.2)\begin{equation} \frac{\sin \beta}{r} + \frac{\mathrm{d} \beta}{\mathrm{d}s} = C_0 + \frac{z}{\ell_c^2}. \end{equation}

A numerical integration of (2) for $\beta$ ranging from 0 to $\pi$ and $\mathcal {C}_0\ell _c$ ranging from 0.5 to 10, provides the drop shape for radii ranging from $R=0.9$ to 4.3 mm, plotted as full lines in figure 1(b). These theoretical shapes are found to match the shapes obtained for water drops kept in place by a needle (see dotted lines), except at the drop north pole, where the needle locally forms a meniscus. Increasing the drop size tends to saturate the puddle height $H$ to its maximum value $H_{max} \approx 2 \ell _c$, that is roughly 5 mm for water at 100$\,^\circ \mathrm {C}$. Drops smaller than the capillary length $\ell _c$ are quasi-spherical while drops larger than $\ell _c$ become flattened owing to gravity. As a consequence, evaporation induces geometric changes, particularly on the drop aspect ratio $2R/H$. The presence of strong internal flows could also in principle deform the liquid interface. The Reynolds number associated with the inner flows with typical velocity $V\sim$ cm s$^{-1}$ and kinematic viscosity $\nu$ writes as $Re= R V/\nu$. For a millimetric drop, we have $Re \approx 10^3$, so that inertia overcomes viscosity. The Weber number $We$, which compares inertia with the resisting capillarity is expressed as $We=\rho _0 V^2R/\sigma _0\sim 10^{-2}$. Capillarity therefore outbalances inertia, which justifies that the drop shape does not deviate from the static ones.

2.3 Drop internal flow structure

The apparent quietness of Leidenfrost drops does not reflect what really happens inside the liquid. Side-view PIV measurements performed in a median plane of a water drop containing tracer particles have revealed strong internal flows, with velocities of a few centimetres per second. As the drop shrinks owing to evaporation, Leidenfrost flows undergo a series of successive symmetry breakings. This is further evidenced by focusing on the drop top surface. A Leidenfrost drop, kept on a concave substrate, is seeded with surface particles that have a greater affinity for the air interface. These hollow glass beads are (i) pre-dispersed in water, (ii) skimmed from the interface (where particles accumulate owing to a wetting or shape defect) and (iii) introduced into a pre-dispensed Leidenfrost drop. Despite the apparent axisymmetry of the experiment, interfacial flow structures emerge, as illustrated by the top views in figure 1(c–e) and visible in the supplementary material and movie SM1). When the drop has a radius $R>2.5$ mm, it hosts multiple vortices, which can be described by the azimuthal wavenumber $m \geq 3$ in a cylindrical coordinate system, by denoting ($r,\, z,\, \theta$) and using a periodic wavenumber expansion in $\theta$ as ${\rm e}^{{\rm i} m \theta },\, m\in \mathbb {N}^*$. Within the range $R=[1.8;2.5]$ mm, a mode $m=2$ clearly appears, while for smaller radius $R<1.5$ mm, the droplet rolls in an asymmetric fashion, corresponding to a mode $m=1$. Inner flows thus successively self-organize into six counter-rotating cells ($m=3$, figure 1c); four counter-rotating cells ($m=2$, figure 1d) and eventually a unique rolling cell ($m=1$, figure 1e). We can notice that, at some instance in the supplementary material and movie SM1, the convective cells loose their organization and coherence. The flow structures seem to be transiently perturbed by the drop oscillation in the slightly curved well. However, the dominant modes reappear within a second as a hint of the robustness of the unstable modes. The transition from $m=2$ to $m=1$ can also be visualized by side views, using PIV techniques, as in Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and QuéréBouillant et al. (Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and Quéré2018b) (images reproduced in figures 6a and 5a) or even indirectly measured, as the drop acceleration $a$ (extracted from top views, for water drops initially at rest) suddenly increases with the mode switching to $m=1$. The experiment proposed by Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and QuéréBouillant et al. (Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and Quéré2018b) is reproduced in the supplementary material, for plate temperatures ranging from $250\,^\circ \mathrm {C}$ to $450\,^\circ \mathrm {C}$. The accelerations $a$ of approximately 80 drops as a function of their radii $R$ exhibit similar jumps from $\sim$1 mm s$^{-2}$ (in the regime where drops are flattened by gravity), up to 60 mm s$^{-2}$ as $R$ decreases to $\sim 1$ mm. This indeed corresponds to entering the self-rotation and propelling regime (Reference Bouillant, Mouterde, Bourrianne, Clanet and QuéréBouillant et al., 2018a). The onset for the transition from $m=2$ to $m=1$, referred to as $R_{2\rightarrow 1}$, is found to weakly depend on the plate temperature and to be within the interval $[1; 1.5]$ mm (see § S6 of the supplementary material).

2.4 Temperature gradient at the drop interface

To test the aforementioned thermally based instability scenario, we need to specify the thermal boundary conditions at the drop surface. The drop base is expected to be maintained at roughly the water boiling point ($T_b = 100\,^\circ \mathrm {C}$), while its apex is cooled down by the ambient air. The temperature field in a Leidenfrost drop is measured using an infrared camera (FLIR A600 series), calibrated on the temperature range $[-40 ;+ 150]\,^\circ \mathrm {C}$, only suitable to see water, and not the brass substrate. Water being opaque to infrared wavelengths, this measurement provides the ‘skin’ temperature $T$. Figure 2 shows that, at $t = 0$, that is when $R = 3.5$ mm, $T$ linearly decreases with height $z$ (in millimetres) as $T = -3.75z + 97.0\,^\circ \mathrm {C}$, reaching a maximum $T_{max} = 97.0\,^\circ \mathrm {C}$ at the drop base ($z = 0$). A temperature difference ${\rm \Delta} T\approx 25\,^\circ \mathrm {C}$ thus develops along the drop. At $t = 16$ s, the gradient is slightly smaller with $T = -3.53z + 97.0\,^\circ \mathrm {C}$, and it keeps decreasing as $R$ decreases. These observations are confirmed by introducing a thermocouple inside the liquid (see figure S1 in the supplementary material). For $t = 74$ s, when $R = 1.8$ mm, $T$ suddenly becomes homogeneous with $T \approx 88\,^\circ$C. As best visible in the supplementary material and movie SM1, this coincides with the moment where the flow switches to a symmetry $m=1$. At this transition, the drop starts to vibrate, which enhances mixing. For $t > 100$ s, $T$ thus becomes roughly homogeneous, with a minimum at the drop centre (along the rolling axis), and a maximum at its periphery since liquid is periodically brought close to the hot plate. In this rolling state, the liquid temperature is $T \approx 80\,^\circ$C, a value smaller than the boiling point of water, as consequences of (i) the intensifying evaporation-driven cooling; (ii) a reduction of the flattened area at the drop base from which water is heated, which scales as $R^4/l_c^2$ (Reference Mahadevan and PomeauMahadevan & Pomeau, 1999); (iii) the temperature homogenization due to the rolling-enhanced mixing. We now try to link the existence of such temperature distributions to the origin and structure of the internal flows.

Figure 2. (a) Infrared side views of an evaporating Leidenfrost water drop deposited on a slightly curved surface of brass heated at 350 $^{\circ }$C. Images, taken with a thermal camera give access to the surface temperature (calibration range from $-40\,^{\circ }$C to $+150\,^{\circ }$C to focus on the water surface). Images are extracted from supplementary material and movie SM2. (b) Surface temperature $T$ of a given water drop along its central vertical axis $z$ (white dotted line), showing the change of temperature profile as the drop radius $R$ decreases.

3.1 Problem formulation

Based on the experimental observations, we develop a minimal model assuming that (i) the liquid adopts at any time the static shape of a non-wetting drop; (ii) the base flow inside a drop is steady; (iii) the temperature at the bottom is fixed to the liquid boiling temperature; (iv) the interaction between the liquid and the surrounding gas is not solved completely but modelled by heat transfer correlation laws applied at the side boundary. We sketch in figure 3(a) the problem, where we represent from the side a static drop provided by (2.2). We denote by $\Omega$ the liquid domain, and by $\partial {\varOmega}$ the boundaries, which are decomposed into $\partial {\varOmega} _F$, the upper free interface and $\partial {\varOmega} _S$, the bottom interface of the drop. We also introduce $\partial {\varOmega} _A$ the vertical centreline of the drop, the axis of symmetry. Hence, the numerical computation is done only on the half-domain $r=[0; R]$, as illustrated in figure 3(b). The bottom interface $\partial {\varOmega} _S$ is assumed to be isothermal at temperature $T_s=100\,^\circ \mathrm {C}$ (the phase change of a pure body occurs at a given constant temperature), while the temperature on the side of the drop $\partial {\varOmega} _F$ needs to be evaluated using the heat transfer balance.

Figure 3. Sketch of the problem. (a) The drop (domain ${\varOmega}$) presents boundaries $\partial {\varOmega}$ including the upper free surface $\partial {\varOmega} _F$, the bottom interface $\partial {\varOmega} _S$ and the axis of symmetry $\partial {\varOmega} _A$. (b) Thermal conditions.

3.2 Governing equations

Under the Boussinesq approximation, the governing Navier–Stokes equation for the velocity fields $\boldsymbol {u}=[u_r,\,u_\theta,\,u_z]^{\rm T}$ and the temperature $T$ in cylindrical coordinates $(r,\,\theta,\,z)$ defined in the domain ${\varOmega}$ of boundaries $\partial {\varOmega}$ reads

(3.1a)$$\begin{align}{\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} = {-}\frac{1}{\rho_0} \boldsymbol{\nabla} p - \frac{\rho(T)}{\rho_0} \boldsymbol{g} + \nu \nabla^2 \boldsymbol{u}} \quad &{\rm in}\ {\varOmega}, \end{align}$$

(3.1b)$$\begin{align}{\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0} \quad &{\rm in}\ {\varOmega}, \end{align}$$

(3.1c)$$\begin{align}{\frac{\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa \nabla^{2} T} \quad &{\rm in}\ {\varOmega}, \end{align}$$

where $\boldsymbol {g}$ is the gravitational acceleration in the $z$ direction, $\nu$ the liquid kinematic viscosity and $\kappa$ its thermal diffusivity. The boundary conditions for $\partial {\varOmega} _F,\, \partial {\varOmega} _S$ and $\partial {\varOmega} _A$ are, respectively,

(3.2a)$$\begin{align} {\boldsymbol{S} \boldsymbol{\cdot} \boldsymbol{n} = \boldsymbol{S}_{{gas}} \boldsymbol{\cdot} \boldsymbol{n}-\sigma(\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{n}) \boldsymbol{n} + (\boldsymbol{I}-\boldsymbol{nn}) \boldsymbol{\cdot} \boldsymbol{\nabla} \sigma} \quad &{\rm on} \ \partial{\varOmega}_F, \end{align}$$

(3.2b)$$\begin{align} {\boldsymbol{u\boldsymbol{\cdot} n} = 0}\quad & {\rm on } \ \partial{\varOmega}_F, \end{align}$$

(3.2c)$$\begin{align} {k \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\nabla} T = {-}h_{c} (T-T_{a}) - n' \mathcal{L}}\quad & {\rm on}\ \partial{\varOmega}_F, \end{align}$$

(3.2d)$$\begin{align}{T= T_0}\quad & {\rm on} \ \partial{\varOmega}_S, \end{align}$$

with suitable boundary conditions on the flow axis of symmetry, $\partial {\varOmega} _A$, detailed in § 3.3 for linear perturbations and reading ${u}_r={u}_\theta = 0$ for the base axisymmetric case. Here, $\boldsymbol {S}$ and $\boldsymbol {S}_{{gas}}$ are the stress tensors in the liquid and the gas, respectively, $\sigma$ is the surface tension, $\boldsymbol {n}$ the normal vector, $k$ the conductivity, $h_c$ the convective heat transfer coefficient in air, $n'$ the evaporation rate and $\mathcal {L}$ the latent heat. Equation (3.2c) expresses the heat flux balance at the interface stemming from the liquid ($\,j_k$) transferred into air ($j_c$) and into evaporative cooling ($j_e$), as schematized in figure 3(b).

Both the liquid density $\rho$ and the surface tension $\sigma$ are assumed to vary linearly with the temperature: $\rho (T) = \rho _0 (1-\beta _w(T-T_0)),\, \ \sigma = \sigma _0 - \sigma _1 (T-T_0) ,$ where $\beta _w$ is the water thermal expansion coefficient, $\beta _w= {\rho _0}^{-1} {\partial \rho }/{\partial T}$ and $\sigma _1$ the surface tension variation with the temperature, $\sigma _1= -{\partial \sigma }/{\partial T}$. We decompose now the temperature $T$ as

(3.3)\begin{equation} T = T_0 - \varTheta(r,z), \end{equation}

where $T_0$ is the water boiling temperature taken as reference ($T_0=100\,^\circ \mathrm {C}$) at $z=0$ and $\varTheta$ the deviation from $T_0$. The reference density $\rho _0$ and the surface tension $\sigma _0$ are defined at $T_0$. With these definitions, the density and the surface tension write

(3.4a,b)\begin{equation} \rho(T) = \rho_0 (1-\beta_w \varTheta ) , \quad \sigma = \sigma_0 - \sigma_1 \varTheta. \end{equation}

Denoting by $p_0$ the hydrostatic pressure, the pressure $p$ inside the liquid decomposes as $p = p_0(z) + p_1(r,\,z)$. The time, velocity and length are scaled with $H^2/\kappa,\, \kappa /H$ and $H$, respectively, where $H$ is the drop height. The temperature is scaled with the temperature difference ${\rm \Delta} T(=|\varTheta _{z=0}-\varTheta _{z=H}|)$ and the pressure is scaled with $\rho _0 \nu \kappa / H^2$, which naturally appears when plugging (3.5a–g) in (3.1a), and comparing the viscous term with the pressure term

(3.5a–g)\begin{equation} r = H\hat{r},\quad z = H \hat{z},\quad t = \frac{H^2}{\kappa} \hat{t}, \quad \boldsymbol{u} = \frac{\kappa}{H} \hat{\boldsymbol{u}}, \quad {\varTheta} = {\rm \Delta} T \hat{\varTheta}, \quad p = \frac{\rho_0 \nu \kappa}{H^2} \hat{p}, \quad \boldsymbol{\nabla} = \frac{1}{H}\hat{\boldsymbol{\nabla}}. \end{equation}

We provide in the supplementary material Table S1 with the physical parameters relative to water, vapour, air and their interface relevant to describing our problem. We introduce the dimensionless numbers, Prandtl, Rayleigh, Marangoni, Biot and Sherwood numbers as

(3.6a–e)\begin{equation} {Pr} = \frac{\nu }{\kappa}, \quad {Ra} = \frac{\beta_w g {\rm \Delta} T H^{3}}{\nu \kappa}, \quad Ma = \frac{\sigma_1 {\rm \Delta} T H}{\rho_0 \nu \kappa}, \quad Bi = \frac{h_{c} H}{k}, \quad Sh = \frac{h_{m} H}{D_{va}}.\end{equation}

The governing (3.1) can be now recast into

(3.7a)$$\begin{align} {Pr^{{-}1} \left( \frac{\partial \hat{\boldsymbol{u}}}{\partial \hat{t}} + \hat{\boldsymbol{u}} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} \hat{\boldsymbol{u}} \right) = {-}\hat{\boldsymbol{\nabla}} \hat{p}_1 + Ra \hat{\varTheta} \boldsymbol{e}_z + \hat{\nabla}^2 \hat{\boldsymbol{u}}} \quad &{\rm in} \ {\varOmega}, \end{align}$$

(3.7b)$$\begin{align} {\hat{\boldsymbol{\nabla}} \boldsymbol{\cdot} \hat{\boldsymbol{u}} = 0}\quad & {\rm in} \ {\varOmega}, \end{align}$$

(3.7c)$$\begin{align}{\frac{\partial \hat{\varTheta}}{\partial t} + \hat{\boldsymbol{u}} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}}\hat{\varTheta} = \hat{\nabla}^{2}\hat{\varTheta}}\quad & {\rm in} \ {\varOmega}, \end{align}$$

with the boundary conditions:

(3.8a)$$\begin{align} {-\hat{p}_1 \boldsymbol{n} + (\hat{\boldsymbol{\nabla}} \hat{\boldsymbol{u}} + (\hat{\boldsymbol{\nabla}} \hat{\boldsymbol{u}})^T)\boldsymbol{\cdot} \boldsymbol{n} = {-} Ma ( \boldsymbol{I}- \boldsymbol{nn})\boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} \hat{\varTheta}} \quad & {\rm on} \ \partial{\varOmega}_F, \end{align}$$

(3.8b)$$\begin{align}{\hat{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{n} = 0} \quad &{\rm on} \ \partial{\varOmega}_F\ {\rm and}\ \partial{\varOmega}_S, \end{align}$$

(3.8c)$$\begin{align}{\boldsymbol{n} \boldsymbol{\cdot} \hat{\boldsymbol{\nabla}} \hat{\varTheta} = {-}Bi \left( \hat{\varTheta} + \frac{T_0 -{T}_a}{{\rm \Delta} T} \right) - \frac{n'\mathcal{L}H}{k}} \quad &{\rm on} \ \partial{\varOmega}_F, \end{align}$$

(3.8d)$$\begin{align}{\hat{\varTheta} = 0} \quad &{\rm on} \ \partial{\varOmega}_S, \end{align}$$

and suitable symmetry conditions on $\partial {\varOmega} _A$ (see § 3.3 for details). The evaporation rate $n'$ is linked to the $Sh$ number, which represents the evaporative heat transfer coefficient. The heat exchanges on $\partial {\varOmega} _F$, which determine $Bi$ and $Sh$ are modelled using the Ranz–Marshall correlation (Reference Incropera, DeWitt, Bergman and LavineIncropera, DeWitt, Bergman, & Lavine, 1996; Reference Ranz and MarshallRanz & Marshall, 1952) as detailed in the supplementary material (§ S2), where the input ambient air temperature $T_a$ is also measured and provided in § S4.

3.3 Stability analysis

The steady toroidal base flow ($\boldsymbol {u}_b,\, p_b,\,\varTheta _b$) is obtained by finding the nonlinear steady state solution of (3.1) satisfying the boundary condition (3.2) using the Newton method. The base flow is solved using the dimensional equations since some dimensionless numbers depend on ${\rm \Delta} T$, which is also an unknown of the thermal problem. (Note that the non-dimensional equation can be resolved by defining $Ra$ and $Ma$ with the known parameters, i.e. $T_s$. However, we kept the classical definitions of $Ra$ and $Ma$ which are with ${\rm \Delta} T$.) We first compute the base flows, estimate ${\rm \Delta} T$ and deduce the dimensionless numbers. This approach differs from thermoconvective studies on a flat plate or in a cylinder. In these situations, there is no radial pressure gradient and the vertical gradient is solely balanced with the density, without inducing any velocity field. The temperature difference ${\rm \Delta} T$ and thermal dimensionless parameters are control parameters. In contrast, when the radial pressure gradient is non-zero (as in the Leidenfrost configuration, owing to the boundary conditions), it induces a steady flow with non-zero velocity. Both ${\rm \Delta} T$ and the dimensionless parameters are solutions of the problem. Assuming infinitesimal perturbations on the base flow with the complex frequency ${\varOmega}$ and the azimuthal wavenumber $m$, the dimensionless flow field, pressure and temperature are decomposed using the normal mode expansion

(3.9)\begin{equation} [\hat{u}_r,\hat{u}_\theta,\hat{u}_z,\hat{p}_1,\hat{\varTheta}](r,\theta,z,t) = [{u}_r,{u}_\theta,{u}_z,{p},{\varTheta}](r,z) \exp({\rm i}m\theta - {\rm i} {\varOmega} t ) + {\rm c.c.}, \end{equation}

where c.c. indicates complex conjugate. The linearized (3.7) becomes

(3.10a)$$\begin{align} {Pr^{{-}1} ( {\rm i} {\varOmega} \boldsymbol{u} + \boldsymbol{u}_b \boldsymbol{\cdot} \nabla_m + \boldsymbol{u} \boldsymbol{\cdot} \nabla_0 \boldsymbol{u}_b ) = {-}\nabla_m {p} + Ra \varTheta \boldsymbol{e}_z + \nabla^2_m \boldsymbol{u}} \quad &{\rm in} \ {\varOmega}, \end{align}$$

(3.10b)$$\begin{align}{{\nabla}_m \boldsymbol{\cdot} \boldsymbol{u} = 0} \quad & {\rm in} \ {\varOmega}, \end{align}$$

(3.10c)$$\begin{align}{{\rm i} {\varOmega} \varTheta + \boldsymbol{u}_b \boldsymbol{\cdot} {\nabla_m} \varTheta + \boldsymbol{u} \boldsymbol{\cdot} {\nabla_0} \varTheta_b = {\nabla}_m^{2}\varTheta} \quad & {\rm in}\ {\varOmega}, \end{align}$$

where $\nabla _m$ represents the derivative in $\theta$ is replaced by ${\rm i}m$. The boundary conditions on $\partial {\varOmega} _F$ are

(3.11a)$$\begin{align} {-{p} \boldsymbol{n} + ({\nabla}_m \boldsymbol{u} + ({\nabla}_m \boldsymbol{u})^T)\boldsymbol{\cdot} \boldsymbol{n} = {-} Ma ( \boldsymbol{I}- \boldsymbol{nn})\boldsymbol{\cdot} {\nabla}_m {\varTheta}} \quad & {\rm on}\ \partial{\varOmega}_F, \end{align}$$

(3.11b)$$\begin{align}{\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} = 0} \quad & {\rm on} \ \partial{\varOmega}_F, \end{align}$$

(3.11c)$$\begin{align}{\boldsymbol{n} \boldsymbol{\cdot} {\nabla}_m {\varTheta} = {-}Bi {\varTheta}} \quad & {\rm on}\ \partial{\varOmega}_F. \end{align}$$

Note that the perturbation temperature is now only affected by the $Bi$ as the other heat transfer coefficients are constant and contribute only in the zero-order base flow equation. One could also linearize the last term in (3.2c) and include it in the stability analysis, but the effect on the stability results is negligible (Reference Yim, Bouillant and GallaireYim, Bouillant, & Gallaire, 2021). The condition prescribed on the drop axis of symmetry $\partial {\varOmega} _A$ (illustrated in figure 3) depends on $m$, in particular on the mode symmetry. The following conditions are thus used as in Reference Batchelor and GillBatchelor and Gill (Reference Batchelor and Gill1962):

(3.12a)$$ m = 0, \quad \quad {u_r = u_\theta = \frac{\partial u_z}{\partial r} = \frac{\partial \varTheta}{\partial r} = 0}\quad\quad {\rm on} \ \partial {{\varOmega}_A},$$

(3.12b)$$ m = 1, \quad\quad {u_z = \varTheta = p = \frac{\partial u_r}{\partial r} = 0 = \frac{\partial u_\theta}{\partial r} = 0} \quad \quad {\rm on} \ \partial{{\varOmega}_A},$$

(3.12c)$$ m \geq 2, \quad\quad {u_r = u_\theta = u_z = \varTheta = p = 0}\quad\quad {\rm on}\ \partial{{\varOmega}_A}.$$

Finally, we prescribe on $\partial {\varOmega} _S$ a Dirichlet condition for the temperature and a free slip condition for the velocity since the bottom surface is not in contact with the plate

(3.13a)$$\begin{align} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{n} = 0 \quad {\rm on} \ \partial{\varOmega}_S, \end{align}$$

(3.13b)$$\begin{align}\varTheta = 0 {\rm on} \ \partial{\varOmega}_S. \end{align}$$

3.4 Numerical methods

All the numerical analyses are performed using FreeFEM++ software (Reference HechtHecht, 2012) for axisymmetric cylindrical coordinates $(r,\, z)$. The velocity, pressure and temperature are discretized with Taylor–Hood P2, P1 and P2 elements, respectively. The typical number of triangles is ${\sim}10^4$. The linear equations and the eigenvalue problem are solved using UMFPACK library and the ARPACK shift–invert method, respectively. Starting with the initial radius $R=3.5$ mm, the nonlinear solution of (3.1) is solved for the given water properties. Once the solution for one radius is found, it is used as the initial guess for the smaller radius. Within five iterations, the root-mean-squared (L2) norm residual becomes smaller than $1 \times 10^{-8}$. The zero normal velocity condition on the free surface is applied using the Lagrange multiplier method (Reference BabuškaBabuška, 1973; Reference Yim, Bouillant and GallaireYim et al., 2021).

4.1 Pure buoyancy-induced flow ($Ma=0$)

4.1.1 Base flow ($Ma=0$)

Let us first consider the case of pure buoyant flows, neglecting thermocapillary effects. The base flow is thus computed following (3.1) while setting the superficial stress on $\partial {\varOmega} _F$ to zero. We restrict this parametric study in $R$ to the limit $R \lesssim 2$ mm, since, according to infrared measurements (figure 2), the drop surface temperature tends to homogenize, suppressing Marangoni surfaces flows. Figure 4 shows the base flow obtained in this limit of $Ma=0$, for drop radii $R=2$ mm, $R=1.3$ mm and $R=0.8$ mm, respectively. In the absence of a surface tension gradient, the base flow exhibits pure thermal convection: the flow rises along the centre axis, which is warmer since it is insulated from the drop interface, and descends along the side of the drop, where the drop is cooled. The inner velocities for purely buoyant flows (${\sim}1$ cm s$^{-1}$) underestimate the experimental observation (${\sim}5$ cm s$^{-1}$) of a similar size of drop (see figure 5(a,c) for the experimental measurements).

Figure 4. Base flows for $Ma=0$: (a) $R=2$ mm, (b) $R=1.3$ mm and (c) $R=0.8$ mm. The temperature and velocity fields are shown in colour and with arrows. (d) Dominant growth rate ${\varOmega} _i$ of different azimuthal modes with decreasing drop radius $R$.

Figure 5. (a,c) Velocity fields within a droplet with $R \sim 0.9$ mm deduced from PIV measurements. (b,d) Corresponding flow fields deduced from the numerical stability analysis in the absence of Marangoni effects ($m=1,\, Ma=0,\, Ra=6.1 \times 10^3$). Velocity arrows are plotted within the lateral (a,b) and horizontal (c,d) planes. The red dashed line in (c) indicates the area of the flattened base of the drop on which the bottom camera focuses. Colour in (b,d) indicates the perturbed temperature field (normalized with its maximum value).

4.2 Reduced Marangoni approximation

A major challenge for computing the base flow in a Leidenfrost drop is to be able to predict the surface tension distribution. When we take the exact surface tension temperature dependence $\sigma _1$ provided in Table S1 (supplementary material), the Marangoni flows are largely overestimated compared with experimental observations (see figure S3 of the supplementary material). It is known, however, that Marangoni effects are very often markedly reduced (Reference Dhavaleswarapu, Migliaccio, Garimella and MurthyDhavaleswarapu, Migliaccio, Garimella, & Murthy, 2010; Reference Hu and LarsonHu & Larson, 2005a, Reference Hu and Larson2006) or even almost absent (Reference Dash, Chandramohan, Weibel and GarimellaDash et al., 2014; Reference Gelderblom, Bloemen and SnoeijerGelderblom, Bloemen, & Snoeijer, 2012; Reference Marin, Gelderblom, Lohse and SnoeijerMarin, Gelderblom, Lohse, & Snoeijer, 2011). Based on measured quantities, we first try to correct the surface tension temperature dependence $\sigma _1$. To that end, we evaluate the Rayleigh $Ra$ and Marangoni $Ma$ numbers, which compare the buoyancy and surface tension stresses in comparison with inertia, respectively, as defined in (3.6a–e). Using the parameters documented in Table S1 of the supplementary material, they are $Ra \sim 1.8 \times 10^5$ and $Ma \sim 3.6 \times 10^5$ for ${\rm \Delta} T = 25\,^\circ \mathrm {C},\, H=4$ mm, $R=3.5$ mm. Both values greatly exceed the expected critical values for the onset for Rayleigh–Bénard and Marangoni instabilities (see Table S2 in the supplementary material) – typically $Ra_c \sim O(10^3)$ and $Ma \sim O(10^2)$, for similar geometries (Reference ChandrasekharChandrasekhar, 1961). As detailed in § S7 of the supplementary material, the effective surface tension variation $\sigma _{1,eff}$ (and thus $Ma_{eff}$) is determined using numerical analysis. We select the $\sigma _{1,eff}$ that best describes the temperature difference ${\rm \Delta} T$ within the liquid (figure 2), as well as the flow velocities, typically 5 cm s$^{-1}$. The temperature profiles reported in figure 2 (left panels, $R=3.5$ mm and $R=3.0$ mm) are best represented by the curve with $\sigma _{1,eff}=4 \times 10^{-5} \sigma _1$ as shown in Figure S7 of the supplementary material for the surface temperature and figure 6 for the velocity field. The surface stress seems to be reduced by a few orders of magnitude. This has also been reported in Reference Savino, Paterna and FavaloroSavino, Paterna, and Favaloro (Reference Savino, Paterna and Favaloro2002), Reference Hu and LarsonHu and Larson (Reference Hu and Larson2002), Reference Hu and LarsonHu and Larson (Reference Hu and Larson2005b) and Reference Hu and LarsonHu and Larson (Reference Hu and Larson2006), where it is ascribed to surface contamination, or to the fact that at a large Marangoni surface stress, it is moderated by dissipation or by transport-limited properties of the fluid, reducing the achievable velocities (see figure 6 of the supplementary material). In the following stability analysis, we thus use this reduced surface tension variation by adjusting $\sigma _1$ to $\sigma _{1,eff}=4 \times 10^{-5} \sigma _1$, which is the only tuneable parameter of our study (all other physical parameters are kept to exact values of the given thermal properties, provided in Table S1).

Figure 6. (a) PIV measurements in a drop with $R=2.5$ mm (Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and QuéréBouillant et al., 2018b). (b) Base flow for $Ma_{eff}=11.3$ ($\sigma _{1,eff} = 4\times 10^{-5} \sigma _1$). The colour map and arrows give the inner temperature and velocity.

4.3 Thermocapillary flow

4.3.2 Stability analysis with effective $Ma$

Figure 7 shows the dominant eigenvalues as a function of decreasing radius (with $\sigma _{1,eff}=4\times 10^{-5} \sigma _{1}$). Compared with the pure buoyant flow, there exist several unstable azimuthal modes ranging from $m=1$ to 3. Figure 7(a,c,e) shows the two least stable growth rates ${\varOmega} _i$ and their corresponding frequencies ${\varOmega} _r$ are shown in figure 7(b,d, f) for azimuthal wavenumbers $m=1,\,2$ and 3. For all modes $m=1,\,2,\,3$, there exist two branches of unstable modes: one with ${\varOmega} _r=0$ (steady) and other with ${\varOmega} _r \neq 0$ (unsteady). The $m=1$ mode (figure 7a,b) is unstable both at large radius $R>2.8$ mm and small radius $R<1.3$ mm, but it shows a window of stability for intermediate values of the radius. The prevailing unstable mode for $R>2.8$ mm is unsteady, with frequency ${\varOmega} _r \,{\sim}\,10\ {\rm s}^{-1}$, while for $R<1.3$ mm, it becomes steady. Mode $m=2$ (figure 7c,d) is only unstable in the intermediate radius range $R=[2.3 ; 1.3]$ mm. At radius $R\sim 2$ mm, the unsteady branch dominates, yet with a small frequency, but the steady branch takes over at smaller $R$. Mode $m=3$ (figure 7e, f) is unstable when $R<2.8$ mm. Its steady branch dominates at large $R$ but it becomes unsteady as the radius gets closer to $R\approx 2.8$ mm. The unsteady branch of mode $m=3$ reaches a maximum growth rate around $R=2.4$ mm.

Figure 7. Growth rates and the corresponding frequencies for (a,b) $m=1$, (c,d) $m=2$ and (e, f) $m=3$ for the two least unstable modes as a function of decreasing radius. The dominant modes for each $m$: (g) growth rate and (h) corresponding frequency.

Finally, figure 7(g,h) collects and retains only the most unstable modes $m$ and their corresponding frequencies. We note in figure 7(g) that the axisymmetric mode $m=0$ is always stable (${\varOmega} _i<0$). The dependence of the dominant azimuthal mode on radius thus becomes explicit: for $R > 2.1$ mm, the mode $m=3$ is the most unstable (first steady and then oscillatory). As the radius decreases further, the $m=3$ mode becomes stable around $R=2$ mm. At $R=2.2$ mm, the mode $m=2$ starts to grow and becomes the only unstable mode in the range of radius $R=[2; 1.3]$ mm, with the maximum growth rate at $R=1.6$ mm. For $R<1.3$ mm, the steady mode $m=1$ takes over as the dominant unstable mode. This trend is very similar to the apparent mode transition in experiments, which is $m\geq 3$ for $R \gtrsim 2.8$ mm, $m=2$ for $R \approx [2.5 ; 1.8]$ mm and $m=1$ for $R \approx 1.5$ mm, as shown in figure 1. A growth rate ${\varOmega} _i \sim 1~{\rm s}^{-1}$ implies that within 2.3 s, the perturbation amplitude becomes 10 times larger than its initial value. For all modes, we verify that the growth rate is much faster than the drop evaporation rate ${\varOmega} _i \gg (1/R){\rm d}R/{\rm d}t$, underlying the quasi-steady assumption in § 2.1.

Figure 8 illustrates the eigenvectors of the most unstable mode at radii $R=3,\,\ 1.5$ and 1 mm. The top panels display side views while the right ones show top views: $xy$ plane cut at the maximum radius $R=R_{max}$. The colours indicate temperature perturbations and the arrows are the in-plane velocity perturbations. For $R=3$ mm (figure 8a,d), the mode $m=3$ dominates. Temperature perturbations are localized along the interface and the central axis of the drop while velocity perturbations are localized where the temperature perturbations are the weakest. The top view shows three counter-rotating vortex pairs, very similar to our observations (figure 1c). For $R=1.5$ mm (figure 8b,e), the mode $m=2$ becomes the most unstable thus dominant mode. Temperature perturbations are maximum near the central axis, while velocity perturbations remain localized where temperature perturbations are weak. The top view shows four vortex cells (or two vortex pairs), in close agreement with the experiments (figure 1d). For $R=1$ mm (figure 8c, f), the displacement mode $m=1$ is the most unstable. Temperature perturbations are null on the central axis, where velocity perturbations are maximum. However, the flow structure does not match the solid-like rolling motion reported in Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and QuéréBouillant et al. (Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and Quéré2018b) and illustrated in the top of figure 1(e). We interpret this at a consequence of the temperature homogenization evidenced at small $R$ in figure 2. The temperature difference at the liquid surface indeed tends to vanish as $R$ becomes smaller than 1.8 mm. We thus expect thermocapillary effects to weaken and eventually vanish. We provide in Figure S8 (supplementary material), the (reduced) Marangoni and Rayleigh numbers with $R$, showing that both effects weaken as the drop shrinks in size, with a more abrupt decrease of thermocapillary effects than thermobuoyant effects, prompting us to use our predictions for the case $Ma=0$ when $R$ becomes sufficiently small. The observed rolling motion is thus better captured by the eigenmode for $Ma=0$ as shown in figure 5 than the one with $Ma \sim O(10)$ in figure 8(e).

Figure 8. Eigenmodes of the most unstable modes at radii (a,d) $R=3$ mm ($Ma_{eff}=12,\, Ra=1.5\times 10^5$), (b,e) $R=1.5$ mm ($Ma_{eff}=9.6,\, Ra=5\times 10^4$) and (c, f) $R=1$ mm ($Ma_{eff}=8,\, Ra=2.4\times 10^4$). Colours indicate temperature perturbations and the arrows show the in-plane velocity perturbations normalized with their maximum real values. The top view is a plane cut at the maximum radius $R=R_{max}$.