Maximum size of drops levitated by an air cushion

Maximum size of drops levitated by an air cushion

Jacco H. Snoeijer, Philippe Brunet and Jens Eggers Department of Mathematics - University of Bristol, University Walk, Bristol BS8 1TW, United Kingdom
Physics of Fluids Group and J. M. Burgers Centre for Fluid Dynamics, University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands
Laboratoire de Mécanique de Lille UMR CNRS 8107, Boulevard Paul Langevin 59655 Villeneuve d’Ascq Cedex, France
September 15, 2019

Liquid drops can be kept from touching a plane solid surface by a gas stream entering from underneath, as it is observed for water drops on a heated plate, kept aloft by a stream of water vapor. We investigate the limit of small flow rates, for which the size of the gap between the drop and the substrate becomes very small. Above a critical drop radius no stationary drops can exist, below the critical radius two solutions coexist. However, only the solution with the smaller gap width is stable, the other is unstable. We compare to experimental data and use boundary integral simulations to show that unstable drops develop a gas “chimney” which breaks the drop in its middle.

I Introduction

Drops levitated on an air cushion have numerous applications, and have generated interest for a long time. For example, in lens manufacture drops of molten glass can be prevented from contact with a solid substrate Duchemin et al. (2005). This is achieved by levitating the glass above a porous mold, through which an air stream is forced. A second example is the so-called ’Leidenfrost’ drop Leidenfrost (1756), a drop of liquid on a plate hot enough to create a film of vapor between the drop and the plate Holter and Glasscock (1952); Goldshtick et al. (1986); Biance et al. (2003); Linke et al. (2006). Since the drop is thermally isolated insulated by the vapor film, it can persist for minutes Biance et al. (2003). Finally, a thin air film is believed to play a crucial role for the “non-coalescence” of a liquid drop bouncing off another liquid surface Couder et al. (2005); Gilet et al. (2007); Amarouchene et al. (2001).

The question we will address in this paper is whether for a given set of parameters, in particular the radius of the drop as it “rests” on the substrate, a stationary solution exists and whether it is stable. Apart from lense manufacture Duchemin et al. (2005), this question is important for the manipulation of corrosive substances Hervieu et al. (2001) or the frictionless displacement of drops Linke et al. (2006). Of particular interest is the maximum drop size that can be sustained, and the limit of very small flow rates. The drop continues to levitate in this limit since the gap between the liquid and the substrate becomes very small, so the lubrication pressure produced by the viscosity of the gas becomes significant. This enables us to employ asymptotic methods, making use of the disparity of scales between the gap size and that of the drop.

Experimentally, it is observed that the stability limit is reached when the radius equals at least a few capillary lengths . This natural length scale for our system is determined by the surface tension , density of the liquid, and acceleration of gravity . At a few capillary lengths, the drop is flattened to a pancake shape. Biance et al. Biance (2003); Biance et al. (2003) observed a critical radius


where is defined in Fig. 1. Beyond this radius, “chimneys” appeared, i.e. bubbles of air trapped below the curved and concave surface of the drop that rise owing to buoyancy, and eventually burst through the center of the drop. This suggests that the critical radius is related to the Rayleigh-Taylor instability of a heavy fluid (the drop) layered above a light fluid (the gas layer). In Biance (2003); Biance et al. (2003) this idea is used to estimate .

While this is close to the experimental value, the argument ignores the gas flow responsible for the levitation force. This flow was taken into account by Duchemin et al. Duchemin et al. (2005), who calculated the static shape of a drop levitated above a curved porous mould, using a combination of numerics and asymptotic arguments. For large drop volume, they found no physical solutions, while for smaller drops multiple solutions were calculated numerically. Duchemin et al. Duchemin et al. (2005) did not present a formal stability analysis of these solutions, but suggested that the limit of stability is related to the appearance of large oscillations on the underside of the drop.

A large number of studies of Leidenfrost drops have focused on the appearance of self-sustained oscillations of the drop Holter and Glasscock (1952); Watchers et al. (1966); Adachi and Takaki (1984); Takaki and Adachi (1985); Tokugawa and Takaki (1994); Strier et al. (2000); Biance et al. (2003); Snezhko et al. (2008). These oscillations can sometimes lead to a morphological bifurcation of the drop which takes the shape of a star Holter and Glasscock (1952); Adachi and Takaki (1984); Takaki and Adachi (1985); Strier et al. (2000); Snezhko et al. (2008). Similar star-shaped drops have been reported in drops vertically vibrated on non-sticky surfaces, and the shape is generally attributed to a parametric instability Yoshiyasu et al. (1996). Our original question was whether oscillations could perhaps be explained even in the limit of viscous drops, which we focus on in this paper. This is not the case, since both our asymptotic results and simulations of the complete dynamics (i.e. beyond linear stability analysis) show that once unstable, a drop breaks up owing to the formation of a chimney.

