The thermo-wetting instability driving Leidenfrost film collapse

TL;DR

This paper develops a physics-based model for the Leidenfrost point, identifying vapor film instability mechanisms and deriving an expression that matches experimental data across various fluids and conditions.

Contribution

It introduces a linear stability analysis revealing the vapor film collapse mechanism and provides an ab-initio formula for the Leidenfrost point applicable to different fluids and surfaces.

Findings

The vapor film instability is driven by van der Waals forces.

The derived formula accurately predicts the Leidenfrost point.

A dimensionless number encapsulates the instability for wetting surfaces.

Abstract

Above a critical temperature known as the Leidenfrost point (LFP), a heated surface can suspend a liquid droplet above a film of its own vapor. The insulating vapor film can be highly detrimental in metallurgical quenching and thermal control of electronic devices, but may also be harnessed to reduce drag and generate power. Manipulation of the LFP has occurred mostly through experiment, giving rise to a variety of semi-empirical models that account for the Rayleigh-Taylor instability, nucleation rates, and superheat limits. However, a truly comprehensive model has been difficult given that the LFP varies dramatically for different fluids and is affected by system pressure, surface roughness and liquid wettability. Here, we identify the vapor film instability for small length scales that ultimately sets the collapse condition at the Leidenfrost point. From a linear stability analysis, it is shown that the main film stabilizing mechanisms are the liquid-vapor surface tension, viscous transport of vapor mass, and evaporation at the liquid-vapor interface. Meanwhile, van der Waals interaction between the bulk liquid and the solid substrate across the vapor phase drives film collapse. This physical insight into vapor film dynamics allows us to derive an ab-initio, mathematical expression for the Leidenfrost point of a fluid. The expression captures the experimental data on the LFP for different fluids under various surface wettabilities and ambient pressures. For fluids that wet the surface (small intrinsic contact angle), the expression can be simplified to a single, dimensionless number that encapsulates the thermo-wetting instability governing the LFP.

Figures

Tables

	${\sigma }_{LV}$	${\mu }_V$	$k_V\left(T_w-T_s\right)$	$h_{LV}$	${\rho }_L$	$g$	$P_V$
$M$	1	1	1	0	1	0	1
$L$	0	-1	1	2	-3	1	-1
$\tau$	-2	-1	-3	-2	0	-2	-2

	${\sigma }_{LV}$	${\mu }_V$	$k_V\left(T_w-T_s\right)$	$h_{LV}$
$M$	1	1	1	0
$L$	0	-1	1	2
$\tau$	-2	-1	-3	-2

Equations

$$\text{partialderivative}{\bar{u}}_{V}x+\text{partialderivative}{\bar{v}}_{V}y=0$$

$${\bar{u}}_V\text{partialderivative}{\bar{u}}_{V}x+{\bar{v}}_V\text{partialderivative}{\bar{u}}_{V}y=-\frac{1}{{\rho }_V}\text{partialderivative}\bar{\Phi }x+\frac{\Delta \rho g}{{\rho }_V}+\frac{{\mu }_V}{{\rho }_V}\text{partialderivative}[2]{\bar{u}}_{V}y$$

$$\text{partialderivative}\bar{\Theta }t+{\bar{u}}_V\text{partialderivative}{\bar{\Theta }}_{V}x+{\bar{v}}_V\text{partialderivative}{\bar{\Theta }}_{V}y=\frac{k_V}{{\rho }_V{c}_{p,V}}\text{partialderivative}[2]{\bar{\Theta }}_{V}y$$

$$\text{at}y=0,{\bar{u}}_V={\bar{v}}_V=0,{\bar{\Theta }}_V=1$$

$$\text{at}y=\bar{\delta },{\bar{u}}_V+{\bar{v}}_V\text{derivative}\bar{\delta }x=0,{\bar{\Theta }}_V=0$$

$${\rho }_V\text{partialderivative}\bar{\delta }t+\text{partialderivative}x{\int }_0^{\bar{\delta }}({\rho }_V{\bar{u}}_{V}dy)=-\frac{k_V}{h_{LV}}\text{partialderivative}{\bar{\Theta }}_{V}y{|}_{y=\bar{\delta }}$$

$$\text{partialderivative}{u}_{V}x+\text{partialderivative}{v}_{V}y=0$$

$$u_V\text{partialderivative}{u}_{V}x+v_V\text{partialderivative}{u}_{V}y=\frac{\Delta \rho g}{{\rho }_V}+\frac{{\mu }_V}{{\rho }_V}\text{partialderivative}[2]u_{V}y$$

$$u_V\text{partialderivative}{\Theta }_{V}x+v_V\text{partialderivative}{\Theta }_{V}y=\frac{k_V}{{\rho }_V{c}_{p,V}}\text{partialderivative}[2]{\Theta }_{V}y$$

$$u_V=\frac{\Delta \rho g{\delta }^2}{2{\mu }_V}(\eta -{\eta }^2)$$

$${\Theta }_V=\frac{T_V-T_s}{T_w-T_s}=1+(-1-\frac{1-c}{2})\eta +\frac{1-c}{2}{\eta }^3$$

$$\frac{1}{3}\frac{c_{p,V}\Delta T}{h_{LV}}c\left(1-\frac{3}{10}(1-c)\right)=1-c$$

$$\delta =2{\left(1-\frac{1}{3}\frac{c_{p,V}\Delta T}{h_{LV}}\right)}^{1/4}{\left(\frac{x\Delta T{\mu }_V{k}_V}{{\rho }_V{h}_{LV}g\Delta \rho }\right)}^{1/4}$$

$$\text{partialderivative}{u}^{\prime }x+\text{partialderivative}{v}^{\prime }y=0$$

$$u^{\prime }\text{partialderivative}ux+u\text{partialderivative}{u}^{\prime }x+v^{\prime }\text{partialderivative}uy+v\text{partialderivative}{u}^{\prime }y=-\frac{1}{{\rho }_V}\text{partialderivative}{\Phi }^{\prime }x+\frac{{\mu }_V}{{\rho }_V}\text{partialderivative}[2]u^{\prime }y$$

