2. Formulation of the computational model

Consider a Leidenfrost drop placed gently on a uniformly heated rigid horizontal surface at a constant $T_{w}$, where $T_{w}$ is kept above $T_{L}$ ($\gg T_{b}$), see figure 1(a). The drop levitates due to the evaporation-driven lubrication pressure which develops in the draining vapour film between the drop surface and the hot rigid wall. When the lubrication force due to vapour pressure in the gap is equal to the weight of the liquid drop, the drop is at equilibrium, as can be seen in figures 1(b) and figure 1(c). We aim to determine the quasi-static equilibrium shape of such a drop and, in particular, the underlying vapour film geometry.

Figure 1. Schematic of an axisymmetric Leidenfrost drop above an isothermal hot rigid flat surface at temperature $T_{w}$. (a) An initially spherical drop of radius $R_0$ placed on a vapour cushion at an initial height $h_{0}$ above a heated flat surface with cylindrical coordinates ($r,z,\theta$) shown, (b) shows an experimental image of a quasi-static Leidenfrost drop floating on a thin vapour film above a flat surface, taken from Quéré (Reference Quéré2013) and (c) a sketch of a quasi-static Leidenfrost drop levitating on a vapour layer of thickness $h(r)$. Numerical solutions of the theoretical model by Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) (see § 2) are obtained by patching the solution for the upper surface of the drop to that for the lower surface bordering the lubrication vapour layer, at $r=R_{p}$. The image shows the maximum droplet radius $R_{max}$, the radius of the neck $R_{neck}$ and the height of the neck measured from the solid surface (minimum vapour film thickness) $h_{neck}$. The evaporation mass flux across the liquid–vapour interface is denoted by $j$.

We consider an axisymmetric problem, an assumption which is discussed later, described in the cylindrical coordinate system $(r, z, \theta )$, where $r$ is the radial coordinate, $z$ is the axial coordinate measured in the opposite direction of gravity and the problem is considered to be independent of the azimuthal coordinate $\theta$. As shown in figure 1, $z=0$ and $r=0$ represent the solid wall and the axis of symmetry, respectively. Even though the liquid and vapour properties are temperature dependent, in this study, for simplicity, we assume they can be approximately evaluated at the boiling temperature $T_{b}$ and at the mean temperature $T_{m}=(({T_{w}+T_{b}})/{2}$) between the wall and liquid, respectively, similar to Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014).

In this section, we formulate our computational model and then, in the next section, show how this relates to the simplified theoretical model of Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014). Our approach, which naturally extends to the consideration of dynamic processes, is to start with an initially spherical liquid drop of radius $R_{0}$ which is released from rest above a solid at a height $h_{0}=0.1R_{0}$, as shown in figure 1(a), falls due to gravity and then evolves dynamically in time towards a quasi-static shape.

2.1. Liquid flow

The flow of liquid inside the drop is governed by the isothermal incompressible Navier–Stokes equations, with incompressibility

(2.1)\begin{equation} \frac{1}{r}\frac{\partial (ru_{l})}{\partial r}+\frac{\partial w_{l}}{\partial z}=0, \end{equation}

and the momentum equations

(2.2)\begin{gather} \rho_{l}\left(\frac{\partial{u_{l}}}{\partial t} + u_{l}\frac{\partial{u_{l}}}{\partial{r}}+w_{l}\frac{\partial{u_{l}}}{\partial{z}}\right) ={-} \frac{\partial{p_{l}}}{\partial r} + \mu_{l} \left(\frac{\partial^{2}u_{l}}{\partial r^{2}}+\frac{1}{r}\frac{\partial u_{l}}{\partial r}+\frac{\partial^{2}u_{l}}{\partial z^{2}}-\frac{u_{l}}{r^{2}}\right), \end{gather}

(2.3)\begin{gather} \rho_{l}\left(\frac{\partial{w_{l}}}{\partial t} + u_{l}\frac{\partial{w_{l}}}{\partial{r}}+w_{l}\frac{\partial{w_{l}}}{\partial{z}}\right) ={-} \frac{\partial{p_{l}}}{\partial z} + \mu_{l} \left(\frac{\partial^{2}w_{l}}{\partial r^{2}}+\frac{1}{r}\frac{\partial w_{l}}{\partial r}+\frac{\partial^{2}w_{l}}{\partial z^{2}}\right)-\rho_{l}{g}, \end{gather}

where henceforth a subscript $l$ denotes the liquid phase and $v$ will denote the vapour. Here, $\rho$ will represent densities; $\mu$ dynamic viscosities; $u$ and $w$ are velocity components in the $r$ and $z$ directions; $p$ is the pressure; and $t$ is the time.

2.2. Vapour flow

In the vapour film we develop a lubrication model, with height $h(r,t)$ considered slowly varying (${\partial h}/{\partial r}\ll 1$) and small compared with the drop size $R_{max}$ ($\epsilon =h/R_{max} \ll 1$) in those parts of the film that contribute significantly to the pressure drop in the film and vapour generation. In addition, we neglect the influence of gravity in the vapour layer and consider inertial forces negligible compared with viscous ones, for which we need $\epsilon Re_{v} \ll 1$ (Aursand, Davis & Ytrehus Reference Aursand, Davis and Ytrehus2018), where the Reynolds number of the thin vapour layer is defined as $Re_{v}=\rho _{v}h^{3}\vert \boldsymbol {\nabla } p_{v}\vert /\mu ^{2}_{v}$ (Biance et al. Reference Biance, Clanet and Quéré2003; Celestini et al. Reference Celestini, Frisch and Pomeau2012; Aursand et al. Reference Aursand, Davis and Ytrehus2018). This establishes the use of lubrication theory in the vapour film (Biance et al. Reference Biance, Clanet and Quéré2003; Sobac et al. Reference Sobac, Rednikov, Dorbolo and Colinet2014; Aursand et al. Reference Aursand, Davis and Ytrehus2018), with incompressibility unchanged

(2.4)\begin{equation} \frac{1}{r}\frac{\partial (ru_{v})}{\partial r}+\frac{\partial w_{v}}{\partial z}=0, \end{equation}