We treat both the liquid drop and the surrounding gas in the inertialess (Stokes) limit. For the asymptotic analysis, we also require the drop to be much more viscous than the gas. The main effect of this assumption is that there is hardly any flow inside the drop, so it can be treated as being in hydrostatic equilibrium at any instant in time. We also prescribe the rate at which gas is injected into the underside of the drop, thus ignoring the possible interplay between drop dynamics and vapor production in the Leidenfrost problem.

Our analysis is similar in spirit to the earlier paper of Duchemin et al. Duchemin et al. (2005), but we only address the simpler case of a flat substrate. As a result, we are able to perform all the calculations analytically (up to a few universal constants, which have to be computed numerically). Our solution curves are in qualitative agreement with those for a curved substrate Duchemin et al. (2005), but now imply a full analytical description. In addition, we are the first to perform a systematic stability analysis of the stationary states. We find the maximum stable radius


where goes to zero in the limit of vanishing gas flow.

For typical experimental flow rates we find that , consistent with the experimental result (1). At the end of the paper we discuss how our analysis relates to the stability argument of Biance (2003); Biance et al. (2003), based on the Rayleigh-Taylor instability.

Ii Problem formulation

ii.1 Geometry and dimensionless parameters

We consider axisymmetric drops of liquid, levitated above a flat surface by gas flowing into the underside of the drop, cf. Fig. 1. We set out to find the shape of stationary drops and their stability, as a function of the gas flow rate and the drop volume. The size of the drop is expressed by the Bond number


where is the volume of the liquid drop, and the unperturbed radius. The dimensionless gas flow rate supporting the drop reads


where is the volume of gas that escapes through the narrow neck region (see Fig. 1) per unit of time, and is the viscosity of the gas. Our analysis will identify the flux as the relevant quantity, which can be calculated by integrating the gas flux entering from underneath up to the neck position . Let us also introduce a slightly different dimensionalization of the flow rate


which will appear naturally in the analysis.

Finally, another parameter is the viscosity ratio between liquid and gas


but which will be considered asymptotically large for most of this paper. Throughout, lengths will be expressed in , velocities in , and stresses in .

Figure 1: Definitions and sketch of the matching regions.

ii.2 Structure of the problem

The problem we set ourselves is to solve the inertialess, axisymmetric fluid flow equations, with a prescribed influx of gas into the underside of the drop. Most of our analytical work assumes in addition that the drop is much more viscous than the gas. The structure of the expected solution is shown in Fig. 1. The gas pressure below the drop has to be sufficiently large in order to support the weight of the drop. In the limit of small dimensionless gas flux , the gap between the drop and the substrate must therefore be small in order to generate enough pressure. The underside of the drop inflates to a gas pocket, whose width is of similar size as the drop itself. The narrow gap is formed in a small neck region only, where a large curvature assures that the gas pressure can be sustained by corresponding surface tension forces. Apart from this viscous neck region the gas pressure is constant, both in the gas pocket as well as to the exterior of the neck.

This leads to the following asymptotic structure of the problem, characterized by the matching between three different regions. In the limit of small flux all viscous effects become localized in a small neck region, situated at a radius from the center. In this region, there exists a balance between viscous and surface tension forces. In addition, the slope of the gap profile turns out to be small in this region, so lubrication theory Batchelor (1967) permits to reduce the flow equations to an ordinary differential equation for . We will call this the inner solution or neck region.

To close the problem, boundary conditions are needed. These are provided by two outer regions on either side of the neck, denoted by (the gas pocket toward the center of the drop), and (the outside of the drop). Both regions are controlled by a balance of gravity and surface tension alone. First, we solve the equations in each of the regions individually. Second, we require that both the slope and the curvature of the profile match smoothly at the boundaries between two regions. This leads to a set of equations which determines stationary drop solutions uniquely. Solutions exist only below a certain critical neck position , in which case we find two branches, one with a small gap width (the lower branch), and an upper branch with a larger gap width.