$$\text{partialderivative}{\Theta }^{\prime }t+u^{\prime }\text{partialderivative}\Theta x+u\text{partialderivative}{\Theta }^{\prime }x+v^{\prime }\text{partialderivative}\Theta y+v\text{partialderivative}{\Theta }^{\prime }y=\frac{k_V}{{\rho }_V{c}_{p,V}}\text{partialderivative}[2]{\Theta }^{\prime }y$$

$$\text{at}y=0u_V^{\prime }=v_V^{\prime }=0,{\Theta }_V^{\prime }=0$$

$$u_V^{\prime }{|}_{\eta =1}+\frac{{\delta }^{\prime }}{\delta }\text{partialderivative}u\eta {|}_{\eta =1}=0,{\Theta }_V^{\prime }{|}_{\eta =1}+\frac{{\delta }^{\prime }}{\delta }\text{partialderivative}{\Theta }_V\eta {|}_{\eta =1}=0$$

$${\rho }_V\text{partialderivative}{\delta }^{\prime }t+\text{partialderivative}x{\int }_0^{\delta }{\rho }_V{u}_V^{\prime }dy+\text{partialderivative}x{\int }_0^{\delta +{\delta }^{\prime }}{\rho }_V{u}_{V}dy=-\frac{k\Delta T}{h_{LV}}\left({\delta }^{\prime }\text{partialderivative}[2]\Theta y{|}_{\delta }+\text{partialderivative}{\Theta }^{\prime }y{|}_{\delta }\right)$$

$$u_V^{\prime }=a_0^{\prime }+a_1^{\prime }\eta +a_2^{\prime }{\eta }^2$$

$$\text{at}\eta =0u_V^{\prime }=0\to a_0^{\prime }=0$$

$$a_1^{\prime }=\left(\frac{\Delta \rho g\delta {\delta }^{\prime }}{2{\mu }_V}-\frac{\text{partialderivative}{\Phi }^{\prime }x{\delta }^2}{2{\mu }_V}\right)$$

$$a_2^{\prime }=\frac{\text{partialderivative}{\Phi }^{\prime }x{\delta }^2}{2{\mu }_V}$$

$$\text{partialderivative}{p}^{\prime }x=-{\sigma }_{LV}\text{derivative}[3]{\delta }^{\prime }x$$

$$\bar{\varphi }=\frac{A}{6\pi {\bar{\delta }}^3}$$

$$\text{partialderivative}{\varphi }^{\prime }x=-\frac{A}{2\pi {\delta }^4}\text{derivative}{\delta }^{\prime }x$$

$$\text{partialderivative}{\Phi }^{\prime }x=\text{partialderivative}{p}^{\prime }x+\text{partialderivative}{\varphi }^{\prime }x=-{\sigma }_{LV}\text{derivative}[3]{\delta }^{\prime }x-\frac{A}{2\pi {\delta }^4}\text{derivative}{\delta }^{\prime }x$$

$${\Theta }_V^{\prime }=b_0^{\prime }+b_1^{\prime }\eta +b_2^{\prime }{\eta }^2+b_3^{\prime }{\eta }^3$$

$$\begin{gathered} \frac{\delta }{4}\text{partialderivative}{b}_1^{\prime }t+\frac{1}{4}(1-2b_3)\text{partialderivative}{\delta }^{\prime }t+(\frac{1}{6}-\frac{2b_3}{15})\frac{\Delta \rho g{\delta }^2}{2{\mu }_V}\text{derivative}{\delta }^{\prime }x+(\frac{1}{12}-\frac{b_3}{20})(\frac{{\delta }^3}{2{\mu }_V})\left({\sigma }_{LV}\text{derivative}[4]{\delta }^{\prime }x+\frac{A}{2\pi {\delta }^4}\text{derivative}[2]{\delta }^{\prime }x\right) \\ +\frac{a_1}{30}(1-2b_3)\text{derivative}{\delta }^{\prime }x+\frac{a_1\delta }{20}\text{derivative}{b}_1^{\prime }x=\frac{3k_V}{{\rho }_V{c}_{p,V}{\delta }^2}{\delta }^{\prime }-\frac{3k_V}{{\rho }_V{c}_{p,V}\delta }{b}_1^{\prime } \end{gathered}$$

$${\rho }_V\text{partialderivative}{\delta }^{\prime }t+\frac{{\rho }_V\Delta \rho g{\delta }^2}{4{\mu }_V}\text{derivative}{\delta }^{\prime }x+\frac{{\rho }_V{\delta }^3{\sigma }_{LV}}{12{\mu }_V}\text{derivative}[4]{\delta }^{\prime }x+\frac{{\rho }_{V}A}{24\pi {\mu }_V\delta }\text{derivative}[2]{\delta }^{\prime }x+\frac{3k_V\Delta T}{h_{LV}{\delta }^2}{\delta }^{\prime }-\frac{2k_V\Delta T}{h_{LV}\delta }{b}_1^{\prime }=0$$

Full text

1. Introduction

As a surface is superheated above the boiling point of an adjacent fluid, vapor bubbles nucleate and grow. The boiling behavior of the liquid phase undergoes a fundamental change at a critical temperature known as the Leidenfrost point. Beyond this point, a film of insulating vapor forms between the liquid and the surface that suppresses heat transfer from the solid material. This heat flux reduction can be highly detrimental in the quenching of metal alloys by extending cooling rates and precluding the desired increase in strength and hardness 1. Alternatively, film boiling may be used to promote drag reduction as well as enable power generation through self-propulsion 2; 3; 4. Thus, modulation of the LFP through fluid choice, surface texture and chemistry for the specific application is crucial 5.

On a fundamental level, the physical mechanism responsible for the LFP is still uncertain. Many theoretical frameworks have been used to characterize the Leidenfrost effect and estimate the LFP, including hydrodynamic instability 6; 1, superheat spinodal limits 3; 9, and the change of liquid wettability on the heated surface with temperature 10; 11. A thermocapillary model has also been proposed that attributes the film instability to fluctuations at micron length scales; however the analysis posits that the thermocapillary effect is the dominant destabilizing term, which does not explain the significant change in LFP on surfaces with different wettabilities 12. For example, the LFP of water can vary from $300\degree\text{C}$ for hydrophilic surfaces to $145\degree\text{C}$ on hydrophobic surfaces 13; 14.