and the momentum equations simplifying to

(2.5a,b)\begin{equation} \frac{\partial p_{v}}{\partial r}= \mu_{v}\frac{\partial^{2} u_{v}}{\partial z^{2}} \quad \mathrm{and}\quad \frac{\partial p_{v}}{\partial z}=0. \end{equation}

The velocity component tangent to the interface ($z=h(r,t)$) is assumed continuous across it, so that $u_{v,\varGamma }=u_{l,\varGamma }\equiv u_{\varGamma }$, where $\varGamma$ indicates a surface property of the liquid–vapour interface and subscripts $v,\varGamma$ and $l,\varGamma$ refer to bulk properties at the interface. By integrating equation (2.5a,b) with the boundary conditions $u_{v}(z=0)=0$ at the solid wall and $u_{v}(z=h(r,t))=u_{\varGamma }$ at the drop surface, the local radial velocity profile in the lubricating film is obtained as

(2.6)\begin{equation} u_{v}(r,z,t)= \frac{1}{2\mu_{v}}\frac{\partial p_{v}(r,t)}{\partial r}\{z^{2}-zh(r,t)\} + u_{\varGamma}(r,t)\frac{z}{h(r,t)}. \end{equation}

This shows that the radial velocity profile underneath the drop is expressed by the superposition of two components, a Poiseuille flow (pressure-driven) component with an immobile interface and a Couette flow driven by the interface's tangential velocity with no pressure gradient.

In the Leidenfrost phenomenon, heat transfer and evaporation both take place in the region of the vapour film. Under the considered lubrication approximation, the effects of convective heat transfer are negligible compared with thermal conduction when $\epsilon Re_{v} Pr_{v}\ll 1$, where the Prandtl number is $Pr_{v}=\mu _{v}c_{p,v}/k_{v}$ (for water vapour $Pr_{v}\approx 0.95$ at $T_{m}=200\,^{\circ }\mathrm {C}$), and $k_{v}$ and $c_{p,v}$ are thermal conductivity and specific heat capacity (see Aursand et al. Reference Aursand, Davis and Ytrehus2018) . Hence, in the lubrication approximation the energy conservation equation simplifies to

(2.7)\begin{equation} \frac{\partial^{2} T_{v}}{\partial z^{2}}=0, \end{equation}

which indicates that the temperature across the vapour film varies linearly between the solid surface temperature $T_{v}\arrowvert _{z=0}=T_{w}$ and the liquid–vapour interface temperature $T_{v}\arrowvert _{z=h}=T_{b}$.

Let the rate of change of the vapour film height in the frame of reference moving horizontally with speed $u_{\varGamma }$ be $w_{\varGamma }$. Then in the laboratory frame

(2.8)\begin{equation} \frac{\partial h}{\partial t}=w_{\varGamma}-u_{\varGamma}\frac{\partial h}{\partial r}. \end{equation}

Without evaporation, the usual kinematic boundary conditions then give $w_{\varGamma }=w_{l,\varGamma }=w_{v,\varGamma }$, where $w_{l,\varGamma }$ and $w_{v,\varGamma }$ are the vertical liquid and vapour velocities at the interface, respectively. In our case, the heat transferred by conduction from the hot solid to the drop is used for evaporative vapour generation at the liquid–vapour interface (as the temperature in the liquid is assumed fixed at its boiling point, so there is no heat flux into it). In this case, mass and energy conservation at $z=h$ yield (keeping in mind that the interface is nearly horizontal)

(2.9)\begin{equation} \rho_{v}(w_{v,\varGamma}-w_{\varGamma})=\rho_{l}(w_{l,\varGamma}-w_{\varGamma})={-}j,\quad \mathrm{and}\quad q_{v}=jL, \end{equation}

where $j$ is the evaporative mass flux through the liquid–vapour interface and $L$ is the latent heat of evaporation. Equation (2.9) implies that the vertical velocities undergo a jump across the interface due to evaporation.

Since $\eta =\rho _{v}/\rho _{l} \sim 10^{-3}$, (2.9) implies that $w_{v,\varGamma }-w_{\varGamma }$ is very large compared with $w_{l,\varGamma }-w_{\varGamma }$; thus, we can assume $w_{\varGamma }=w_{l,\varGamma }$. Furthermore, the low density ratio $\eta$ establishes that the timescales for viscous and thermal diffusion in the vapour film are very small compared with the liquid in the drop, which also marks a quasi-static process in this Leidenfrost phenomenon. In (2.9), the heat flux $q_{v}$ is modelled by Fourier's law expressed as

(2.10)\begin{equation} q_{v}={-}k_{v}\left(\frac{\partial T}{\partial z}\right){\vert}_{z=h}=\frac{k_{v}\Delta T}{h}, \end{equation}

where $\Delta T=T_w-T_b$. Hence, (2.9) with (2.8) at $z=h$ can be rewritten as

(2.11)\begin{equation} w_{v,\varGamma}=\frac{\partial h}{\partial t}+u_{\varGamma}\frac{\partial h}{\partial r}-\frac{k_{v}\Delta T}{L\rho_{v}h}. \end{equation}

By substituting (2.6) into (2.4) and integrating equation (2.4) over the film thickness in the vertical direction with boundary conditions $w_{v}\arrowvert _{z=0}=0$ and $w_{v}\arrowvert _{z=h}=w_{v,\varGamma }$ (obtained from (2.11)), and applying the Leibniz integral rule, we get

(2.12)\begin{equation} \frac{1}{r}\frac{\partial}{\partial r}\left \{r\left(-\frac{h^{3}}{12\mu_{v}}\frac{\partial p_{v}}{\partial r} + u_{\varGamma}\frac{h}{2}\right) \right \}+\frac{\partial h}{\partial t}-\frac{k_{v}\Delta T}{L\rho_{v}h}=0. \end{equation}