Our stability analysis of the two branches is based on the observation that the relevant dynamic variable is the position of the neck, which can shift easily. The maximum stable neck radius does not coincide with , but is significantly smaller, located on the lower branch.

Iii Inner solution: neck region

iii.1 Lubrication approximation

We consider incompressible, axisymmetric flow in the gas layer, so that mass conservation gives


Here is the depth averaged horizontal velocity of the air layer, while is the rate at which air volume is injected per unit area below the drop. The main focus of the paper will be on stationary states and their stability. Stationary drop profiles are found by taking , and integrating (7) to


where is the flux in the lubrication layer. In the case where the injection source is localized at , the flux is simply constant.

To get a closed equation for in the neck region, we solve for . As our results will confirm, the neck region is shallow, , meaning that we can use the lubrication approximation Batchelor (1967) to analyze the flow, see Fig. 1. Owing to the large viscosity ratio between the drop and the surrounding gas, the liquid drop acts as a no-slip boundary, and the flow in the gas layer is well approximated by


Since the Reynolds number is very small in typical experiments, we use the Stokes approximation Batchelor (1967) to relate this velocity to the pressure. As there is almost no flow inside the drop, the liquid will be at hydrostatic equilibrium, . Furthermore, the pressure jump across the interface equals the curvature times the surface tension, so one obtains the lubrication pressure inside the gas layer as


In what follows we will show that the width of the neck region scales as and thus is asymptotically small in the limit of vanishing flux. We are therefore permitted to neglect the axisymmetric contribution to the curvature in the neck region. Using the horizontal component of the Stokes equation, , we find


Now (8),(11) provide a closed equation for the stationary interface profile :


The right hand side of (12) represents the viscous stress in the flow, and will only become important when is small, i.e. in a small neck region around , where we may set . This gives




A crucial observation is that there is no need to know the precise form of how the gas is injected, but one only requires the flux across the neck. This of course provides a great simplification for the Leidenfrost problem, where evaporation rates are related in a complicated way to the temperature profile inside the drop.

iii.2 Similarity solution for neck region

As gravity is unimportant in the thin neck region, (13) can be further simplified to


Since we are interested in the limit of small flux, we look for similarity solutions




In the limit , the solutions have to match onto a sessile drop of constant curvature. Since


one needs that for the curvature to remain finite as . Together with (17) this fixes and , and hence we have


The form of the similarity function will be determined from the matching below. The fact that justifies the assumptions made so far. First, we find that in the limit , justifying the use of lubrication theory. Similarly, , so that both gravity and the axisymmetric curvature can indeed be neglected in the neck region.

The asymptotic behavior of (17) is quadratic for both ,


Physically, the values of the asymptotic curvatures set the pressure in the corresponding outer regions.

Since (17) is of third order, solutions can be specified by three independent parameters, one of which can be absorbed into a shift of . Therefore, the two asymptotic curvatures uniquely determine the solution. As a consequence, the slopes are dependent variables. To perform the matching, we require the function


whose existence is assured by the above argument. Since (17) is invariant under the transformation and one must have


This function is computed numerically and is plotted in Fig. 2. We show below that stationary solutions correspond to the intersection of with another function , shown in the same figure. It can be seen that the matching breaks down at a critical neck radius , beyond which stationary solutions cease to exist.

Figure 2: Solid line: the function relates the slope to curvature of the inner solution (Eq. (23) with , corresponding to ). Dotted lines: the function provides the matching condition between inner region and gas pocket region (Eq. (51) with ). The three curves correspond to values (below critical), (critical), and (above critical).

Iv Outer solutions

Having seen that viscous effects are localized in the neck region, the rest of the drop is at static equilibrium. Hence, the pressure is constant both in the gas pocket between the drop and the substrate, as well as to the exterior of the neck. These pressures are not equal, however, since one requires a pressure difference to drive the flow across the neck. In Fig. 1 we therefore distinguish two outer regions, denoted by + and - respectively. Since , the outer solutions can be obtained from


where is the curvature of the interface. The constants determine the pressure difference across the neck


and will follow from the matching.

iv.1 Outer ”drop” region: +

Below we will find that the profile of the ”drop” region requires as , in order to match to the neck smoothly. This corresponds to a perfectly non-wetting sessile drop (Fig. 3). When matching the curvature we also require as . Owing to the vanishing slope near we are allowed to write in (24). Hence, one finds .