In this work, we introduce a stability analysis of the vapor film at the nanoscale regime. The dominant destabilizing term arises from the van der Waals interaction between the bulk liquid and the substrate across a thin vapor layer. On the other hand, liquid-vapor surface tension driven transport of vapor and evaporation at the two phase interface stabilize the film. The competition between these mechanisms gives rise to a comprehensive description of the LFP as a function of both fluid and solid properties. For fluids that wet the surface, such that the intrinsic contact angle is small, a single dimensionless number (eqn. 48) can be derived that encapsulates the instability determining the Leidenfrost point.

It is noted that the literature has proposed different names for the critical temperature associated with a droplet levitating on a heated plate (Leidenfrost point) versus the critical temperature for vapor film formation in pool boiling (minimum film temperature). The Leidenfrost point has been shown to be equivalent to the minimum film boiling temperature for saturated liquids on isothermal surfaces 2. In this work, the term LFP will be used for both cases as a matter of convenience, with the understanding that no undercooling is applied to the liquid phase for pool boiling scenarios unless explicitly stated.

2. Film Instability

There are many approaches to examine the stability of a vapor film adjacent to a superheated wall in two dimensions. Models have been developed with a base solution imposing static equilibrium, where the interface is at the saturation temperature corresponding to the imposed, far field liquid pressure 12; 16. Here, we consider the thickness of the vapor film to be in dynamic equilibrium, as in a vertical plate configuration 17 or vapor under a droplet 5. This appears to be a more general analysis since under experimental settings, droplet levitation occurs over a film that is continuously replenished by evaporation and depleted through escape of the buoyant vapor phase. Similarly in a horizontal setup for pool boiling, bubbles pinch off the film, necessitating a nonzero rate of evaporation to sustain a constant mean film thickness 17.

Fig. 1 shows a chosen “model” problem of film boiling on a vertical plate 17; 3. The vapor forms a laminar layer at the wall, with evaporation at the liquid interface sustaining the buoyant transport of vapor mass away from the base of the plate. The surrounding liquid is saturated and motionless with its properties fixed at the saturation temperature. The properties of vapor are assumed to be constant at the superheated wall temperature.

2.1. Governing Equations

The mass, momentum and energy conservation equations in the vapor domain are given by the following boundary-layer equations:

[TABLE]

[TABLE]

[TABLE]

where the parameters $\mu$, $\rho$, $k$, g and $c_{p}$ represent the the dynamic viscosity, density, thermal conductivity, gravitational acceleration and specific heat of the fluid, respectively. The term $\Delta\rho=\rho_{\scriptscriptstyle L}-\rho_{\scriptscriptstyle V}$ represents the density difference, the subscripts $L$ and $V$ denote the liquid and vapor field, and the temperature has been normalized as $\bar{\Theta}=\frac{\bar{T}-T_{s}}{T_{w}-T_{s}}$, the difference between the temperature field and the saturation temperature $T_{s}$ at the interface over the difference between the wall temperature $T_{w}$ and $T_{s}$. Note that we take the liquid phase to be motionless due to its much greater viscosity relative to the vapor ($u_{\scriptscriptstyle L}(x,t)=0$), as well as isothermal at saturation temperature due to its larger thermal diffusivity ($\Theta_{\scriptscriptstyle L}=0$). This simplification allows for an analytical solution for the base state per Burmeister 17 .

The generalized pressure term $\bar{\Phi}$ takes into account both the fluid pressure arising due to surface tension as well as due to van der Waals interactions. To first order in the base solution, these terms are negligible since the liquid-vapor interface is assumed to be locally parallel to the wall 19; 20; this implies $\bar{\Phi}\approx 0+\Phi^{\prime}$, where the overbar variables represent the general solution, the unbarred variables denote the base solution, and the primed variables give the perturbed solution. Additionally, the momentum and temperature equations are modeled as steady in the base solution and only exhibit a time varying term in the linearized equations for the perturbations.

The boundary conditions at the superheated wall and the liquid-vapor interface at $\bar{\delta}(x)$ are given by:

[TABLE]

[TABLE]

Eqn. 4 enforces the temperature and no interfacial slip at the impermeable wall, while eqn. 5 ensures that the tangential component of velocity is continuous and that the temperature is at saturation along the liquid-vapor interface. As with Burmeister, we neglect the "blowing" of vapor toward the plate by assuming $v(x,t)\approx 0$, such that the preceding set of boundary conditions are sufficient to fully specify the problem. Otherwise, the normal component of velocity would also need to be fixed at the interface; this leads to a cubic rather than a parabolic estimate to the velocity field.

The vapor generated due to phase change at this interface is balanced by the streamwise rate of change of the vapor flow in the film and the growth of the film thickness in time:

[TABLE]

where $h_{\scriptscriptstyle LV}$ is the latent heat of vaporization. In the base state, the time variation of the film thickness is taken to be negligible, $\derivative{\bar{\delta}}{t}\approx 0+\derivative{\delta^{\prime}}{t}$.

2.2. Base Flow

Along with the boundary conditions (eqn. 4 and 5), the base flow equations are:

[TABLE]

[TABLE]

[TABLE]

The steady, developed solution can be found approximately by using an integral expansion method, which is described in full detail by Burmeister 17. After introducing the normalized variable $\eta=y/\delta$, we can determine the velocity ($u_{\scriptscriptstyle V}$) and temperature field ($\Theta_{\scriptscriptstyle V}$) in the base solution:

[TABLE]

[TABLE]

Note that the velocity profile is locally parabolic due to the buoyancy force, while the velocity variation in the streamwise direction under mass conservation is encapsulated in the $\delta(x)^{2}$ dependence and occurs on a much larger length scale. The value of $c$ can be found by solving the quadratic expression:

[TABLE]