Note that, in regions where the lubrication assumptions are not valid, $h$ is large, the last (evaporation) term vanishes, so the expression in braces is constant and then $\partial p_v/\partial r\approx 0$, which is still correct, so the equation remains valid. From (2.12), one readily obtains the pressure distribution $p_{v}(r,t)$ in the film expressed as a differential equation along the liquid–vapour interface

(2.13)\begin{equation} \frac{\partial^{2}p_{v}}{\partial r^{2}}+\frac{\partial}{\partial r}\{\ln{(rh^{3})}\}\frac{\partial p_{v}}{\partial r}=\frac{6\mu_{v}}{rh^{3}}\left \{{-}2r\frac{k_{v}\Delta T}{L\rho_{v}h}+\frac{\partial}{\partial r}\left(rhu_{\varGamma}\right )+2r\frac{\partial h}{\partial t}\right \}. \end{equation}

The corresponding boundary conditions are given by

(2.14)\begin{equation} \frac{\partial p_{v}}{\partial r}{\bigg \vert}_{r=0}=0\quad \mathrm{and}\quad p_{v}\arrowvert_{r=r_{o}}=p_{o}. \end{equation}

Here, $r_{o}$ is defined at the boundary of the vapour film where $p_{o}=0$ is equal to the atmospheric pressure (i.e. pressures are measured relative to their atmospheric value). The exact choice of $r_{o}$ is discussed below. Finally, in the lubrication approximation, the shear stress on the drop surface is given by $\tau _{v}(r,t)=\mu _{v}(\left.{\partial u_{v}}/{\partial z}\right |)_{z=h}$. This leads to

(2.15)\begin{equation} \tau_{v}=\frac{h}{2}\left(\frac{\partial p_{v}}{\partial r}\right)+\mu_{v}\frac{u_{\varGamma}}{h} \end{equation}

with, as expected, terms from both Poiseuille and Couette components.

3. Theoretical model

This section reviews the theoretical framework proposed by Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) and discusses its connections to the computational model just derived. Further details of the numerical implementation are provided in Appendix A. As previously, the surface of an axisymmetric Leidenfrost drop at equilibrium is divided into two parts, separated by a patching point ($R_{p}, h_{p}$), with the upper and lower drop surfaces treated differently, as shown in figure 1(c). We choose the patching point to lie at the position where ${\rm d}\tilde {z}/{\rm d}\tilde {r}=1$. This is different from the choice we have made in our computational model (see § 2); thus, agreement between the two models will also indicate insensitivity to the choice of the patching point. The upper surface of the drop is determined from the Young–Laplace equation and the lower surface is, in addition, affected by a lubrication pressure from the thin vapour layer, with the two solutions matched at the patching point in order to predict the equilibrium shape of the drop.

First, the shape of the upper surface of the drop at equilibrium is governed by a balance between the interfacial capillary pressure and hydrostatic pressure variations in the liquid,

(3.1)\begin{equation} \sigma\kappa-\rho_{l}g z_1=\sigma\kappa_{apex}, \end{equation}

where $z_1=z-z_{top}$ is the vertical coordinate measured from the apex of the drop ($z=z_{top}$) and $\kappa _{apex}$ is the curvature of the drop surface at the apex. This equilibrium condition neglects internal flow in the drop, which is taken into account in our computational model. In dimensionless form (with respect to the capillary length),

(3.2)\begin{equation} \tilde \kappa =\tilde{z}_1 + \tilde{\kappa}_{apex}. \end{equation}

Given $\tilde {\kappa }_{apex}$, (3.2) can be solved to obtain the shape of the upper surface, in the form of the dependences ${\tilde {z}_1}({\tilde {s}})$ and ${\tilde {r}}({\tilde {s}})$, where $\tilde {s}$ is the arc length measured from the apex (Duchemin, Lister & Lange Reference Duchemin, Lister and Lange2005). Rather than specifying the size of the drop $\tilde {R}_{max}$, it is more convenient to choose ${\tilde {\kappa }}_{apex}$ instead and $\tilde {R}_{max}$ is then found from the solution. Note that at this point $\tilde {z}_{top}$ remains unknown, in other words, the shape is determined up to a vertical shift.

Second, for the lower surface we include the vapour pressure, and then the equilibrium condition (again neglecting the liquid flow) is

(3.3)\begin{equation} p_v+\rho_l gh+\sigma\kappa=\mathrm{const.} \end{equation}

We next use (2.12), where we again neglect the liquid flow and therefore $u_{\varGamma }$, which gives

(3.4)\begin{equation} \frac{\partial h}{\partial t}=\frac{k_{v}\Delta T}{L\rho_{v}h}-\frac{1}{r}\frac{\partial}{\partial r}\left \{r\left(-\frac{h^{3}}{12\mu_{v}}\frac{\partial p_{v}}{\partial r}\right) \right \}. \end{equation}

Using (3.3) to eliminate $p_v$ and switching to dimensionless variables, we get

(3.5)\begin{equation} \frac{\partial \tilde{h}}{\partial \tilde{t}}=\frac{E}{\tilde{h}}-\frac{1}{12\tilde{r}}\frac{\partial}{\partial \tilde{r}}\left[\tilde{r} \tilde{h}^{3}\frac{\partial}{\partial \tilde{r}}\left(\tilde{h}+\tilde{\kappa} \right)\right], \end{equation}

where ${\tilde {t}}=t/t_c$ with $t_{c}=\mu _{v}l_{c}/\sigma$. Note that (3.5) does not describe the true time evolution of the film thickness, since (3.3) is only valid in equilibrium; in fact, true dynamics depends on liquid viscosity and becomes infinitely slow as $\mu _l\to \infty$. However, the steady state of (3.5) is correct (indeed, putting $\partial \tilde {h}/\partial \tilde {t}=0$ turns it into the main equation of the theory of Sobac et al. Reference Sobac, Rednikov, Dorbolo and Colinet2014); the time derivative is retained for computational reasons and the steady-state profile is obtained as a result of time evolution in the long-time limit. It can be seen that (3.5) depends on a single dimensionless parameter

(3.6)\begin{equation} E=\frac{k_{v}\mu_{v}\Delta T}{\sigma L\rho_{v}l_{c}}, \end{equation}