Figure 3: The outer solution of the ”drop” region corresponds to a perfectly non-wetting sessile drop. The size of the drop, characterized by , sets the curvature for the inner solution.

To deal with the overhang of the sessile drop, it is convenient to solve the profile in terms of the arclength along the interface. We define as the angle with the horizontal and rewrite (24) as


with boundary conditions


where is the value of the arclength at the top. Two of these five boundary conditions serve as the definitions of and , so that the remaining three boundary conditions fix the solution uniquely. The equations have been solved numerically.

Each value of thus gives a solution with a different , some of which are shown in Fig. 3. The numerically obtained relation between and is depicted in Fig. 4. For the maximal neck radius , introduced below, one finds .

The value of effectively sets the volume of the drop. Namely, the weight of the sessile drop is carried by the pressure exerted by the substrate on the contact area . The difference between the liquid and the gas pressures at is simply , so we find


where is defined as the real volume divided by , i.e.


Note that to obtain the real liquid volume, one has to subtract the volume of the gas pocket. However, goes to zero in the limit of vanishing gas flow, as shown below.

Figure 4: The outer solution fixes the value of as a function of the neck radius (which controls the drop volume). The maximal radius gives .

iv.2 Outer ”gas pocket” region: -

In the ”gas pocket” region, the profile is no longer multivalued and we can express the curvature as


The solution is then specified by (24) with boundary conditions


In Appendix A we show that the solution can be written as an expansion


where is a Bessel function of the first kind. Using furthermore that the curvature has to match the curvature of the inner solution , and thus , we can further simplify to


We see that the thickness scale of the gas pocket is set by the value of . In the limit of vanishing flux we expect this thickness to tend to zero, making a small parameter. To find solution branches, it is crucial to go beyond linear order and to find the term of order in (40). The only quantity that is needed to perform the matching to (22), coming from the inner solution, is the slope . This calculation is done in Appendix A.

At this point we can already infer an upper bound on the possible values of . Figure 5 shows the outer gas pocket solution (with normalized amplitude) for various values of . The outer solutions are defined on the domain where , hence the maximum possible neck radius is achieved when has its first minimum at the maximal radius (vertical line). The corresponding solution is drawn with a heavy line.

Figure 5: Outer solutions for the gas pocket region (amplitudes normalized to unity). Thin solid lines correspond to ; the dashed line, corresponding to , illustrates that would have to become negative to realize a neck radius larger than . The heavy line shows the maximum possible , corresponding to the first minimum of .

A second remark is that vanishes at , so must become zero at this radius. At even smaller radii turns negative, which yields negative values of . However, the inner solution cannot reach large for negative , which means that the matching procedure described here does not work. Dealing with this problem requires an additional matching region between the inner and outer (gas pocket) solution, the introduction of which lies beyond the scope of this paper. We will simply stay away from and instead focus on radii close to the maximal value , as detailed below.

V Matching the asymptotic regions

v.1 Matching conditions

We can now match the asymptotic regions by expressing (20),(21) in their original variables and expanding the outer solutions around ,


Therefore the matching conditions become


The conditions on the curvature were already taken into account when computing the outer profiles from (24). Typical values for are of order unity, while the slope requires as . This is why for the first outer solution we considered a perfectly non-wetting drop.

The conditions are more subtle. The thickness of the gas pocket goes to zero asymptotically so that both and will be small. In this case the selection of the solution explicitly requires the slope condition, which we express as


Together with (22) this closes the matching problem:


This equation indeed contains the three matching regions: implicitly depends on through the outer solution, is determined by the inner solution, while follows from the outer solution.

v.2 Bifurcation: critical radius

For a given value of the flux , we have reduced the problem to finding the intersections of the functions and . This is sketched in Fig. 2, showing and for and several values of . Depending on the value of , there can be two intersections, one intersection (when the curves are tangent), or no intersection. Each intersection corresponds to a stationary drop solution. This can be translated into a bifurcation diagram showing vs (Fig. 7). For small radii there are two branches of solutions, corresponding to the two intersections, which merge at . No stationary drop solutions exist for .

We analyze the bifurcation in the limit of vanishing flux, . We will show that


so that the critical neck radius approaches the maximal radius in the limit of vanishing flux. To analyze the vicinity of the critical point we introduce