and $\Delta T=T_{w}-T_{s}$. For typical Jakob numbers around $Ja=\frac{c_{p,\scriptscriptstyle V}\Delta T}{h_{\scriptscriptstyle LV}}=\frac{1}{10}$, $c$ can be found from a simplified linear equation $c=1-\frac{1}{3}\frac{c_{p,\scriptscriptstyle V}\Delta T}{h_{\scriptscriptstyle LV}}$. Using this approximation, the film thickness $\delta$ is described by:

[TABLE]

Note that to first order, the velocity and temperature fields as well as the film thickness are steady.

2.3. Linearized Equations

The base solutions for the velocity, temperature and film thickness are perturbed, yielding the following linearized equations:

[TABLE]

[TABLE]

[TABLE]

The boundary conditions at the wall are:

[TABLE]

At the perturbed interface location $\delta+\delta^{\prime}$, the tangential velocity and temperature conditions (eqn. 5) after applying the locally parallel approximation $\derivative{\delta}{x}\approx 0$ give:

[TABLE]

The phase change equation at the interface (eqn. 6) is linearized as:

[TABLE]

Analogous to the base solution, the perturbed velocity is expanded in powers of $\eta$.

[TABLE]

[TABLE]

The terms $a_{1}^{\prime}$ and $a_{2}^{\prime}$ can be found as functions of the fluid properties and the generalized pressure gradient $\partialderivative{\Phi^{\prime}}{x}=\partialderivative{p^{\prime}}{x}+\partialderivative{\phi^{\prime}}{x}$ from the momentum equation eqn. 15 evaluated at the wall ($\eta=0$) and the tangential velocity condition (eqn. 18).

[TABLE]

[TABLE]

The pressure gradient arises from the liquid-vapor surface tension $\sigma_{\scriptscriptstyle LV}$ at the two phase interface due to capillary pressure induced by local nonzero curvature:

[TABLE]

This implies that positive curvature corresponds to the center of curvature lying in the vapor domain, such that the vapor bulges into the liquid. Here, we also introduce the disjoining pressure term $\bar{\phi}$, which describes the van der Waals interaction between the fluid and the substrate:

[TABLE]

The streamwise derivative of this term is negligible in the base state under the locally parallel interface approximation. The Hamaker constant $A$ is typically positive, denoting attractive interactions between dipoles 21. The perturbed component is:

[TABLE]

This gives an expression for the perturbed, generalized pressure term evaluated at the liquid-vapor interface.

[TABLE]

Next, the expanded perturbed temperature is:

[TABLE]

From the wall boundary condition (eqn. 17) and energy conservation equation (eqn. 16) at $\eta=0$, we find that $b_{0}^{\prime}=b_{2}^{\prime}=0$. Similarly, the temperature condition (eqn. 18) and energy conservation (eqn. 16) at the interface $\eta=1$ leads to an expression for $b_{1}^{\prime}$ in terms of $\delta^{\prime}$:

[TABLE]

The time evolution for the perturbed $\delta^{\prime}$ as a function of $b^{\prime}_{1}$ follows from the linearized phase change expression (eqn. 19) and the expressions for $u_{\scriptscriptstyle V}^{\prime}$ and $\Theta_{\scriptscriptstyle V}^{\prime}$ (eqn. 20 and 28):

[TABLE]

The perturbation equations 29 and 30 give two homogeneous conditions for $\delta^{\prime}$ and $b_{1}^{\prime}$. The perturbations can now be expressed in terms of normal modes:

[TABLE]

[TABLE]

To avoid introducing new notation, we will represent the amplitudes without subscripts $\delta_{a}^{\prime}\rightarrow\delta^{\prime}$ and $b_{1a}^{\prime}\rightarrow b_{1}^{\prime}$. Here, $k$ is the wave number and $\omega$ the time rate of growth of the perturbation. We combine eqn. 29 and 30 to obtain a single equation with the coefficient $\delta^{\prime}$. To simplify the representation, we introduce the following dimensionless parameters:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

This leads to the general, characteristic equation for the temporal growth rate $i\omega$ of the perturbation after eliminating $b_{1}^{\prime}$ using eqn. 29 and 30:

[TABLE]

The marginal state occurs when the real part $\text{Re}(i\omega)=0$, separating zones of stability ($\text{Re}(i\omega)<0$), where the perturbation amplitude decays in time, from regions of instability ($\text{Re}(i\omega)>0$), where the base state becomes unstable (Fig. 2 a). Note that only three dimensionless numbers $Ja$, $\pi_{LP}$ and $\pi_{LB\sigma_{\scriptscriptstyle LV}}$ govern the stability of the perturbed solution. With the inclusion of van der Waals interactions, the buoyancy terms described by $\pi_{LB\sigma_{\scriptscriptstyle LV}}$ become negligible at nanoscale, as will be discussed shortly.

This analysis incorporated time variation and convective transport in the energy equation. We can obtain a simpler expression for the stability problem by neglecting these two terms:

[TABLE]

where the buoyancy term $\pi_{LB\sigma_{\scriptscriptstyle LV}}$ is explicitly shown to give the dimensionless traveling wave velocity, signifying that the base flow acts only to convect the perturbation and does not affect its growth. The full derivation of eqn. 38 is provided in the Supporting Information.

The diffusive expression (eqn. 38) is a good estimate to the full stability equation (eqn. 37) for small Jakob numbers, as demonstrated in Fig. 2 c. Since the dimensionless parameters are calculated from the vapor properties at the superheated wall temperature, the Jakob number is small ($Ja\lessapprox\frac{1}{10}$) for the vapor phase of most fluids, implying that the thermal energy imparted by the heated solid is predominantly consumed through the latent heat of phase change rather than as sensible heat in raising the temperature of the vapor. Physically, dropping the time varying term in the energy equation implies that the heat conduction time scale is much longer (quasi-steady) than that for the perturbation growth in the phase change equation.

The Leidenfrost point corresponds to the lowest, critical $\pi_{LP,\;\text{crit}}$ on the marginal stability curve, below which the flow becomes unconditionally stable for all values of $k^{\prime\prime}$. For eqn. 37, a good approximation for the critical $\pi_{LP}$ can be derived by noting that due to how we scaled the dimensionless parameters, $\pi_{LP,\;\text{crit}}$ occurs at $k^{\prime\prime}=1$. This leads to an algebraic equation for $\pi_{LP,\;\text{crit}}$ :

