Computational modelling of Leidenfrost drops

Abstract The Leidenfrost effect, where a drop levitates on a vapour film above a hot solid, is simulated using an efficient computational model that captures the internal flow within the droplet, models the vapour flow in a lubrication framework and is capable of resolving the dynamics of the process. The initial focus is on quasi-static droplets and the associated geometry of the vapour film formed beneath the drop, where we are able to compare with experimental analyses and assess the range of validity of the theoretical model developed in Sobac et al. (Phys. Rev. E, vol. 103, 2021, 039901). The computational model also allows us to explore parameter space, varying both the drop size and viscosity of the liquid, with computational results in excellent agreement with the theoretical model for high-viscosity liquids. Interestingly, for large water drops, discrepancies between the computational model and experiments occur, and possible reasons for this observation are provided. Our predictions reveal features including a regime with a dimpleless bottom surface of the drop and a minimum in the vapour layer thickness as a function of the drop size. Finally, the capability to simulate dynamics is revealed by computations that predict and track the vapour ‘chimney’ instability for large drops.

across a very hot pan, but is also crucial to numerous drop-based technologies, where the Leidenfrost effect can prevent efficient heat transfer, such as the spray cooling of high performance electronic devices (Kim 2007) and spray combustion (Moreira, Moita & Panã 2010). Moreover, recent attention has been given to fascinating new topics, such as the self-propulsion of drops on ratchet surfaces (Linke et al. 2006) and hydrodynamic drag reduction (Vakarelski et al. 2012). Fundamental studies of Leidenfrost drops, and numerous applications, are reviewed in great detail by Quéré (2013) and Ling & Mudwar (2017).
At the most basic level, one would like to understand how and when the vapour film is able to retain a steady cushion for the drop and this naturally leads to a study of the geometry of the vapour film, which is complex experimentally due to its thin and almost hidden nature. In Biance et al. (2003) and Burton et al. (2012) these challenges were overcome in order to experimentally record the equilibrium shapes of Leidenfrost water drops and measure the geometry of the vapour film for different drop sizes. It was demonstrated that the film thickness increases with increasing drop radii, being around 30-100 μm for drop radii of 1.6-7 mm. Using lubrication theory, and assuming an approximately uniform film, Biance et al. (2003) found good agreement with their experimental measurements, with two regimes identified based on whether the size of the drop, defined by its maximum radius R max , is larger or smaller than the capillary length l c = (σ/ρ l g) 1/2 (dimensionlesslyR max = R max /l c , with˜henceforth denoting dimensionless quantities), where σ is the surface tension of the liquid-vapour interface, g is gravitational acceleration and ρ l is the density of the liquid. Notably, for water, Burton et al. (2012) showed that the vapour layer cannot be considered uniform and observed that the minimum height of this layer, at the 'neck', near its edge, varies from 5-100 μm for drop sizes of radii 0.5-10 mm, thus motivating the development of more complex models for the film.
Variations in the profile of the vapour layer are accompanied by different drop shapes: for small drops (R max < 1), surface tension dominates gravity so that the drop shape becomes quasi-spherical, with a slightly flattened base near the solid surface. Interestingly, for very small Leidenfrost drops (radii 1-30 μm), gap thicknesses actually increase with drop radius, as predicted and confirmed experimentally in Celestini, Frisch & Pomeau (2012), eventually leading to 'take off', but this regime is yet to be considered in detail computationally/theoretically.
For large drops (R max > 1), gravity dominates so that a flat puddle shape is formed with a concave pocket-like geometry of the vapour layer (Burton et al. 2012). For a sufficiently large puddle drop, a vapour chimney forms underneath the puddle and eventually bursts at the upper surface of this drop. The onset of a chimney instabilityR max ≈ 3.84 is analytically estimated by Biance et al. (2003) by examining the Rayleigh-Taylor instability (Taylor 1950) of the liquid-vapour interface and then verified experimentally. Notably, Snoeijer, Brunet & Eggers (2009) also found a critical radiusR max ≈ 4.0 using a model that considered the static shapes of drops levitated by a lubricating film of air injected with constant velocity through the substrate, under isothermal conditions. For large drops, spontaneous symmetry-breaking oscillations can be observed that create the appearance of star-shaped drops (Ma, Liétor-Santos & Burton 2017;Brunet & Snoeijer 2011), but these are beyond the scope of this article.
For the study of quasi-static Leidenfrost drops, Sobac et al. (2014) proposed a theoretical model that balances the evaporation-driven viscous lubrication, hydrostatic and capillary pressures in the vapour film and matches this to a (Young-Laplace type) capillary and hydrostatic balance for the upper surface of the drop, under the assumption of axisymmetry. There are two main assumptions of the model: (i) the process is assumed to be quasi-static, with the evaporation time of the drop much longer than the viscous and thermal relaxation times in the vapour, so that the change of the drop radius can be neglected, and (ii) the liquid motion inside the drop is neglected. The theoretical predictions show very good agreement with experiments, for both the shapes of Leidenfrost water drops and the geometry of the vapour film. Despite the agreement, we will revisit the theory of Sobac et al. (2014) in § 2, as we have seen that the original numerical solutions required reassessment (Chakraborty, Chubynsky & Sprittles 2020); a fact which motivated our present study and led to an Erratum being published in Sobac et al. (2021) after the authors discovered a typo in their code.
Recently, there has been interest in the dynamics of drops impacting on hot surfaces to predict the transition between the boiling and Leidenfrost regime; a question of key importance for cooling technologies (see Ling & Mudwar 2017). Similarly to the static case, above the dynamic Leidenfrost temperature the drop does not touch the substrate (Tran et al. 2012;Shirota et al. 2016). Here, we will focus on developing a reliable computational model for the case of quasi-static Leidenfrost drops, with some dynamics captured for the chimney instability, and extend this in future work to study impact events.
In the aforementioned investigations, Leidenfrost drops are typically of millimetre size while the underlying vapour films are two orders of magnitude smaller, revealing a multiscale problem that makes direct computational approaches challenging. Numerical studies of the dynamics of droplets impacting on a hot surface above T L have been reported using the level set/arbitrary Lagrangian Eulerian methods (Ge & Fan 2005) and direct simulations of level set/ghost fluid methods (Villegas et al. 2017, and references therein). By using the interface capturing techniques and solving the full Navier-Stokes-Fourier equations in the whole computational domain, the drop's shape evolution and the thickness of thin vapour layer during the impact process were examined. These approaches require, on one hand, high density mesh in both the liquid and surrounding medium, especially in the thin vapour layer close to the hot surface, and on the other hand, very small time steps to produce stable and accurate numerical solutions, rendering such approaches computationally expensive and extremely challenging to produce accurate results, particularly as one approaches the transition between contact and bouncing .
Despite a wealth of experimental data, the computational modelling of quasi-static Leidenfrost drops over a range of parameters remains lacking from the literature. Furthermore, little research has been concerned with the interior flow and hence the effect of liquid viscosity on the process. In the present work, a robust hybrid/multiscale computational model is proposed based on coupling the lubrication approximation for the vapour flow to the Navier-Stokes equations for the flow within the drop. Numerical simulations are conducted using finite-element software and extend the work of Chubynsky et al. (2020), who studied the isothermal impact of a droplet on a solid surface to capture the suspension of a millimetre-sized drop by a nanoscale air film. The notable advantage of our approach is that only the liquid domain requires meshing with finite elements, with the vapour manifesting itself through the drop's boundary conditions, resulting in an efficient numerical method.
The paper is organized as follows. In § § 2 and 3, we present two modelling approaches designed to provide insight into the Leidenfrost phenomenon: the hybrid/multiscale computational model and the theoretical framework of Sobac et al. (2014). In § 4, we present and discuss the computational results for the various shape regimes of quasi-static Leidenfrost drops over a wide range of parameters that enable comparisons with previous 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, θ) 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é (2013) and ( experimental analyses (Biance et al. 2003;Burton et al. 2012). Interestingly, in § 5 we discover an unexpected effect of liquid viscosity and divergence from experiments, which motivates further study. Next, in § 6, the limits of applicability of the lubrication approximation are established and a method to go beyond this is considered. Finally, in § 7, we show that our model is able to capture the dynamic chimney instability, before making concluding remarks in § 9.

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 ( 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. We consider an axisymmetric problem, an assumption which is discussed later, described in the cylindrical coordinate system (r, z, θ), 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 θ. 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. (2014).
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. (2014). 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.

Liquid flow
The flow of liquid inside the drop is governed by the isothermal incompressible Navier-Stokes equations, with incompressibility 1 r ∂(ru l ) ∂r + ∂w l ∂z = 0, (2.1) and the momentum equations where henceforth a subscript l denotes the liquid phase and v will denote the vapour. Here, ρ will represent densities; μ dynamic viscosities; u and w are velocity components in the r and z directions; p is the pressure; and t is the time.

Vapour flow
In the vapour film we develop a lubrication model, with height h(r, t) considered slowly varying (∂h/∂r 1) and small compared with the drop size R max ( = h/R max 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 Re v 1 (Aursand, Davis & Ytrehus 2018), where the Reynolds number of the thin vapour layer is defined as (Biance et al. 2003;Celestini et al. 2012;Aursand et al. 2018). This establishes the use of lubrication theory in the vapour film (Biance et al. 2003;Sobac et al. 2014;Aursand et al. 2018), with incompressibility unchanged (2.4) and the momentum equations simplifying to The velocity component tangent to the interface (z = h(r, t)) is assumed continuous across it, so that u v,Γ = u l,Γ ≡ u Γ , where Γ indicates a surface property of the liquid-vapour interface and subscripts v, Γ and l, Γ 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 Γ at the drop surface, the local radial velocity profile in the lubricating film is obtained as . (2.6) 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 , and k v and c p,v are thermal conductivity and specific heat capacity (see Aursand et al. 2018) . Hence, in the lubrication approximation the energy conservation equation simplifies to which indicates that the temperature across the vapour film varies linearly between the solid surface temperature T v | z=0 = T w and the liquid-vapour interface temperature Let the rate of change of the vapour film height in the frame of reference moving horizontally with speed u Γ be w Γ . Then in the laboratory frame (2.8) Without evaporation, the usual kinematic boundary conditions then give w Γ = w l,Γ = w v,Γ , where w l,Γ and w v,Γ 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) 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 η = ρ v /ρ l ∼ 10 −3 , (2.9) implies that w v,Γ − w Γ is very large compared with w l,Γ − w Γ ; thus, we can assume w Γ = w l,Γ . Furthermore, the low density ratio η 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 Hence, (2.9) with (2.8) at z = h can be rewritten as (2.11) 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 | z=0 = 0 and w v | z=h = w v,Γ (obtained from (2.11)), and applying the Leibniz integral rule, we get 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 ∂p v /∂r ≈ 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 The corresponding boundary conditions are given by (2.14) 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 τ v (r, t) = μ v ( ∂u v /∂z|) z=h . This leads to with, as expected, terms from both Poiseuille and Couette components.

Boundary conditions
In order to solve the Navier-Stokes equations (2.1)-(2.3), we need the free surface boundary conditions: normal and shear stresses are enforced on the entire drop surface. For this purpose, the drop surface is divided into two parts, as considered also in the theory of Sobac et al. (2014), separated at r = r o . The top of the drop (part 1) is modelled conventionally with the ambient pressure equal to the atmospheric pressure (p o = 0), the normal stress is equal to the Laplace pressure, and the shear stress is zero. The bottom of the drop (part 2), adjacent to the vapour layer, has (i) the normal stress set to the sum of Laplace pressure and the lubrication pressure p v (r, t), where p v (r, t) is obtained by solving (2.13) and (ii) the shear stress τ v (r, t) given by (2.15). Note that p v (r, t) and τ v (r, t) depend on the radial velocity of the liquid (u Γ = u l,Γ ) at the free surface, and ∂h/∂t in (2.8) depends on the vertical velocity of the liquid (w Γ = w l,Γ ) at the interface, which are obtained simultaneously by solving the Navier-Stokes equations for the drop. We choose r o to be at the maximum radial extent of the drop surface where the vapour pressure is nearly equal to the atmospheric pressure. The computed results are insensitive to the exact value of r o providing at that point the value of h is much larger than that in the vapour layer at its thinnest point; see the supplemental material of Chubynsky et al. (2020) for more details. Note also that the conventional kinematic boundary condition for the evolution of the free surface of the drop in terms of the liquid velocity at the surface applies everywhere on the surface, since in parts of the surface where evaporation is not negligible, u l,Γ = u Γ and w l,Γ = w Γ is still assumed.

Computational implementation
Our hybrid/multiscale computational model is solved in Comsol Multiphysics (COMSOL Ltd., Cambridge, UK; version 5.4) using finite elements, as described in Appendix A.

Theoretical model
This section reviews the theoretical framework proposed by Sobac et al. (2014) 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 dz/dr = 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, where z 1 = z − z top is the vertical coordinate measured from the apex of the drop (z = z top ) and κ 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), Givenκ apex , (3.2) can be solved to obtain the shape of the upper surface, in the form of the dependencesz 1 (s) andr(s), wheres is the arc length measured from the apex (Duchemin, Lister & Lange 2005). Rather than specifying the size of the dropR max , it is more convenient to chooseκ apex instead andR max is then found from the solution. Note that at this pointz 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 We next use (2.12), where we again neglect the liquid flow and therefore u Γ , which gives Using (3.3) to eliminate p v and switching to dimensionless variables, we get 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 μ l → ∞. However, the steady state of (3.5) is correct (indeed, putting ∂h/∂t = 0 turns it into the main equation of the theory of Sobac et al. 2014); 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 which is the (dimensionless) evaporation number. In most cases, we find the condition E 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,κ 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 inr and therefore requires four boundary conditions: dh/dr = dκ/dr = 0 at the axis of symmetryr = 0, and matching of dh/dr andκ atr =R p with the corresponding values obtained from the upper surface of the drop. After obtaining a steady-state solution, the continuity ofh(r) itself simply amounts to translating vertically the upper surface of the drop, thus determiningz 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. (2014) 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 μ l is large.

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. 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. (2012), having been theoretically predicted for a drop suspended by an upward air flow from a permeable solid in Duchemin et al. (2005) for curved substrates and Snoeijer et al. (2009) for flat ones, and having been seen in a variety of other drop/gas-film interactions, such as drop impact (Chandra & Avedisian 1991;Thoroddsen, Etoh & Takehara 2003). In contrast, the dimpleless regime in figure 2(a), also commented on in Sobac et al. (2021), is a new phenomenon.

Geometry of the vapour film
To enable a comparison with the experimental results in Biance et al. (2003) and Burton et al. (2012), 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 h = h centre − h neck and the film's neck radius R neck on the drop size R max Here, simulations are carried out for a higher-than-water dynamic liquid viscosity of μ l = 0.3 Pa s.  varying over the range of 0.65 mm ≤ R max 10.0 mm. This is done for two different wall temperatures T w = 300 • C and 370 • C, corresponding to E = 6.259 × 10 −7 and 1.01 × 10 −6 , at high liquid viscosity μ l = 0.  In figure 3, we compare the result of our computational model with experiments and the theoretical model. Results for h and R neck show very good agreement with the experimental data of Burton et al. (2012) when R max 1.0 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. 2012) or due to the motion of the drop laterally (Bouillant et al. 2018). Qualitatively, it is noteworthy that experimental data for h for 245 and 320 • 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. (2012) 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. (2003) for minimum film thickness h neck , at T w = 300 • C, particularly for larger drop sizes. The reason is not clearly understood, but a probable explanation is that Biance et al. (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. (2014), 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 centrẽ h centre = h centre /l c and at the neckh neck = h neck /l c as a function of the dimensionless drop size 0.0016 ≤R max 4.0 for fixed E = 1.01 × 10 −6 (T w = 370 • C). In this case, which is used in many forthcoming calculations, the vapour density is ρ v = 0.47 kg m −3 , its viscosity μ v = 1.70 × 10 −5 Pa s, its conductivity is k v = 0.0347 Wm −1 K −1 and the latent heat of vaporization is L = 2.26 × 10 6 J kg −1 . Calculations are done using both the dynamic viscosity of water μ l,water = 0.00034 Pa s at T b = 100 • C and the high dynamic viscosity μ l,high = 0.3 Pa 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. 2012) and are now considered by regime. 4.2.1. Dimpleless quasi-spherical drops:R max <R max,DL−D = 0.3 In this regime, figure 4 shows that for both viscosities computational and theoretical curves forh centre merge with those forh neck indicating a dimpleless regime. Interestingly, in this regime, there is non-monotonic behaviour, with both heights initially decreasing with smallerR max , reaching a minimum atR max,hmin ≈ 0.14, after which film heights increase, in qualitative agreement with the experimental discoveries made in Celestini et al. (2012). Computational and theoretical results are in excellent agreement above the minima (0.14 ≤R max ≤ 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 ofh min andR max,hmin are only relatively weakly sensitive to the variation in the evaporation number E or T w .

Dimpled quasi-spherical drops:R max,DL−D <R max < 1
Here, the dimpleless regime gives way to the well-known dimpled regime with figure 4 showing that film heights increase withR max and that forh neck the computational, for both viscosities, and theoretical models agree well with experimental trends from Burton et al. (2012).

4.2.3.
Dimpled puddles: 1 <R max <R max,Ch ≈ 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.

Chimney instability:R max →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 radiusR 0 , that is fixed and notR max , the bifurcation occurs at the point along the solution curve whereR 0 , rather thanR 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. (2014) 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 atR max,Ch , of course, exists on theh 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.

Internal flow
Having observed an unexpected dependence ofh centre on liquid viscosity, in figure 5 a more detailed study is shown for a range of liquid viscosities. This reveals thath centre is independent of liquid viscosity forR max < 1 and above this the curves vary smoothly between the two viscosities previously presented. Notably, figure 5(b,c) shows that the computedh neck andR neck depend only weakly on the liquid viscosity in contrast toh centre .
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 mm for three different viscosities (μ l = 0.3 Pa s, 0.001 Pa s and 0.00034 Pa s) with at fixed E = 1.01 × 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 mm (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. 2012) and theory (Sobac et al. 2014); see also figure 3. In contrast, the shapes of the low-viscosity drop with R max = 3.84 mm (R max = 1.54) shown in figure 6(b) and the water drop with R max = 3.46 mm (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 2011) 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 (R max 1) are reported by Maquet et al. (2016), 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 . Variation ofh centre with the drop radiusR max for E = 1.01 × 10 −6 (T w = 370 • 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. 2021) (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. (2012); the prefactor for large R max is known (Sobac et al. 2021), while for small R max it is estimated based on the computational data.
(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:R max <R max,hmin ≈ 0.14 Consider now the dimpleless small-drop regime, in which film heights start to increase with decreasing drop sizes (lowerR 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 expressionh centre ∼R −1/2 max in the limit of large R max , see Sobac et al. (2021), which highlights the increase ofh centre with decreasingR 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 inh centre with a subsequent increase for bigger drop sizes (i.e. largerR max ). This increase is thus clearly due to deformability of the drop surface.
However, figure 7 also shows that asR max decreases andh 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. (2012) for h centre R max suggest that the characteristic film thickness also follows the power law dependenceh centre ∼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,h(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.

Dynamic chimney instability
When the drop size is larger than the threshold for the chimney instability (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 ofR max,Ch ≈ 4.0 for both models and all viscosities are close to the predictions in the literature ofR max,Ch = 3.84 in Biance et al. (2003) and 3.95 in Snoeijer et al. (2009) and Burton et al. (2012), where this instability is attributed to the Rayleigh-Taylor mechanism, with shapes similar to those observed for inviscid drop computations in Bouwhuis et al. (2013). In particular, careful calculations for the theoretical model with E = 1.014 × 10 −6 corresponding to water with T w = 370 • give the maximum value ofR 0 on the solution curveR 0,Ch ≈ 2.558, corresponding toR max,Ch ≈ 3.944, slightly below the maximum value ofR max , 3.973, and in the simulations of the computational model for the same temperature and μ l = 0.3 Pa 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. (2012). 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.

Discussion
One of the most surprising findings of our work is the poor agreement of the computational model with experimental analyses for large water drops, where the model predicts a strong internal flow that creates a prolate drop profile which is not observed experimentally. As a starting point, we consider whether the magnitude of the velocities predicted by the computational model (≈0.4 m s −1 ) is commensurate with simple estimates and, after, we discuss the possible reasons for the discrepancy between computations and experiments. 8.1. Estimates for internal flow Having in mind the situation in figure 6(c), with the corresponding vapour film profile in figure 6(d), we consider the approximation where the film has a constant thickness h and radius R f . Notably, the assumption of a constant film height cannot be made for the high-viscosity shapes as in figure 6(a), where Laplace pressure gradients compensate vapour pressure gradients. Even though the flow in the drop will drive Couette flow in the film, we assume that Poiseuille flow still dominates so the Couette term (containing u Γ ) in (2.12) can be neglected in that equation (this assumption will be checked later). Since the stationary case is considered, the time derivative term vanishes. Then the solution for the pressure is and, by integrating this over the circle of radius R f , the total pressure force is obtained This force supports the weight of the drop, ρ l gV, where V is its volume, which gives Using E = 10 −6 , liquid parameters for water at 100 • C, and, according to figure 6(c), R f = 3 mm and the volume corresponding to that of a sphere of radius 3.5 mm, this gives a film height of 0.076 mm, consistent with the wimple profile in figure 6(d).
There are two ways in which the vapour film can drive flow in the drop: first, due to shear stress on the liquid-vapour interface from the Poiseuille flow in the film and, second, due to pressure gradients that exist in the film, which are nearly continuous across the vapour-liquid interface (as it is almost flat).

Shear-stress-induced flow
The shear stress on the liquid-vapour interface from the Poiseuille flow in the film is given by This stress is continuous across the interface and drives flow in the drop. The horizontal velocity of this flow drops rapidly in a boundary layer above the interface. Assuming that the thickness of this layer is much smaller than R f and ignoring the pressure gradients, the equation for the horizontal component of the flow velocity is Here, z is measured from the liquid-vapour interface. The solution of this equation vanishing at z = ∞ and with the boundary condition (8.4) at z = 0, to ensure continuity of shear stress across the interface, is where U = (3T 2 /2ρ l μ l ) 1/3 and the thickness of the boundary layer b = (6μ l /ρ l U) 1/2 . Note that if the Reynolds number is defined as Re = ρ l UR 2 f /μ l , (as UR f is a velocity scale) then b/R f ∼ Re −1/2 , similar to the Blasius boundary layer. At the interface, u Γ = Ur and the maximum value at r = R f is Using the same parameters as above gives u max Γ = 0.5 m s −1 , similar to the value found in COMSOL simulations in figure 6(c), giving us confidence in the prediction of such large speeds.
We now check the self-consistency of our assumptions. The ratio of the Couette and Poiseuille terms in (2.12) is where (8.4) was used. Using the same parameters as above and μ v = 1.7 × 10 −5 Pa s gives ≈ 0.12, sufficiently small for the approximation to be appropriate. The value of Re is ≈ 4200 1, so the boundary layer is indeed thin (as also seen in figure 6c).

Pressure-driven flow
We can now estimate the pressure-driven flow by combining Bernoulli's equation p v + ρ l u 2 /2 = const. along the free surface (a streamline) with (8.1), where we again assume a uniform film thickness, to find that Then using our expression for the film height (8.3), we obtain Interestingly, this is only dependent on geometry and gravity and totally independent of any liquid or vapour properties! The maximum value at r = R f , for R f = 3 mm and V as for a sphere of radius 3.5 mm is 0.5 m s −1 , the same as for the shear-driven component. Again, this gives confidence that the velocities observed in figure 6(c) are plausible, although the pressure-driven component is more sensitive to the interface shape and is likely overestimated by this analysis, suggesting that the shear stress component may be dominant.

Potential reasons for the discrepancy between computations and experiments
Experimentally, it is seen that larger drops are far more complex than their smaller counterparts, with the appearance of instabilities that can initiate cell-like vortices within apparently axisymmetric-shaped drops (Bouillant et al. 2018), drive vertical (Bouillant et al. 2018) and/or star shaped oscillations (Ma et al. 2017) of the drop surface, and finally rupture the drop due to the previously discussed chimney mechanism. Understanding the physical mechanisms for the formation of these structures is not so simple, as for larger drops additional physical effects come into play. Here, we discuss the two main assumptions in our model and assess, if relaxed, their potential to drive dominant flow features.

Temperature gradients in the droplet
Experiments show temperature variations T within Leidenfrost droplets of a few Kelvin (Wciślik 2016;Bouillant et al. 2018;Yim et al. 2020;Ma et al. 2017) which have the potential to generate thermal Marangoni effects and/or buoyancy-induced convection. These effects will drive internal flows that could alter the shape of the droplet and thus the geometry of the vapour layer, if included in our computational model. To capture these effects numerically requires the energy equation to be solved in the liquid and the dynamics of the air/vapour phase mixture to be computed and coupled to the drop through appropriate boundary conditions. Such an approach has been considered in, for example, (Pan et al. 2000) for non-Leidenfrost drops, but is beyond the scope of this article, where we focus instead on estimates for the scale of the various additional effects.
Marangoni effects -theoretical estimates for the flow driven by these forces are known to overpredict by many orders of magnitude experimental findings (Hu & Larson 2005;Yim et al. 2020), and this has been attributed to impurities/contamination in the liquid. Furthermore, Marangoni forces would enhance the internal flow shown in figure 6(c), with fluid pulled up the free surface from regions of low surface tension at the base, where the drop is hottest, towards the apex, where it is coldest and the surface tension is larger. Therefore, it is likely that including Marangoni effects would lead to an increased discrepancy between theory and experiments.
Buoyancy-driven flow -opposes the internal flow from figure 6(c), with fluid transported up along the axis of symmetry from the low density (hot) base towards the (colder) apex (e.g. see figure 1 in Yim, Bouillant & Gallaire 2021). To investigate further the potential for buoyancy effects to nullify the internal flow seen in figure 6(c), we estimate the typical velocity scales associated with this effect considering an inertially dominated regime, as the liquid Reynolds number based on flow speeds (∼10 cm s −1 ) in figure 6 is greater than 10 2 . To do so, we determine how the buoyancy force (per unit volume) ρ 0 βg T, where β = 7 × 10 −4 K −1 is the thermal expansion coefficient of water at the boiling temperature and ρ 0 is a characteristic density, can drive a characteristic acceleration U 2 /R to obtain a velocity scale of U = √ βRg T which based on T = 10 K gives U ∼ cm s −1 for the drop in figure 6(c), which is consistent with values seen in Yim et al. (2021). The buoyancy-driven flow is not insignificant, but it is still an order of magnitude below those velocities seen in figure 6, so it seems unlikely that the inclusion of these effects would drastically alter the drop shape.

Non-axisymmetric flow
Whilst Marangoni-and buoyancy-induced effects seem relevant, neither appear to overcome the observed discrepancies. Therefore, it seems most likely that, as experiments (Bouillant et al. 2018) indicate, in this regime the flow inside the drop loses its axial symmetry and becomes more complex, while our computational model assumes axisymmetric motion. It is possible that the azimuthally averaged flow is much weaker than axisymmetric simulations indicate, which naturally reduces the ability of the vapour flow to drive a huge internal flow and alter the geometry of the film. Recent articles have probed this non-axisymmetric flow by conducting linear stability analyses about a base flow (Yim et al. 2020(Yim et al. , 2021, but at present these do not capture vapour-driven internal flow. Clearly, in future work, there is scope for us to (i) use our model to establish an axisymmetric base state and then perturb it in the azimuthal direction and (ii) extend our model to capture three dimensional flow.

Concluding remarks
A computational model has been developed that allows us to go beyond the theory of Sobac et al. (2014) by taking into account internal flow and drop dynamics. Comparisons with Sobac et al. (2014) identified discrepancies that were subsequently corrected in Sobac et al. (2021) 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. (2012), 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. 2012;Shirota et al. 2016). Importantly, as for drop impact in isothermal conditions , 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.
observed experimentally in such regimes (Bouillant et al. 2018), but these aspects are beyond the scope of this article.
Convergence of all features of the drop were confirmed under spatial and temporal refinement studies, with results presented in the article graphically indistinguishable from those obtained under further refinement.
A.2. Numerical approach to solving the theoretical model The approach for the upper surface mimics that presented in Duchemin et al. (2005), but with the full axisymmetric expression for the curvature used here and iterations over the value of the curvature at the apex used to match a requiredR max . The equations are solved using a fourth-order Runge-Kutta method, starting from the apex of the drop and ending at the patching pointr =R p .
In order to compute the film shape, the governing one-dimensional equation (3.5) is discretized using a spatially second-order finite-difference method with time stepping based on the second-order Adams-Bashforth method. Here, the simulations are performed using the dimensionless grid size 10 −3 and the dimensionless time step 10 −5 ; chosen after conducting grid and time-step sensitivity tests. To solve (3.5), as mentioned in the main text, we require four boundary conditions: dh/dr = dκ/dr = 0 at the axis of symmetrỹ r = 0, and matching of dh/dr andκ atr =R p with the corresponding values obtained from the upper surface of the drop. Then, the procedure for the film shape is to start with a constant initial thickness of the vapour filmh(r, 0) =h 0 and iterate towards a steady state. After obtaining a steady-state solution satisfying ∂h/∂t = 0, ensuring the continuity of h(r) itself simply amounts to translating vertically the upper surface of the drop to reveal the complete equilibrium shape of the drop. The location of the patching point was varied from dz/dr = 1 to 10 with a negligible influence on the results.
Notably, our numerical results agree with those in Sobac et al. (2021).

Appendix B. Going beyond lubrication theory for small spherical drops
Consider an evaporating rigid sphere of radius R at distance h above an infinite hot planar solid surface. While we no longer assume that the gap between the sphere and the surface is thin compared with the size of the sphere, we retain a few other assumptions of the lubrication approximation we have used, namely, Stokes flow of the vapour and uniformity of its properties in space, and convective heat transfer is still neglected. In addition, since gas flow has to be considered in an infinite space, boundary conditions at infinity are important and we choose them as no flow and T = T w . Define dimensionless cylindrical coordinates (r,z) in units of R and the dimensionless temperatureT via It obeys the Laplace equation and the boundary conditions areT = 0 on the sphere andT = 1 on the solid surface and at infinity. The solution has onlyh = h/R as the parameter T =T(r,z;h).