At the same time we will find that . This means that as the limit of going to zero is reached, and , which implies according to Fig. 4. These two data fix the solution of (17) uniquely, and lead to the asymptotic profile shown in Fig. 6. From its minimum, one finds that


in agreement with the scaling found by Duchemin et al. (2005).

Figure 6: The inner solution obtained from numerical integration of (17), with and . It will follow that these values correspond to the critical solution. The minimum value determines the thickness of the neck (49).

We now analyze the first correction to the solution as increases, but in the limit where . This can be done by considering the corresponding limit of the functions and , cf. Fig. 2. Namely, the function approaches a constant, which is found numerically to be


on the other hand, the asymptotic form of becomes


The first term of (51) is found by expanding (40) for close to :


where we used the property . We need to keep the term as it can become of the same order as . For details we refer to Appendix A, where we show that .

The matching condition (cf. (46)) is now reduced to a horizontal line intersecting a cubic function:


Solving for , one finds


which has been plotted for different values of in Fig. 7. Thus for a given value of one finds two solution branches, which end at the critical value


as claimed before. Note that the smallness of the power makes non-negligible. For typical experimental values of the flux the critical point is thus substantially shifted with respect to the asymptotic value .

Plugging this back into (54), one finds the value of at the critical point: . But it follows from (40) that , and thus the maximum gap width is to leading order


This concludes the analysis of the stationary solutions, which are described by (54). At a given flow rate , the critical neck radius is given by (55), which approaches the maximal value in the limit .

Figure 7: The bifurcation diagram (), derived from (54). Curves correspond to different values of the flux, , revealing the weak dependence on . The dashed lines represent perturbations which are tangent to the solution curve. They represent marginal perturbations, separating stable from unstable solutions.

Vi Stability

We now turn to the important question of which part of the solution branches shown in Fig. 7 are stable. Essentially, we find that the lower branch is linearly stable, while the upper is linearly unstable. Surprisingly, however, the marginal point is not exactly at the maximum radius, but slightly before.

vi.1 Stability limit

Once more we make use of the fact that there is no flow inside the drop, so that the drop shape adjusts quasi-statically to variations in the neck region. We therefore consider infinitesimal variations in the neck position, and assess the corresponding change in levitation force . Since the pressure difference across the neck acts for , this force reads


A marginal perturbation occurs whenever the resulting levitation force is unchanged, , so that it still equilibrates the weight of the drop. Hence, we find the marginal condition


where the prime denotes the derivative with respect to . In order to produce the same levitation force, an increase in thus has to be compensated by a decrease of . Had the pressure stayed constant, would be larger than the weight of the drop leading to the formation of a chimney and thus to instability. Similarly, pressures smaller than the marginal condition leads to a stable situation, giving the stability criterion


In the limit of small , the pressure difference (25) is simply the difference of the curvatures:


so that stability requires


The derivative can be read off from Fig. 4, and is negative. Clearly, this has a stabilizing effect. The sign of can be inferred from the bifurcation diagram. The lower branch has a stabilizing contribution, while the upper branch is be destabilizing. The location of the marginal point, however, depends on the numerical values of the three terms.

Taking the derivative of (34), we find


Moreover, for vanishing flux , hence we may replace , giving the stability criterion


Near the maximal radius the criterion for stability becomes


Indeed, the upper branch with is unstable, but the marginal point is not at the maximum radius , but slightly before. This is indicated in Fig. 7 by the dashed lines, which each have a slope .

vi.2 Linear stability analysis

We now include dynamics into the stability analysis. We note first that an infinitesimal variation of the neck position, , also induces a variation of the curvature, , and of the flux, . These three parameters are related through mass conservation of the liquid and the gas. The analysis is closed by a third equation coming from matching the dynamic inner region to the hydrostatic outer regions.

The volume of the liquid can be calculated by the volume enclosed by the sessile drop solution, , minus the volume of the gas pocket , i.e.


This is exact up to asymptotically small corrections due to the inner region. The volume is (numerically) determined by the value of , while can be computed analytically using (40),


The expression simplifies at , because of the property . Since the liquid volume is strictly conserved, , one finds near


where the constant has been defined by (64). Relation (67) expresses the fact that when increases, increasing , the volume of the gas pocket has to increase by a similar amount to keep the liquid volume constant. This is achieved by an increase of .

Mass conservation of the gas is described by continuity (7), which can be integrated to