[TABLE]

Eqn. 39 was verified against a numerical solution to the full stability equation (eqn. 37) with $\text{Re}(i\omega)=0$, and was found to give the same solution for $\pi_{LP,\;\text{crit}}$ up to machine precision for all parameter sets tested (Fig. 2 b). For small Jakob numbers, Fig. 2 c shows that the critical $\pi_{LP}$ can also be estimated from the diffusive expression (eqn. 38).

[TABLE]

The implication of eqn. 40 as a good estimate is threefold. First, it means that the boundary layer approximation is not required if film boiling is assumed to be diffusion dominated in both energy and momentum via lubrication theory. Secondly, the change in $\delta(x)$ along the streamwise direction is only captured in the base solution and becomes insignificant in the perturbation equations. Lastly, perturbation growth is independent of the base flow, which carries vapor generated at the liquid-vapor interface out of the local control volume. This explains why the LFP is not found to be strongly dependent on the configuration of the experimental set up. Any configuration eventually takes the vapor out of the film by some buoyancy driven force (even if it is horizontal or upside down - vapor eventually finds its way up). The strength of that driving force would indeed depend on configuration, but it does not matter for the perturbation solution.

2.4. Stabilizing Terms

The diffusive approximation to the critical $\pi_{LP}$ (eqn. 40) reveals that the main stabilizing terms are the liquid-vapor surface tension $\sigma_{\scriptscriptstyle LV}$, the evaporative phase change that replenishes local vapor mass $\frac{k_{\scriptscriptstyle V}\Delta T}{h_{\scriptscriptstyle LV}\rho_{\scriptscriptstyle V}}$, and the viscous shear $\mu$ that reduces mass transport away from any given point in the vapor field.

The liquid-vapor surface tension acts as a restoring force against oscillatory modes imposed on the basic, locally-parallel solution. Positive curvature of the liquid-vapor interface with its center in the vapor region (curving into the liquid), induces high pressure locally with an adjacent low pressure zone due to the negative curvature of the continuous two phase interface. This creates a pressure gradient that attempts to restore the basic state by dampening all possible oscillatory frequencies.

Similarly, a perturbed interface that bulges into the vapor steepens the thermal gradient in the vapor film, triggering an increase in the rate of evaporation locally that restores the base state and vice versa. Viscous shear is larger for smaller film thicknesses, therefore reducing mass transport away from a local bulge into the vapor domain and enhancing transport away from a bulge into the liquid field; this also acts to dampen perturbed modes.

2.5. van der Waals Interaction

The heterogeneous Hamaker constant $A_{SVL}$ is used to characterize the van der Waals dispersion forces between an uncharged surface and an adjacent liquid separated by vacuum. It is incorporated into this analysis via the generalized pressure gradient (eqn: 27). From the diffusive expression of the perturbation stability (eqn. 38), it is shown that these dispersion forces between the liquid and solid substrate across the vapor film is the only destabilizing term for attractive interactions $A>0$. The film is unconditionally stable if the interaction is purely repulsive $A<0$. The former case holds in general for a liquid separated from a solid by a vacuum or an intermediate gas phase 22; 23.

The relationship between the heterogeneous Hamaker constant and the contact angle $\theta$ of the substrate has been derived using Lifshitz theory 24; 25 :

[TABLE]

where $H_{SVL}$ is the equilibrium contact separation between the solid substrate (S) and the liquid (L) across vacuum (V) and takes on values in the order of magnitude of $1$ nm for most materials. Eqn. 41 can thus be used to account for the effect of surface wettability on the stability of the perturbed solution. As the van der Waals interaction only plays a significant role for film thicknesses that have reached the same order of magnitude as $H_{SVL}$, we approximate the ratio $\frac{H_{SVL}}{\delta}\approx 1$.

In this nanoscale regime, the neutral curve described by the full perturbation solution (eqn: 37) is insensitive to $\pi_{LB\sigma_{\scriptscriptstyle LV}}$, which encapsulates the buoyancy force on the vapor film and is on the order of $1\mathrm{e}{-14}$. The diffusive expression of the perturbation equation (eqn. 38) has an explicit dependence on $\pi_{LB\sigma_{\scriptscriptstyle LV}}$ only in the imaginary part of the temporal growth rate, such that the marginal state predicted is completely agnostic to changes in $\pi_{LB\sigma_{\scriptscriptstyle LV}}$. This implies that the stability criterion (eqn. 39 or 40) can be applied to capture the Leidenfrost point on plates of arbitrary orientation, as the direction and magnitude of the gravitational field does not play a significant role in the instability mechanism examined. The further implications of the nanoscale regime on this analysis is presented in the Supporting Information.

3. Verification by experiment and simulation

To compare against experimental data, eqn. 40 for the critical $\pi_{LP}$ can be rewritten in terms of material properties:

[TABLE]

Substituting in eqn. 41, we obtain the corresponding expression for the LFP in terms of the intrinsic contact angle on the substrate.

[TABLE]

where we have defined two new dimensionless parameters that we will show to be significant:

[TABLE]

Eqn. 43 provides an explicit relationship between the intrinsic contact angle of a fluid on a substrate and the Leidenfrost point for the system. Since each fluid property ($k_{\scriptscriptstyle V}(T),\mu_{\scriptscriptstyle V}(T),\sigma_{\scriptscriptstyle LV}(T)$, etc.) is calculated at the superheated wall temperature, the left hand side of eqn. 43 is in general a nonlinear function of temperature. The temperature at which eqn. 43 is satisfied corresponds to the predicted LFP; this can be found numerically with the temperature and pressure dependent fluid properties available from databases like NIST and tabulations from literature 26; 27; 28.

3.1. Surface dependence of the LFP

The LFP for water has been demonstrated to vary dramatically with changes in the liquid wettability of the solid surface 13; 14. Fig. 3 shows that the diffusive prediction of the LFP (eqn. 43) accurately captures the relationship between the LFP and the contact angle as delineated by experiments 29; 14; 30; 13; 31. Physically, larger contact angles indicate a hydrophobic substrate, which exhibits less attractive van der Waals interaction with the bulk liquid and presents a smaller destabilizing effect to the vapor film; the LFP thus decreases to near the boiling point. Without considering van der Waals interaction between the liquid and substrate surfaces, such a relationship cannot be explained or predicted from first principles.