which is the (dimensionless) evaporation number. In most cases, we find the condition $E \ll 1$ which represents slow evaporation, corroborating the consideration of a quasi-static process. Note that the exact expression for the curvature, not assuming the slow variation of film thickness,

(3.7)\begin{equation} \tilde{\kappa}=\frac{\dfrac{\partial^{2} \tilde{h}}{\partial \tilde{r}^{2}}}{\left [1+\left (\dfrac{\partial \tilde{h}}{\partial \tilde{r}} \right)^{2} \right]^{3/2}}+\dfrac{\dfrac{\partial \tilde{h}}{\partial \tilde{r}}}{\tilde{r} \left [1+\left (\dfrac{\partial \tilde{h}}{\partial \tilde{r}} \right)^{2} \right]^{1/2}} \end{equation}

is used in (3.5), which makes that equation valid in parts of the surface where the lubrication approximation is not valid (but the terms calculated using that approximation are negligible).

Equation (3.5) is fourth order in $\tilde {r}$ and therefore requires four boundary conditions: ${\rm d}\tilde {h}/{\rm d}\tilde {r}={\rm d}{\tilde {\kappa }}/{\rm d}{\tilde {r}}=0$ at the axis of symmetry $\tilde {r}=0$, and matching of ${\rm d}\tilde {h}/{\rm d}\tilde {r}$ and $\tilde {\kappa }$ at $\tilde {r}=\tilde {R}_{p}$ with the corresponding values obtained from the upper surface of the drop. After obtaining a steady-state solution, the continuity of $\tilde {h}(\tilde {r})$ itself simply amounts to translating vertically the upper surface of the drop, thus determining $\tilde {z}_{top}$ and revealing the complete equilibrium shape of the drop.

Whilst this formulation still requires a numerical solution, we refer to it as the ‘theoretical model’ or ‘theoretical solutions’ to distinguish it from the solutions to our computational model, which additionally solve for the internal flow via the Navier–Stokes equations.

To summarize, the theoretical model of Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) can be obtained from our computational model by neglecting the flow inside the drop and thus the results of the models are expected to coincide when this flow is negligible, i.e. when the liquid viscosity $\mu _l$ is large.

4. Leidenfrost drop shapes and identification of regimes

We are now in a position to compare our computational model with experimental analyses conducted on water drops and utilize it to establish the limits of applicability of the theoretical model. As described, the theoretical model neglects flow within the droplet and so to imitate this state with our computational model we align all parameters with those used in experiments (i.e. for water drops of a given size on a surface at a given temperature) with the exception of liquid viscosity which is set at a high value. This has the effect of ‘turning off’ the connection between the internal flow and the vapour film dynamics. The effect of internal flow can then be isolated and is investigated in § 5.

4.1. Shapes of Leidenfrost drops

Figure 2 shows the shapes of quasi-static Leidenfrost drops at equilibrium predicted by our computational model for four different sizes $\tilde {R}_{max}=0.26, 0.41, 2.69, 3.72$ at a fixed evaporation number $E=6.26\times 10^{-7}$ (all parameters are given in the caption) of a high-viscosity liquid ($\mu _{l}=0.3\,\mathrm {Pa}\,\mathrm {s}$). Superimposed are the theoretical solutions showing excellent agreement is obtained. Four distinct regimes are observed upon increasing the drop size:

1. (i) Dimpleless quasi-spherical drop – where the drop is approximately spherical as $\tilde {R}_{max}<1$ and the curvature does not change sign at the bottom of the drop, seen in figure 2(a) for $\tilde {R}_{max}=0.26$.

2. (ii) Dimpled quasi-spherical drop – where the drop is approximately spherical as $\tilde {R}_{max}<1$ but now a vapour pocket is formed under the drop, which is seen in figure 2(b) at $\tilde {R}_{max}=0.41$.

3. (iii) Puddle-like drop – where gravity flattens the drop as $\tilde {R}_{max}>1$, as seen in figure 2(c) for $\tilde {R}_{max}=2.69$.

4. (iv) Chimney instability – where stable quasi-static shapes no longer exist. In figure 2(d) we show a case $\tilde {R}_{max}=3.72$ close to where this instability occurs.

Figure 2. Comparison of the computational (black line) and theoretical (yellow dashed line) equilibrium shapes of Leidenfrost drops at $T_{b}=100\,^{\circ }\mathrm {C}$ above a hot rigid surface at $T_{w}=300\,^{\circ }\mathrm {C}$ for (a) $\tilde {R}_{max}=0.26$, (b) $\tilde {R}_{max}=0.41$, (c) $\tilde {R}_{max}=2.69$ and (d) $\tilde {R}_{max}=3.72$. At $T_{m}=200\,^{\circ }\mathrm {C}$, the input parameters (Biance et al. Reference Biance, Clanet and Quéré2003) are $\rho _{v}=0.5\,\mathrm {kg}\,\mathrm {m}^{-3}$, $\mu _{v}=1.63\times 10^{-5}\,\mathrm {Pa}\,\mathrm {s}$, $k_{v}=0.032\,\mathrm {W}\,\mathrm {m}^{-1}\,\mathrm {K}^{-1}$ and the latent heat of evaporation $L=2.26\times 10^{6}\,\mathrm {J}\,\mathrm {kg}^{-1}$ at $T_{b}=100\,^{\circ }\mathrm {C}$ and thus the calculated value of evaporation number $E=6.259\times 10^{-7}$. Here, simulations are carried out for a higher-than-water dynamic liquid viscosity of $\mu _{l}=0.3\,\mathrm {Pa}\,\mathrm {s}$.