The second term on the right hand side can be identified as the rate of change of gas pocket volume , which we will write as . This change absorbs part of the injected air, decreasing the flux passing across the neck. Considering the radius somewhere inside the neck region, , the equation can be simplified to (using (11) and ):


where the variation of the flux reads


The matching condition (54) closes the dynamical system, taking into account the dependencies (67) and (70). The marginal case corresponds to a curve tangent to any of the lines shown in Fig. 7. Since in addition the slope of such a tangent curve must be according to (67), this uniquely fixes a point on any of the lines at constant . The critical tangent curve is drawn dashed in Fig. 7. Below this point, on the lower branch, solutions are stable, above they are unstable.

Formally, the growth rate of perturbations is computed by writing


Now using (67), (70), and the first variation of (54) one finds


The partial derivatives are to be evaluated from (54).

This indeed gives the same stability boundary as (63), which was based on a global force balance (note that ). The maximum stable radius is found by the condition , yielding


as an equation for at the stability boundary. This value of is inconsistent with the asymptotic estimate considered so far, indicating that the point where the solution exchanges stability is at a distance slightly larger from the critical point . This means that is smaller than expected (further down the lower branch, cf. Fig. 7). Thus the second term on the right of (73) is small compared to the other two, and we obtain


If we evaluate (54) in the same limit, we finally obtain


Thus for vanishing flux the maximum radius of stable solutions approaches , but with an even smaller power than . This scaling implies that , as seen in Fig. 7.

Vii Numerical tests

vii.1 Non-linear dynamical behavior

We begin with a simulation of the full axisymmetric Stokes problem, using a boundary integral method Rallison and Acrivos (1978), which has the advantage that it tracks the interface with high precision. The idea is to regard the interface as a continuous distribution of point forces, which point in the direction of the normal and whose strength is proportional to the mean curvature. Since for Stokes flow one knows the Green function giving the velocity field resulting from a point force, one can write the velocity anywhere in space as an integral over the free surface. In an axisymmetric situation, the angle integral can be performed, so the remaining integration is one-dimensional.

External flow sources can simply be added; in the present case we take the gas flow as a point source of strength situated at the origin on the solid plate which bounds the flow. For this a simple exact solution is available Happel and Brenner (1983). Likewise for the Green function one must take into account the presence of a no-slip wall. This is possible using the method of images Blake (1972), and the resulting boundary integral formulation has been applied successfully to the motion of drops relative to a wall Pozrikidis (1990). If, as in our case, the viscosity of the drop is different from that of the surrounding, one must account for the stress mismatch across the interface. This can be done at the cost of introducing another integral over the velocity on the interface into the equation, which turns the equation for the velocity field into an integral equation. After solving this equation for a given interface shape, the thus computed velocity field can be used to advance the interface.

We follow closely an earlier implementation of the boundary integral method, used for example the coalescence of two drops inside another fluid Eggers et al. (1999). The only significant difference is that the free space Green function has been replaced by those for half space, bounded by a wall. We tested the code by comparing to an exact solution of a sphere moving perpendicular to a wall Bart (1968). This is realized in the limit of a very small drop, or of very large drop viscosity, so that there is hardly any deformation. The agreement was good, but significant deviations occurred when the gap between the wall and the drop was smaller than 5 % of the drop radius. At present, we do not know the origin of this numerical problem, which prohibits us from investigating the asymptotic limit of very small gap spacings. Instead, we report on simulations at moderate gap spacings, which show the nonlinear stages of chimney formation, not captured by our linear stability analysis.

Figure 8: Boundary integral simulation of a drop with parameters , Bo=4.2, and . The drop relaxes toward a stable state, which is drawn as the heavy line.

Figure 8 shows a viscous drop which is slightly smaller then the stability boundary. Starting from a configuration shown as the light curve, it relaxes toward a stationary stable state (heavy line). For a Bond number which is just slightly larger, the same initial condition leads to a rising gas bubble in center of the drop, see Fig. 9. A thin film forms between the rising gas bubble and the top of the drop, which drains slowly.

Figure 9: The same as Fig. 8, but with a slightly larger Bo = 4.4. The air bubble under the center of the drop lifts up to form a chimney. The time interval between the profiles is = 3000, in units of .

As seen in Fig. 8, the neck radius is , giving . On the basis of our asymptotic theory (75), a rough estimate of the stability boundary gives .