Further evidence of the significant role played by van der Waals forces in governing the LFP arises from X-ray imaging of the vapor film collapse 33. Images spanning the film lifespan between formation and collapse showed that film collapse on the macroscopic level is preceded by submicron length scale vapor film thicknesses where the bulk liquid appears to wet the substrate. Although instabilities on the micron scale and above perturb the liquid-vapor interface and induce frequent local contact between the liquid and solid, only when the vapor film becomes unstable on the smallest length scales where van der Waals interactions dominate will the film completely collapse.

Lastly, we note that our theoretical analysis predicts that the main effects governing the LFP operate in the nanoscale regime, which is accessible by molecular dynamics (MD). Figure 4 shows one of the boiling heat transfer simulations performed using LAMMPS 34 software to numerically determine the LFP and the corresponding intrinsic contact angle of the substrate. Details of the MD implementations are provided in the Supporting Information. Note that a vapor film forms when the liquid water adjacent to the bottom plate is heated above the LFP, whereas liquid contact with the solid surface is preserved below the LFP due to the attractive heterogeneous van der Waals interactions. The relationship between the contact angle and the Leidenfrost point of the SPC/E (extended simple point charge) water model is in good agreement with the diffusive prediction of the LFP (eqn. 43). These simulations show that vapor film stability is ultimately determined at the proposed nanometric length scale where fluid-surface van der Waals interactions cannot be discounted and where the effect of gravity driven instabilities is nonexistent.

3.2. Fluid dependence of the LFP

For most experimentally available data on the Leidenfrost point, the contact angle of the fluid on the substrate material is low, around $\theta\approx 20\degree$. Nonetheless, the Hamaker constant must be found to determine the equilibrium separation $H_{SVL}$. Although the assumption $\frac{H_{SVL}}{\delta}\approx 1$ is made, the base film thickness $\delta$ still needs to be incorporated into our instability expression via $\pi_{2}$ (eqn: 45). We can find the homogeneous Hamaker constant of the fluid (acentone, ethnanol, benzene, etc.) and the substrate (gold, aluminum, copper), and take the geometric mean to obtain the heterogeneous value 23; 35. From the relationship between the surface energy and homogeneous Hamaker constant, we can obtain the homogeneous contact separations via:

[TABLE]

[TABLE]

Figure S1 shows that it is possible to determine either the heterogeneous Hamaker constant given the Leidenfrost point for a fluid on a solid substrate, or vice versa with knowledge of the homogeneous Hamaker constants of both species.

In general, experimental data on the homogeneous Hamaker constants may not be available for a fluid or substrate of interest. Here, we note an avenue for simplification: it is observed that the dimensionless quantity $\pi_{2}=\frac{h_{\scriptscriptstyle LV}\rho_{\scriptscriptstyle V}\delta}{\sigma_{\scriptscriptstyle LV}}$ in the diffusive expression is around $0.06$ for most fluids at their respective Leidenfrost temperatures. This suggests that there exists a functional dependence $H_{SVL}=F(\frac{h_{\scriptscriptstyle LV}\rho_{\scriptscriptstyle V}}{\sigma_{\scriptscriptstyle LV}})$. Additionally, most experimental setups in the film boiling regime feature fluids that wet the surface in contact, such that their intrinsic contact angle are small ($\theta\approx 20\degree$) 28.

From the diffusive expression (eqn. 43) valid for low Jakob numbers, the above approximations leads to a simplified, dimensionless prediction to the Leidenfrost point for fluids/substrate systems with low, intrinsic contact angles:

[TABLE]

This dimensionless quantity also arises by application of the Buckingham’s Pi Theorem to the system, as discussed in the Supporting Information. Fig. 5 shows that the temperature at which this equality is satisfied captures the experimental data on the LFP for a variety of different fluids, including cryogens and liquid metals. The single dimensionless number describes the terms that stabilize the vapor film, including surface tension, phase change and viscous transport, while the critical value corresponding to the LFP denotes the destabilizing effect of attractive van der Waals interaction between the bulk liquid and solid substrate. Larger values of $\pi_{1}$ above the critical imply the system is the film boiling regime, since the stabilizing terms dominate.

3.3. Pressure dependence of the LFP

Experimental work has shown that the LFP depends on the ambient pressure applied, such that the Leidenfrost temperature gradually increases from near the boiling point towards the critical point of the fluid 46. Figure 6 demonstrates that the diffusive expression (eqn. 43) captures the LFP of various fluids for both subatmospheric and superatmospheric pressures up to the critical point.

4. Conclusion

The dynamic stability of a vapor film on a heated vertical wall under the effects of gravity were considered. The only possible instability at nanoscale was driven by attractive van der Waals interaction between the bulk liquid and the substrate, which could be stabilized by the liquid-vapor surface tension, evaporative phase change and viscous transport. The marginal or neutral state can be found analytically (eqn. 39) for the most general case, or simplified for small Jakob number flows to a diffusive approximation (eqn. 40). The resulting theoretical solution for the LFP captures the variation of experimental data with surface wettability, fluid properties and pressure.

A single, dimensionless number $\pi_{1}$ is found to encapsulate the physical instability mechanism of the Leidenfrost phenomenon for wetting fluids. The value of $\pi_{1}$ with respect to the critical denotes regimes in which the vapor film is stable or unstable, providing a useful characterization of both the thermodynamic state and the physical means by which transition to the pool boiling regime occurs.

This insight into the nanoscale mechanisms inducing the transition from film to nucleate boiling enables control of the phase adjacent to the surface 48. It would be of interest to extend the instability mechanism towards surface roughness, which experiment has shown to effect dramatic changes in the LFP beyond what can be explained by variation in surface wettability 5; 14; 49. In addition, a theoretical treatment of the Nukiyama temperature corresponding to the critical heat flux may reveal the mechanism underlying transition boiling and provide a comprehensive understanding of the entire boiling curve under a unified, physical framework.

Author contributions