It is noticeable from figure 2(b–d) that the region of the deformed drop adjacent to the steady vapour film takes a dimple shape (an approximately parabolic shape with curvature opposite to that of the undeformed sphere), with a maximum vapour film thickness $h_{centre}$ at the centre of the vapour layer $r=0$ and a minimum film thickness $h_{neck}$ defining the neck $r=R_{neck}$; see also figure 1. This dimpled regime is well known, and was first observed experimentally discovered for Leidenfrost droplets in Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012), having been theoretically predicted for a drop suspended by an upward air flow from a permeable solid in Duchemin et al. (Reference Duchemin, Lister and Lange2005) for curved substrates and Snoeijer et al. (Reference Snoeijer, Brunet and Eggers2009) for flat ones, and having been seen in a variety of other drop/gas-film interactions, such as drop impact (Chandra & Avedisian Reference Chandra and Avedisian1991; Thoroddsen, Etoh & Takehara Reference Thoroddsen, Etoh and Takehara2003). In contrast, the dimpleless regime in figure 2(a), also commented on in Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2021), is a new phenomenon.

4.2. Geometry of the vapour film

To enable a comparison with the experimental results in Biance et al. (Reference Biance, Clanet and Quéré2003) and Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012), we consider the geometry of the vapour film using the parameters from these experiments (with the exception of liquid viscosity). In particular, in figure 3 we present the dependence of the vapour film thickness at the neck $h_{neck}$, the difference in vapour film thickness $\Delta h=h_{centre}-h_{neck}$ and the film's neck radius $R_{neck}$ on the drop size $R_{max}$ varying over the range of $0.65\,\mathrm {mm} \le R_{max} \lesssim 10.0\,\mathrm {mm}$. This is done for two different wall temperatures $T_{w}=300\,^{\circ }\mathrm {C}\ \mathrm {and}\ 370\,^{\circ }\mathrm {C}$, corresponding to $E=6.259\times 10^{-7}\ \mathrm {and}\ 1.01\times 10^{-6}$, at high liquid viscosity $\mu _{l}=0.3\,\mathrm {Pa}\,\mathrm {s}$. In general, it is observed that the computed $h_{neck}$, $\Delta h$ and $R_{neck}$ all increase with $R_{max}$; $\Delta h$ and $R_{neck}$ are nearly independent of wall temperature $T_{w}$; and $h_{neck}$ somewhat increases with $T_{w}$.

Figure 3. Comparison of the computational model with experiments for the geometry of the vapour layer underneath a Leidenfrost drop. (a) The vapour film thickness at the neck $h_{neck}$ vs the drop size $R_{max}$, (b) the difference in the vapour film thickness $\Delta h=h_{centre}-h_{neck}$ as a function of $R_{max}$ and (c) the neck radius $R_{neck}$ vs $R_{max}$, see figures 1 and 2. The solid line represents the results of our computational model for $T_{w}=300\,^{\circ }\mathrm {C}\ \mathrm {and} 370\,^{\circ }\mathrm {C}$ and $\mu _{l}=0.3\,\mathrm {Pa}\,\mathrm {s}$, compared with the symbols corresponding to the experimental data reported by Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012) for $T_{w}=245\,^{\circ }\mathrm {C}, 320\,^{\circ }\mathrm {C}\ \mathrm {and}\ 370^{\circ }\mathrm {C}$ (for water drops), the dashed-dot lines are the solutions of the theoretical model; the vertical dotted line denotes the computed critical value of $R_{max,DL-D}$ below which ‘dimpleless’ (DL) regime with a nearly spherical drop shape can be seen, with the ‘dimpled’ (D) regime for higher $R_{max}$. In addition, the experimental data of Biance et al. (Reference Biance, Clanet and Quéré2003) for $h_{neck}$ (filled circle symbols) at $T_{w}=300\,^{\circ }\mathrm {C}$ (for water) are plotted in (a).

Notably, the computational model was able to identify differences with the solutions presented in Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) that, upon seeing our results (in Chakraborty et al. Reference Chakraborty, Chubynsky and Sprittles2020), the authors were able to correct in Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2021). In particular, for smaller values of $R_{max}$, where the dimpleless regime sits, the original solution of the theoretical model from Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) diverges from the corrected solution (Sobac et al. Reference Sobac, Rednikov, Dorbolo and Colinet2021).

In figure 3, we compare the result of our computational model with experiments and the theoretical model. Results for $\Delta h$ and $R_{neck}$ show very good agreement with the experimental data of Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012) when $R_{max}\gtrsim 1.0\,\mathrm {mm}$. For smaller values of $R_{max}$, we see discrepancies appearing between our predictions and experiments. This is worthy of further attention and could be caused by increased scatter from the indirect measurements as the film height shrinks (Burton et al. Reference Burton, Sharpe, van de Veen, Franco and Nagel2012) or due to the motion of the drop laterally (Bouillant et al. Reference Bouillant, Mouterde, Bourrianne, Lagarde, Clanet and Quéré2018). Qualitatively, it is noteworthy that experimental data for $\Delta h$ for 245 and $320\,^{\circ }\mathrm {C}$ rapidly decrease with decreasing $R_{max}$, indicating an approach to the dimple-less regime we have predicted. However, in figure 3(a) data for $h_{neck}$ from the experiment of Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012) go into the dimpleless regime, where the neck no longer exists in our computational model; the reason for this discrepancy remains unclear. Furthermore, our results in figure 3(a) show less good agreement with experimental data of Biance et al. (Reference Biance, Clanet and Quéré2003) for minimum film thickness $h_{neck}$, at $T_w = 300\,^{\circ }\mathrm {C}$, particularly for larger drop sizes. The reason is not clearly understood, but a probable explanation is that Biance et al. (Reference Biance, Clanet and Quéré2003) measure an effective vapour film thickness, and never actually discuss the height at the neck, and it is inferred here, as in Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014), that this value is closely related to the local measurement $h_{neck}$.