vii.2 Lubrication approximation

To test the bifurcation scenario in more detail we resort to direct numerical simulation of the lubrication equation. Due to the overhang of the drop, we separate the upper part of the drop and the lower part of the drop at the maximum radius, , defined by the point . The upper part is solved as described in Sec. IV A, and for the lower part of the drop we use


This describes both the inner and outer regions in the lower part of the interface, while we have conveniently taken the rate of injection to be constant for all . Boundary conditions for this 3rd order equation are


where is the curvature at the point where the upper and lower solutions are patched. A one-parameter family of solutions is obtained through variation of the upper part of the drop. It was shown in Duchemin et al. (2005) that this procedure provides drop solutions that are quantitatively accurate.

The numerically obtained drop profiles are conveniently characterized by the position of the neck, , and the gap below the center of the drop, . Numerical results for the solution branches are shown as solid lines in Fig. 10 for two values of the flux. is a typical experimental value encountered for Leidenfrost drops, while illustrates the convergence toward the asymptotic limit. As predicted, there is a critical radius beyond which no stationary solutions exist. The asymptotic predictions shown in Fig. 7 have been translated to the dashed lines of Fig. 10. These are obtained from (54), where was computed from using (40). Good quantitative agreement is achieved for small enough values of the flux.

Finally, we determined the critical radius for a range of values of the flux . Figure 11 shows how the numerical values (dots) indeed approach the asymptotic prediction (solid line) in the limit of vanishing flux. Due to the very small powers , the convergence towards is extremely slow (horizontal line). As a consequence, the correction with respect to this asymptotic value will be significant for typical experimental values of the flux.

Figure 10: Bifurcation diagram vs for and . Smaller yield larger radii. Solid lines were obtained from numerical solution of the lubrication equation (77). Dashed lines correspond to asymptotic theory (54).
Figure 11: Critical radius as a function of the flux . The numerical values obtained from the lubrication equation (dots) indeed approach the theoretical prediction (55) (solid line). The dashed line indicates the asymptotic value .


Viii Discussion

Owing to the smallness of the neck region (49), we can make the simplification that the pressure inside the gas pocket below the drop is constant (Fig. 1). This pressure is larger than the atmospheric pressure and provides the force required to levitate the drop. Matching the pressure difference across the neck with the viscous flow then provides the bifurcation diagram of Fig. 7, yielding a critical neck radius . In the limit of vanishing flux, the critical radius approaches . This value arises because it is the first minimum of the function characterizing the shape of the gas pocket, which is the Bessel function . For larger , the gas pocket shape would need to become negative, which is of course not allowed.

Experimentally, the size of the drop is measured by looking at the drop from above. This measurement provides the maximum radius rather than the neck radius, cf. Fig. 1. For large puddles the difference between and approaches . For drop sizes relevant here we confirmed numerically that . Combined with (55), we thus find


for the boundary of stability, expressed in terms of the capillary length. Typical experimental values of can be extracted using (49), and typical experiments yield , obtained from diffraction data Biance et al. (2003). This gives , and thus , to be compared to reported experimental values of Biance (2003); Biance et al. (2003). A similar estimate of is obtained from considering the latent heat of evaporation Biance et al. (2003). Furthermore, our boundary integral simulations show that the nonlinear dynamics for larger drops lead to the formation of chimneys, as observed experimentally. We are therefore confident that the analysis in terms of Stokes flow provides an accurate description of this instability.

Let us now return to the argument put forward in Biance (2003); Biance et al. (2003), relating chimney formation to the Rayleigh-Taylor instability. The later occurs when a layer of fluid is suspended above another fluid of lower density, so that the system tends to destabilize due to buoyancy forces. Surface tension opposes this effect, so that the instability occurs at long wavelengths only. Biance et al. Biance (2003); Biance et al. (2003) propose that levitated drops remain stable as long as axisymmetric perturbations that fit inside the drop are stable with respect to this buoyancy driven instability.

For an infinitely extended liquid film, one finds that are axisymmetric eigenmodes, with the stability criterion . While the Bessel function does not have a well-defined period, the maximum drop size was estimated in Biance (2003) by the first minimum of the mode with , occurring at . In hindsight, our results justify this choice of taking the minimum of as the stability boundary, provided that it is identified with the neck radius, rather than with . With this connection, our results reduce to the Rayleigh-Taylor argument in the limit of vanishing gas flow, showing that the balance between buoyancy and surface tension provides the right mechanism. The effect of the gas flow is to slightly reduce the range of stable solutions (81).