N.A.P. and T.Y.Z. conceived and planned the research, performed the analyses, and wrote the manuscript.

Competing interests

The authors have no competing financial interests or other interests that might be perceived to influence the results and/or discussion reported in this paper.

Corresponding author

To whom correspondence should be addressed. E-mail: n-patankar@northwestern.edu

Supporting Information (SI)

Buckingham’s Pi Theorem

In general, the thermodynamic variables that affect the physical onset of film boiling include the thermal conductivity of the vapor ($k_{V}$), temperature ($T$), viscosity of the vapor ($\mu_{V}$), liquid-vapor surface tension ($\sigma_{LV}$), the specific latent heat of vaporization ($h_{LV}$), the liquid density $\rho_{L}$, the vapor pressure $P_{V}$, and the gravitational constant $g$ 1; 2; 3. This leads to four dimensions, namely mass $M$, length $L$, time $\tau$ and temperature $\phi$. Taking the product $k_{V}(\Delta T)=k_{V}(T_{w}-T_{s})$ as a single variable, the number of dimensions can be reduced to three, forming the dimensional matrix 4:

The repeating variables are chosen to be $\sigma_{LV}$, $k_{V}(T_{w}-T_{s})$ and $\rho_{L}$, which span the space of $M$, $L$, and $\tau$. Given seven variables and a rank three matrix, four dimensionless numbers are expected: $\Pi_{1}=\frac{\sigma^{2}}{k_{V}(T_{w}-T_{s})\mu_{V}}$, $\Pi_{2}=\frac{\sigma_{LV}^{2}h_{LV}}{k_{V}^{2}(T_{w}-T_{s})^{2}}$, $\Pi_{3}=\frac{\sigma_{LV}^{5}g}{k_{V}^{4}(T_{w}-T_{s})^{4}\rho_{L}}$, $\Pi_{4}=\frac{\sigma^{2}P_{V}}{k_{V}^{2}(T_{w}-T_{s})^{2}\rho_{L}}$. We note immediately that for the dimensionless numbers $\Pi_{1}$ and $\Pi_{2}$, the repeating variable $\rho_{L}$ is absent. This suggests that density or mass is not relevant to the subset of variables $k_{V}(T_{w}-T_{s})$, $\sigma_{LV}$, $h_{LV}$ and $\mu_{V}$; the dimensionless combination of these variables may be responsible for the static equilibrium behavior of a fluid in the different boiling regimes.

Thus, to reduce the number of independent variables describing the LFP, the gravity, pressure and density of either phase are postponed for consideration in the present analysis. Prior work has shown that accounting for gravity and the density ratio between liquid and vapor accurately captures the length scales of the droplet and vapor film as well as their evolution in time specific to the Leidenfrost regime 5. The present work focuses instead on the characteristics of a fluid that determines its LFP and distinguishes film boiling from nucleate or transition boiling.

By reducing the variable space, we obtain the dimensional matrix:

The rank of the matrix is two as the rows are linearly dependent; since density, pressure and gravity are neglected, mass is not relevant to the problem. With four variables and a rank two matrix, two dimensionless numbers are expected from the analysis. The repeating variables are chosen to be $\sigma_{LV}$ and $k_{V}(T_{w}-T_{s})$, which yield the dimensionless numbers:

[TABLE]

[TABLE]

The dimensionless number $\pi_{1}$ represents the ratio of the liquid-vapor surface tension to the viscosity, thermal conductivity and temperature of the vapor phase. The LFP is posited to increase for fluids with larger liquid-vapor surface tension based on the model for the maximum superheat limit 3 and the Taylor type instability 1. The same correlation between surface tension and the LFP has been demonstrated experimentally across a broad range of fluids 2. From the superheat limit viewpoint, the relationship arises since a lower liquid-vapor surface tension decreases the energy barrier and increases the rate of vapor nucleation, which reduces the critical temperature necessary to form a vapor film.

From the Taylor instability perspective, the role of surface tension appears contradictory. Larger surface tension is theorized to stabilize the liquid-vapor interface by suppressing low wavelength perturbations and reducing the growth of high wavelength disturbances 1. It is therefore expected that higher surface tension would lead to a lower LFP by promoting the formation of a stable vapor film. However, the predicted expression for the Leidenfrost temperature originating from Taylor instability theory proposes that the LFP should increase with the surface tension, based on a fitted expression for bubble radii 1. While accurate for n-pentanes and carbon tetrachloride, it has been shown experimentally that the LFP predicted from the instability analysis is less accurate for water, cryogenic fluids and liquid metals 2.

The experimentally corroborated relationship between liquid-vapor surface tension and the LFP of a fluid means that the Leidenfrost regime should occur for small values of the dimensionless number $\pi_{1}$, corresponding to low surface tension, high vapor viscosity, high vapor thermal conductivity and high temperature; nucleate or transition boiling occurs for large values of $\pi_{1}$. For consistency, this suggests that higher vapor thermal conductivity should trigger the Leidenfrost phenomenon by conveying more heat to the liquid-vapor interface and boosting the rate of vaporization. A higher vapor viscosity should also induce film boiling by stabilizing the vapor layer from perturbations and penalizing the escape of the vapor phase from the film.

Physically, $\pi_{1}$ compares two velocities $\frac{k_{V}(T_{w}-T_{s})}{\sigma}$ and $\frac{\sigma}{\mu_{V}}$; the units of these quantities are $\frac{L}{\tau}$. The former $\frac{k_{V}(T_{w}-T_{s})}{\sigma_{LV}}$ gives a velocity scale of the vapor phase in the direction normal to the liquid vapor interface, which is approximately parallel to the heated surface. The latter $\frac{\sigma_{LV}}{\mu_{V}}$ gives a velocity scale of the vapor phase in the direction parallel to the liquid vapor interface. If the normal component of the velocity $\frac{k_{V}(T_{w}-T_{s})}{\sigma_{LV}}$ exceeds the tangential component $\frac{\sigma_{LV}}{\mu_{V}}$, a film of vapor can be sustained beneath the droplet. The transition to the Leidenfrost regime therefore occurs at $\pi_{1}$ around order one.