To further study the geometry of the vapour film and establish how/if the liquid viscosity influences it, in figure 4 we compute the dimensionless vapour film thickness at the centre $\tilde {h}_{centre}=h_{centre}/l_{c}$ and at the neck $\tilde {h}_{neck}=h_{neck}/l_{c}$ as a function of the dimensionless drop size $0.0016 \le \tilde {R}_{max} \lesssim 4.0$ for fixed $E=1.01\times 10^{-6}$ ($T_{w}=370\,^{\circ }\mathrm {C}$). In this case, which is used in many forthcoming calculations, the vapour density is $\rho _v=0.47\,{\rm kg}\,{\rm m}^{-3}$, its viscosity $\mu _v=1.70\times 10^{-5}\,{\rm Pa}\,{\rm s}$, its conductivity is $k_v=0.0347\,{\rm Wm}^{-1}\,{\rm K}^{-1}$ and the latent heat of vaporization is $L=2.26\times 10^{6}\,{\rm J}\,{\rm kg}^{-1}$. Calculations are done using both the dynamic viscosity of water $\mu_{l,water}=0.00034\,\mathrm {Pa}\,\mathrm {s}$ at $T_{b}=100\,^{\circ }\mathrm {C}$ and the high dynamic viscosity $\mu _{l,high}=0.3\,\mathrm {Pa}\,\mathrm {s}$ (used in previous results). The theoretical model's results are also plotted and have no dependence on internal flow and hence liquid viscosity. These are both compared with experimental data for Leidenfrost water drops (Burton et al. Reference Burton, Sharpe, van de Veen, Franco and Nagel2012) and are now considered by regime.

Figure 4. The vapour film thickness at the centre $\tilde {h}_{centre}$ and at the neck $\tilde {h}_{neck}$ vs the drop size $\tilde {R}_{max}$ for $E=1.01\times 10^{-6}$ ($T_{w}=370\,^{\circ }\mathrm {C}$). The red and blue solid/dashed lines correspond to the results of the computational model for $\mu _{l,water}=0.00034\, \mathrm {Pa}\,\mathrm {s}$ (water) and $\mu _{l,high}=0.3\, \mathrm {Pa}\,\mathrm {s}$ (high viscosity), respectively, compared with the theoretical solution. Regimes marked on the plot correspond to (i) dimpleless quasi-spherical drops, (ii) dimpled quasi-spherical drops, (iii) dimpled puddles and (iv) chimney instabilities. The experimental data for $\tilde {h}_{centre}$ are obtained from the data for $\tilde {h}_{centre}-\tilde {h}_{neck}$ and $\tilde {h}_{neck}$ given by Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012).

4.2.1. Dimpleless quasi-spherical drops: $\tilde {R}_{max} < \tilde {R}_{max,DL-D}=0.3$

In this regime, figure 4 shows that for both viscosities computational and theoretical curves for $\tilde {h}_{centre}$ merge with those for $\tilde {h}_{neck}$ indicating a dimpleless regime. Interestingly, in this regime, there is non-monotonic behaviour, with both heights initially decreasing with smaller $\tilde {R}_{max}$, reaching a minimum at $\tilde {R}_{max,hmin} \approx 0.14$, after which film heights increase, in qualitative agreement with the experimental discoveries made in Celestini et al. (Reference Celestini, Frisch and Pomeau2012). Computational and theoretical results are in excellent agreement above the minima ($0.14 \le \tilde {R}_{max} \le 0.3$) but discrepancies appear below it, where the accuracy of the lubrication approximation is reduced as the film height becomes comparable to the drop size; going beyond the lubrication approximation will be discussed in § 6. Notably, it can be observed (not shown here) that our results of $\tilde {h}_{min}$ and $\tilde {R}_{max,hmin}$ are only relatively weakly sensitive to the variation in the evaporation number $E$ or $T_{w}$.

4.2.2. Dimpled quasi-spherical drops: $\tilde {R}_{max,DL-D}<\tilde {R}_{max} <1$

Here, the dimpleless regime gives way to the well-known dimpled regime with figure 4 showing that film heights increase with $\tilde {R}_{max}$ and that for $\tilde {h}_{neck}$ the computational, for both viscosities, and theoretical models agree well with experimental trends from Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012).

4.2.3. Dimpled puddles: $1<\tilde {R}_{max} < \tilde {R}_{max,Ch} \approx 4.0$

Here, results for water viscosity from the computational model diverge from both experiments and the theoretical solutions. Specifically, for water drops the computational model predicts a complex non-monotonic behaviour and much smaller values of $h_{centre}$ are observed than for the high viscosity computations, the theoretical model and experiments, which all align. This will be probed further in § 5.

4.2.4. Chimney instability: $\tilde {R}_{max} \to \tilde {R}_{max,Ch}$

In figure 4, we show two branches of the theoretical solution, with a stable lower branch meeting an upper unstable one at a ‘turning point’ or ‘fold’ (shown by an arrow) indicative of a saddle-node bifurcation point. (Strictly speaking, since it is the drop volume, or, equivalently, its initial radius $\tilde {R}_0$, that is fixed and not $\tilde {R}_{max}$, the bifurcation occurs at the point along the solution curve where $\tilde {R}_0$, rather than $\tilde {R}_{max}$, reaches its maximum, which is slightly below the turning point on the plot.) This point defines the onset of the chimney instability whose dynamics will be considered further by our computational model in § 7. The advantage of the theoretical model of Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) is that it can also reveal the unstable branch, whereas obtaining this in the computational model is difficult and sometimes impossible if the unstable solution corresponds to a self-intersecting curve. A similar turning point at $\tilde {R}_{max,Ch}$, of course, exists on the $\tilde {h}_{neck}$ curve, as well as the low-viscosity curves, for which the instability threshold is nearly the same; for these curves the unstable branches are not shown, but the existence of the turning point is also not apparent from the behaviour of the stable branches, as the slopes of the curves rise and approach infinity in a very narrow region near the threshold.

5. Internal flow

Having observed an unexpected dependence of $\tilde {h}_{centre}$ on liquid viscosity, in figure 5 a more detailed study is shown for a range of liquid viscosities. This reveals that $\tilde {h}_{centre}$ is independent of liquid viscosity for $\tilde {R}_{max}<1$ and above this the curves vary smoothly between the two viscosities previously presented. Notably, figure 5(b,c) shows that the computed $\tilde {h}_{neck}$ and $\tilde {R}_{neck}$ depend only weakly on the liquid viscosity in contrast to $\tilde {h}_{centre}$.