Appendix A Gas pocket solution

In this appendix we expand the gas pocket solution for small amplitudes and compute the constant . We consider the equation


with boundary conditions


This is equivalent to solving


with boundary conditions


where We expand in ,


This yields a hierarchy of equations


with boundary conditions


The first equation gives , which can be inserted into the right hand side of the equation for .

To compute the constant , we require the ratio . In terms of the expansion


where we used the properties and . Comparing to (51), we simply find


We obtained this value numerically by solving the ODE for , for which we numerically obtained .


  • Duchemin et al. (2005) L. Duchemin, J. R. Lister, and U. Lange, J. Fluid Mech. 533, 161 (2005).
  • Leidenfrost (1756) J. G. Leidenfrost, De aquae communis nonnullis qualitibus tractatus, part2 (Duisburg, 1756).
  • Holter and Glasscock (1952) N. J. Holter and W. R. Glasscock, J. Acoust. Soc. Am. 24, 682 (1952).
  • Goldshtick et al. (1986) M. A. Goldshtick, V. M. Khanin, and V. G. Ligai, J. Fluid Mech. 166, 1 (1986).
  • Biance et al. (2003) A. L. Biance, C. Clanet, and D. Quéré, Phys. Fluids 15, 1632 (2003).
  • Linke et al. (2006) H. Linke, B. J. Alemán, L. D. Melling, M. J. Taormina, M. J. Francis, C. C. Dow-Hygelund, V. Narayanan, R. P. Taylor, and A. Stout, Phys. Rev. Lett. 96, 154502 (2006).
  • Couder et al. (2005) Y. Couder, S. Protière, E. Fort, and A. Boudaoud, Nature 437, 208 (2005).
  • Gilet et al. (2007) T. Gilet, N. Vandewalle, and S. Dorbolo, Phys. Rev. E 76, 035302(R) (2007).
  • Amarouchene et al. (2001) Y. Amarouchene, G. Cristobal, and H. Kellay, Phys. Rev. Lett. 95, 164504 (2001).
  • Hervieu et al. (2001) E. Hervieu, N. Coutris, and C. Boichon, Nucl. Eng. Design 204, 167 (2001).
  • Biance (2003) A. Biance, Gouttes inertielles: de la caléfaction a l’étalement (2003).
  • Watchers et al. (1966) L. H. J. Watchers, H. Bonne, and H. J. van Nouhuis, Chem. Eng. Sci. 21, 923 (1966).
  • Adachi and Takaki (1984) K. Adachi and R. Takaki, J. Phys. Soc. Jap. 53, 4184 (1984).
  • Takaki and Adachi (1985) R. Takaki and K. Adachi, J. Phys. Soc. Jap. 54, 2462 (1985).
  • Tokugawa and Takaki (1994) N. Tokugawa and R. Takaki, J. Phys. Soc. Jap. 63, 1758 (1994).
  • Strier et al. (2000) D. E. Strier, A. A. Duarte, H. Ferrari, and G. Mindlin, Physica A 283, 261 (2000).
  • Snezhko et al. (2008) A. Snezhko, E. Ben Jacob, and I. S. Aranson, New J. Phys. 10, 043034 (2008).
  • Yoshiyasu et al. (1996) N. Yoshiyasu, K. Matsuda, and R. Takaki, J. Phys. Soc. Jap. 65, 2068 (1996).
  • Batchelor (1967) G. K. Batchelor, An introduction to Fluid Dynamics (Cambridge University Press, 1967).
  • Rallison and Acrivos (1978) J. M. Rallison and A. Acrivos, J. Fluid Mech. 89, 191 (1978).
  • Happel and Brenner (1983) J. Happel and H. Brenner, Low Reynolds Number Hydrodynamics (Martinus Nijhoff Publishers, 1983).
  • Blake (1972) J. R. Blake, Proc. Camb. Phil. Soc 70, 303 (1972).
  • Pozrikidis (1990) C. Pozrikidis, J. Fluid Mech. 215, 331 (1990).
  • Eggers et al. (1999) J. Eggers, J. R. Lister, and H. A. Stone, J. Fluid Mech. 401, 293 (1999).
  • Bart (1968) E. Bart, Chem. Engin. Sci. 23, 193 (1968).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description