The second dimensionless number $\Pi_{2}$ represents the ratio of the liquid-vapor surface tension and the latent heat of vaporization to the thermal conductivity and temperature of the vapor phase. The Leidenfrost regime should occur for small values of $\Pi_{2}$, corresponding to low surface tension, high vapor thermal conductivity, high temperature and low latent heat of vaporization. With all else constant, a lower latent heat can trigger the Leidenfrost effect by increasing the rate of vaporization to sustain a stable film 3.

Physically, $\Pi_{2}$ compares the specific energies $h_{LV}$ and $\frac{k^{2}_{V}(T_{w}-T_{s})^{2}}{\sigma_{LV}^{2}}$. The energy scale $h_{LV}$ gives the specific latent heat needed to sustain vaporization at steady state. The latter $\frac{k^{2}_{V}(T_{w}-T_{s})^{2}}{\sigma_{LV}^{2}}$ represents the energy scale for heat conducted to the liquid-vapor interface. If the conducted heat $\frac{k^{2}_{V}(T_{w}-T_{s})^{2}}{\sigma_{LV}^{2}}$ that reaches the interface exceeds the specific latent heat $h_{LV}$, a vapor film can be sustained beneath the liquid. The transition to the Leidenfrost regime may thus occur for $\Pi_{2}$ around order one.

Molecular dynamics

Molecular dynamics simulations were implemented using LAMMPS 34 software. Liquid water was initially equilibrated at saturation temperature in the canonical ensemble with constant pressure imposed by a piston constrained to move only in the direction orthogonal to the bottom surface 6. The solid plate was constructed using a graphene sheet, and the interaction between the solid substrate and the SPC/E water molecules was governed by the 6-12 Lennard Jones pair potential. After equilibration, the liquid adjacent to the bottom surface was heated to a target temperature around the LFP, whereas the liquid adjacent to the piston was maintained at constant, saturation temperature to simulate nonequilibrium heat transfer conditions 7.

Relationship between the Hamaker constant and the equilibrium separation

The heterogeneous contact separation have been estimated as the arithmetic or geometric mean of the homogeneous values 24, although these means may not provide a good approximation for the actual $H_{SVL}$ in general. Nonetheless, for our theory to be physically consistent, the heterogeneous contact separation corresponding to the experimental Leidenfrost temperatures must be between the two bounding homogeneous values. Figure S1 shows that for the alkane family, this condition $H_{SVS}<\delta\approx H_{PVL}<H_{LVL}$ is satisfied, and the arithmetic and geometric means provide a reasonable estimate to the actual heterogeneous value.

Drummond has shown that as the chain length of the alkane species increases, the contact angle increases correspondingly. This suggests that longer chain alkanes in the liquid phase have unfavorable energetic interactions with a given substrate (greater liquid-solid surface energy $\sigma_{SL}$) compared to small chain alkanes on the same solid material. The heterogeneous, equilibrium distance $H_{PVL}$ therefore tends to increase with the straight chain length of the alkane species, moving from near the solid separation $H_{SVS}$ toward the liquid value $H_{LVL}$.

4.1. Derivation of the diffusive equation

We begin with the general formulation of the vertical film boiling setup, taking into account only the diffusive terms in the steady momentum and energy conservation equations:

[TABLE]

[TABLE]

[TABLE]

where $f_{d}$ is a buoyancy driven force dependent on the system configuration. As before, the boundary conditions at the superheated wall and the liquid-vapor interface at $\bar{\delta}(x)$ are given by:

[TABLE]

[TABLE]

while the phase change equation at the interface is written as:

[TABLE]

The base solution for the velocity field as found using an integral expansion is identical to our prior solution to the boundary layer equations:

[TABLE]

The temperature field in the base solution is linear, as expected from a diffusive approximation:

[TABLE]

The corresponding perturbed solutions are:

[TABLE]

[TABLE]

The variation of the film thickness in the streamwise direction is found from eqn. S8 after assuming time invariance in the base state:

[TABLE]

Note that compared with the solution for the film thickness obtained from the boundary layer equations, eqn. S15 is missing the term $(1-\frac{1}{3}\frac{c_{PV}\Delta T}{h_{LV}})^{(}1/4)$. This is expected given only conduction is consumed in the thin film limit.

Linearizing the phase change conservation condition (eqn: S8) at the interface we obtain:

[TABLE]

The change in $\delta^{\prime}$ is obtained after substituting in the expressions for $u^{\prime}$ and $\Theta^{\prime}$ (eqn: S11 and S12 respectively):

[TABLE]

Expressing in terms of normal modes:

[TABLE]

where for brevity we take $\delta^{\prime}_{a}\rightarrow\delta^{\prime}$. This gives us:

[TABLE]

After substituting in the set of dimensionless numbers introduced previous, we recover the diffusive expression (eqn. 38).

4.2. Nanoscale implications

Given that the phenomenon is localized in the nanoscale regime, it is important to discuss the applicability of the Navier-Stokes equations in the context of O(1) in the Knudsen number (Kn). Hadjiconstantinou showed that the second order Knudsen layer correction is qualitatively robust well beyond $Kn\approx 0.4$, such that the underlying Navier-Stokes constitutive laws captures the behavior of arbitrary flows in spite of superimposed kinetic corrections in the flow field 8.

Here, the base solutions for velocity (eqn. S9) and temperature (eqn. S10) for the diffusive equation above are represented by the lowest order polynomials possible, given the buoyancy force in the momentum equation and no heat source in the energy equation. An examination of the system from a time and spatially averaged statistical standpoint would give us similar profiles with an effective slip and temperature jump at the interfaces, giving rise to effective viscosities and thermal conductivities that converge to the continuum values as the film thickness is increased.

It would be of great interest and rigor to frame the film stability problem in the context of a linearized Boltzmann analysis to accommodate larger Kn flows; however, the excellent agreement between the present, simplified model with experiment suggests the dominant effects governing the Leidenfrost film stability is captured despite the absence of statistical analyses describing intermolecular collisions. Additionally, the granularity of linear Boltzmann equation may obscure the effect of surface tension and van der Waals interactions, which have been shown to be important factors in pinpointing the LFP.