Figure 5. (a) The film height at the centre of the vapour pocket $\tilde {h}_{centre}$, (b) the neck height of the vapour film $\tilde {h}_{neck}$ and (c) the neck radius $\tilde {R}_{neck}$ as a function of the drop size $\tilde {R}_{max}$ obtained from COMSOL simulations for $E=1.01\times 10^{-6}$ at the wall temperature $T_{w}=370\,^{\circ }{\rm C}$ and different liquid viscosities. The square symbols for $\tilde {h}_{neck}$ and $\tilde {R}_{neck}$ represent the experimental data (Burton et al. Reference Burton, Sharpe, van de Veen, Franco and Nagel2012).

To study further the influence of liquid viscosity and internal flow, in figure 6 we show the results of our computational model for drops with an initial radius $R_{0}=3.75\,\mathrm {mm}$ for three different viscosities ($\mu _{l}=0.3\,\mathrm {Pa}\,\mathrm {s},\ 0.001\,\mathrm {Pa}\,\mathrm {s}\ \mathrm {and}\ 0.00034\,\mathrm {Pa}\,\mathrm {s}$) with at fixed $E=1.01\times 10^{-6}$. For a highly viscous drop as shown in figure 6(a), there is very little internal flow, the shape of the drop becomes a puddle, a dimple is formed and the computed drop size becomes $R_{max}=4.63\, \mathrm {mm}$ ($\tilde {R}_{max}=1.85$). This result and the corresponding geometry of the vapour pocket below the drop (see figure 6d) are in excellent agreement with experiment (Burton et al. Reference Burton, Sharpe, van de Veen, Franco and Nagel2012) and theory (Sobac et al. Reference Sobac, Rednikov, Dorbolo and Colinet2014); see also figure 3.

Figure 6. The steady equilibrium shape of a drop with initial radius $R_{0}=3.75\,\mathrm {mm}$ and the velocity field inside the drop including the velocity vectors of flow for (a) $\mu _{l}=0.3\,\mathrm {Pa}\,\mathrm {s}$, (b) $\mu _{l}=0.001\,\mathrm {Pa}\,\mathrm {s}$ and (c) $\mu _{l}=0.00034\,\mathrm {Pa}\,\mathrm {s}$ (water at $T_{b}=100\,^{\circ }\mathrm {C}$) and the fixed evaporation number $E=1.01\times 10^{-6}$ placed on a very hot surface at $T_{w}=370\,^{\circ }\mathrm {C}$. The colour denotes the corresponding velocity magnitude inside the drop. (d) The air film profiles for the cases presented in (a–c).

In contrast, the shapes of the low-viscosity drop with $R_{max}=3.84\, \mathrm {mm}$ ($\tilde {R}_{max}=1.54$) shown in figure 6(b) and the water drop with $R_{max}=3.46\, \mathrm {mm}$ ($\tilde {R}_{max}=1.39$) shown in figure 6(c) display almost no dimple; instead a ‘wimple’ like shape of the vapour film is observed in figure 6(d). This shape has also been observed in a previous investigation (Chan, Klaseboer & Manica Reference Chan, Klaseboer and Manica2011) for the case of a drop suspended by a thin air film above a solid surface and similar theoretical results for the case of a large drop ($\tilde {R}_{max} \gtrsim 1$) are reported by Maquet et al. (Reference Maquet, Sobac, Darbois-Texier, Duchesne, Brandenbourger, Rednikov, Colinet and Dorbolo2016), who studied the levitation of ethanol drops on a vapour layer above the surface of a heated liquid pool. Moreover, it can be seen from figure 6(b,c) that the drop attains a distinctly different prolate shape.

The observed situation is rather unusual, the more accurate computational model suddenly starts giving worse agreement with experiments than the less accurate theoretical model, which neglects internal flow. The complexity of the process is such that there are many potential reasons for this situation, but the weakest assumptions in our model are (i) that the liquid drop is at a constant (boiling) temperature and (ii) that the system is axisymmetric; these possibilities are discussed in § 8.

6. Beyond lubrication theory: $\tilde {R}_{max}<\tilde {R}_{max,hmin} \approx 0.14$

Consider now the dimpleless small-drop regime, in which film heights start to increase with decreasing drop sizes (lower $\tilde {R}_{max}$) according to both the computational and theoretical models, which are underpinned by lubrication theory, see figure 7. To make further analytic progress, one can assume the drop takes a rigid-sphere shape, and derive an asymptotic expression $\tilde {h}_{centre} \sim \tilde {R}^{-1/2}_{max}$ in the limit of large $R_{max}$, see Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2021), which highlights the increase of $\tilde {h}_{centre}$ with decreasing $\tilde {R}_{max}$. This plot also shows at what point deformability of the sphere becomes important, as the rigid-sphere theoretical result starts to deviate from that making no such assumption, which occurs near where there is a minimum in $\tilde {h}_{centre}$ with a subsequent increase for bigger drop sizes (i.e. larger $\tilde {R}_{max}$). This increase is thus clearly due to deformability of the drop surface.

Figure 7. Variation of $\tilde {h}_{centre}$ with the drop radius $\tilde {R}_{max}$ for $E=1.01\times 10^{-6}$ ($T_{w}=370\,^{\circ }\mathrm {C}$ for water) for small- and medium-sized drops. Different approaches are used: the theoretical (blue) and computational (red) approaches used in the rest of the paper, a similar theoretical approach assuming that the drop is a rigid sphere (Sobac et al. Reference Sobac, Rednikov, Dorbolo and Colinet2021) (orange), the same approach with a different patching point (magenta) and, finally, an approach going beyond lubrication, as described in the text (turquoise). The short dashed lines denote a scaling exponent $-1/2$ obtained (in both limits) from the scaling laws derived in Celestini et al. (Reference Celestini, Frisch and Pomeau2012); the prefactor for large $R_{max}$ is known (Sobac et al. Reference Sobac, Rednikov, Dorbolo and Colinet2021), while for small $R_{max}$ it is estimated based on the computational data.

However, figure 7 also shows that as $\tilde {R}_{max}$ decreases and $\tilde {h}_{centre}$ increases, eventually the two become comparable so that the lubrication approximation becomes invalid. In this regime, there is a discrepancy between the computational and theoretical solutions, but also a dependence of the theoretical solution on the exact position of the matching point $R_{p}$ (see figure 1).

In fact, in this non-lubrication regime, scaling arguments in Celestini et al. (Reference Celestini, Frisch and Pomeau2012) for $\tilde {h}_{centre}\gg \tilde {R}_{max}$ suggest that the characteristic film thickness also follows the power law dependence $\tilde {h}_{centre} \sim \tilde {R}^{-1/2}_{max}$, with the pre-factor now unknown. In this regime, one can likewise make the rigid-sphere assumption for the drop, but, on the other hand, the vapour flow around it needs to be considered in full. The consideration is simplified significantly by making additional assumptions of various degrees of realism, including Stokes flow of the vapour and its uniformity in space (even though it is, in fact, a vapour–air mixture of a varying composition). Under such assumptions, conveniently, $\tilde {h}(\tilde {R}_{max})$ still depends on just a single parameter $E$. For details, see the Appendix B. The result is plotted in figure 7 and approaches the asymptotics in both limits, as expected.

For small drops, we see that lubrication theory loses its accuracy as the vapour film loses its high aspect ratio. To capture these regimes reliably one must solve the full Navier–Stokes–Fourier equations in the ambient phase and account for the interplay of air and water vapour. In fact, this is not so complex as the multiscale nature of the problem is lost in this regime and the equations for the physical processes are well known.

7. Dynamic chimney instability

When the drop size is larger than the threshold for the chimney instability ($\tilde {R}_{max,Ch}$), the quasi-static state is lost and dynamic evolution of the interface commences, with a vapour pocket under the drop growing towards the upper surface. The values we find for this instability of $\tilde {R}_{max,Ch} \approx 4.0$ for both models and all viscosities are close to the predictions in the literature of $\tilde {R}_{max,Ch}=3.84$ in Biance et al. (Reference Biance, Clanet and Quéré2003) and $3.95$ in Snoeijer et al. (Reference Snoeijer, Brunet and Eggers2009) and Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012), where this instability is attributed to the Rayleigh–Taylor mechanism, with shapes similar to those observed for inviscid drop computations in Bouwhuis et al. (Reference Bouwhuis, Winkels, Peters, Brunet, van der Meer and Snoeijer2013). In particular, careful calculations for the theoretical model with $E=1.014\times 10^{-6}$ corresponding to water with $T_w=370^{\circ }$ give the maximum value of $\tilde {R}_0$ on the solution curve $\tilde {R}_{0,Ch}\approx 2.558$, corresponding to $\tilde {R}_{max,Ch}\approx 3.944$, slightly below the maximum value of $\tilde {R}_{max}$, 3.973, and in the simulations of the computational model for the same temperature and $\mu _l=0.3\, \mathrm {Pa}\,\mathrm {s}$ very similar values of 2.569 and 3.970, respectively, are obtained.

Physically, the vapour pocket eventually bursts, i.e. penetrates the upper surface. Reassuringly, in our computational model, the chimney is formed in around a second, in agreement with observations in Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012). As previously described, our computational model can capture the dynamics and its prediction and simulation of this event highlight these capabilities, although the code is terminated prior to film rupture/pocket bursting as this topological change requires further work. Nevertheless we are able to capture the development of the instability in figure 8, which illustrates a quick formation of a dimpled shape followed by a slow growth of the instability. For such large drops the negative curvature of the bottom surface cannot be sufficiently large to balance the difference in hydrostatic pressure between the top and bottom of the drop, thus, the thin liquid ‘film’ at the top gradually drains until it breaks up.

Figure 8. Simulation of the chimney instability for a drop with an initial drop radius $R_0=6.45$ mm ($\tilde {R}_0=2.58$), which is above the threshold for the chimney instability (calculated using the theoretical model to be at $\approx \tilde {R}_0=2.56$). Here, $E=1.01\times 10^{-6}$ ($T_w=370\,^{\circ }\mathrm {C}$ for water) and viscosity $\mu =0.3\,\mathrm {mPa}\,\mathrm {s}$. (a) The evolution of the shape of the drop from an initially spherical drop, through to a puddle and then into a chimney instability and (b) the flow field with the colormap showing the velocity magnitude and the arrows indicating the direction of flow.

9. Concluding remarks

A computational model has been developed that allows us to go beyond the theory of Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) by taking into account internal flow and drop dynamics. Comparisons with Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2014) identified discrepancies that were subsequently corrected in Sobac et al. (Reference Sobac, Rednikov, Dorbolo and Colinet2021) and show that for a large range of parameters this theory does an excellent job of predicting the characteristics of quasi-static Leidenfrost drops. In particular, both the model and theory allowed us to discover that small drops exhibit a transition to a dimpleless drop which, as far as we are aware, has not yet been recovered experimentally, being on the edge of what is measured in Burton et al. (Reference Burton, Sharpe, van de Veen, Franco and Nagel2012), see figure 3(a), and thus should be a stimulus for further analyses. The limits of applicability of the model have been discussed in § 8, and directions for future work proposed.

Finally, our model demonstrated, via the chimney instability, its capability to capture dynamic processes that will be the focus of future work. In particular, studying drop impact Leidenfrost is of most interest, as it represents a highly practically relevant case and many fascinating high-accuracy experimental results have recently been obtained there (Tran et al. Reference Tran, Hendrik, Prosperetti and Lohse2012; Shirota et al. Reference Shirota, van Limbeek, Sun, Prosperetti and Lohse2016). Importantly, as for drop impact in isothermal conditions (Chubynsky et al. Reference Chubynsky, Belousov, Lockerby and Sprittles2020), the developed framework exhibits accuracy beyond what may be expected, as for most drop sizes the lubrication pressure only manifests itself when the vapour film is long and thin, i.e. when the lubrication approximation is indeed valid.