2.1 Region I: the cusp singularity

At very high speeds, the interface is bent so severely that it appears to meet the plate tangentially, at an apparent contact angle $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}$ . We describe this situation using the asymptotic solution of the contact line problem of Benney & Timson (Reference Benney and Timson1980) for $M=0$ and $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}$ , and assuming a no-slip boundary condition. This is appropriate for our case, as we are on an intermediate scale excluding the contact line region. Since the interface is parallel to the plate at the contact line, the contact line singularity is weakened and the local energy dissipation remains finite, as we will see. As a result, the paradox discovered by Huh & Scriven (Reference Huh and Scriven1971) does not exist, although the curvature does still diverge.

For $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}$ , since viscous stresses dominate gravitational effects about the contact line, the liquid phase is described by the Stokes equation (1.2) with $g=0$ . This is equivalent to solving the biharmonic equation for the streamfunction $\unicode[STIX]{x1D713}$ , represented in polar coordinates $(r,\unicode[STIX]{x1D719})$ , defined in figure 6(a),

(2.1a-c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}^{2}\unicode[STIX]{x1D713}=0,\quad u_{r}=\frac{1}{r}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}},\quad u_{\unicode[STIX]{x1D719}}=-\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}{\unicode[STIX]{x2202}r}, & & \displaystyle\end{eqnarray}$$

subject to the boundary conditions

(2.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}rcl@{}}\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}: & \quad & \unicode[STIX]{x1D713}=0,\quad u_{r}=Ca,\\ \unicode[STIX]{x1D719}=g(r): & \quad & \unicode[STIX]{x1D713}=0,\quad \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{n}=-\unicode[STIX]{x1D705},\quad \boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D748}\boldsymbol{\cdot }\boldsymbol{t}=0.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Here we have written velocities in units of the capillary speed $\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D702}$ , stress in units of $\unicode[STIX]{x1D6FE}/l_{c}$ , and $\unicode[STIX]{x1D705}$ , the curvature of the interface, in units of $1/l_{c}$ . This corresponds to the no-slip boundary conditions at the plate, the dynamical stress conditions at the liquid–vacuum interface, and a further condition that the liquid–vacuum interface and the liquid–solid interface are streamlines, where the dividing streamline is taken to be $\unicode[STIX]{x1D719}=0$ . Benney & Timson (Reference Benney and Timson1980) commit a sign error in front of the curvature term which does not affect their method of calculation, but of course invalidates their results. Ngan & Dussan V. (Reference Ngan and Dussan V.1984) noticed the sign error, but claim that the corrected results lead to conclusions which are ‘physically meaningless’. Mahadevan & Pomeau (Reference Mahadevan and Pomeau1999) claim to have a found a singularity-free solution, and incorrectly conclude that the interface shape should be regular in the case of a $180^{\circ }$ contact angle. Here we hope to set the record straight, and test our conclusions by direct comparison with our numerical simulations.

Figure 6.

(a) Sketch of the cusp solution for a $180^{\circ }$ contact angle. (b) Solutions of (2.9).

Following Benney & Timson (Reference Benney and Timson1980), we find a similarity solution for this problem, where the free surface, to leading order as one approaches the contact line, has the power-law form

(2.3) $$\begin{eqnarray}\displaystyle h(x)=ax^{q}. & & \displaystyle\end{eqnarray}$$

In order to produce a $180^{\circ }$ contact angle, we must have $q>1$ for consistency.

In the classical calculation of Huh & Scriven (Reference Huh and Scriven1971), (2.1) and (2.2) (with the exception of the normal stress balance) are solved in a wedge, such that $h(x)=x\tan \unicode[STIX]{x1D703}$ instead of (2.3). This gives a unique solution for the limit $r\rightarrow 0$ , and the normal stress balance will in general not be satisfied. The normal stress balance is then used to calculate corrections to the wedge-shaped interface in a perturbative fashion.

In the present calculation, we use the additional freedom of choosing the exponent $q$ in order to satisfy the normal stress balance as well; the value of $q$ then results from a solvability condition, which makes this an example of self-similarity of the second kind (Eggers & Fontelos Reference Eggers and Fontelos2015), as opposed to the Huh–Scriven problem, which is of the first kind. The constant $a$ in (2.3) sets the amplitude of the solution. In principle, it has to be determined by matching to an outer solution; we will find it below by comparing to the numerical solution.

In polar coordinates $(r,\unicode[STIX]{x1D719})$ , see figure 6, the surface becomes $\unicode[STIX]{x1D719}=g(r)\approx ar^{q-1}$ . The solution for $\unicode[STIX]{x1D713}$ , with a uniform flow $-Ca\,\boldsymbol{e}_{x}$ from the downwards velocity of the plate, has the similarity form

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}=-Ca\,r\sin \unicode[STIX]{x1D719}+r^{q}F(\unicode[STIX]{x1D719})+O(r^{2q}). & & \displaystyle\end{eqnarray}$$

This will ensure that, to leading order, the interface is the streamline $\unicode[STIX]{x1D713}=0$ , as long as $F(g(r))\simeq F(0)=a\,Ca$ is satisfied. For $\unicode[STIX]{x1D713}$ to be a solution of (2.1), $F$ has the form (if $q\neq 2$ )

(2.5) $$\begin{eqnarray}\displaystyle F(\unicode[STIX]{x1D719})=A\cos q\unicode[STIX]{x1D719}+B\sin q\unicode[STIX]{x1D719}+C\cos (q-2)\unicode[STIX]{x1D719}+E\sin (q-2)\unicode[STIX]{x1D719}, & & \displaystyle\end{eqnarray}$$

where $A,~B,~C$ and $E$ are unknown constants. To leading order, the curvature is $\unicode[STIX]{x1D705}\approx aq(q-1)r^{q-2}$ , and using $F(0)=A+C=a\,Ca$ , (2.2) can be written as a system of four homogeneous equations for $F$ , for which the determinant is

(2.6) $$\begin{eqnarray}\displaystyle \text{Det}=8q(q-1)^{2}(q-2)(\cos q\unicode[STIX]{x03C0})^{2}(2Ca+\tan q\unicode[STIX]{x03C0}). & & \displaystyle\end{eqnarray}$$

For a non-trivial solution to exist, this determinant must vanish. We are interested in solutions such that, in the limit $Ca\rightarrow 0$ , the profile converges towards a static meniscus, so that $q=2$ . This results in a solution that is regular at the contact line, characterized by finite curvature. Indeed, repeating the above calculation for $q=2$ , in which case

(2.7) $$\begin{eqnarray}\displaystyle F(\unicode[STIX]{x1D719})=A\cos 2\unicode[STIX]{x1D719}+B\sin 2\unicode[STIX]{x1D719}+C+E\unicode[STIX]{x1D719}, & & \displaystyle\end{eqnarray}$$

the pressure becomes logarithmically singular:

(2.8) $$\begin{eqnarray}\displaystyle p=-\frac{4Ca\,a}{\unicode[STIX]{x03C0}}\ln r. & & \displaystyle\end{eqnarray}$$

The logarithmic behaviour of the pressure is reminiscent of the logarithmic pressure behaviour of a closing ‘hinged plate’ and a plate in contact with a constant surface stress (Moffatt Reference Moffatt1964). This contradicts Mahadevan & Pomeau’s (Reference Mahadevan and Pomeau1999) claim that the $q=2$ solution is singularity-free for any $Ca$ . In fact, the normal stress balance is satisfied to leading order only if $Ca=0$ .

Instead, to satisfy all boundary conditions (2.2) we make (2.6) vanish by choosing

(2.9) $$\begin{eqnarray}\displaystyle \tan (q\unicode[STIX]{x03C0})=-2Ca, & & \displaystyle\end{eqnarray}$$

the branches of which are shown in figure 6. The condition that $q=2$ for vanishing $Ca$ singles out the branch shown in red, since we expect $q$ to decrease with increasing $Ca$ , such that the interface curvature increases with increasing speed. Solving for $q$ , we find

(2.10) $$\begin{eqnarray}\displaystyle q=\frac{\arctan (-2Ca)+2\unicode[STIX]{x03C0}}{\unicode[STIX]{x03C0}}\in \left(\frac{3}{2},2\right]. & & \displaystyle\end{eqnarray}$$

Higher branches correspond to subdominant solutions, while lower branches are unphysical. Our conclusions, to be confirmed by comparison to numerical simulation below, are different from those of Mahadevan & Pomeau (Reference Mahadevan and Pomeau1999) and of Benilov & Vynnycky (Reference Benilov and Vynnycky2013), who find $q\geqslant 2$ . For large $Ca$ , the power law converges towards $q=3/2$ , which is the generic cusp (Eggers & Suramlishvili Reference Eggers and Suramlishvili2017) found for the problem without a solid plate (Jeong & Moffatt Reference Jeong and Moffatt1992).

Since the energy dissipation density $\unicode[STIX]{x1D716}$ behaves like the square of a velocity gradient, it is seen from power counting (and confirmed by direct calculation) that $\unicode[STIX]{x1D716}\propto r^{2(q-2)}$ , so the total dissipation is finite for $q>1$ . This confirms that the usual contact line singularity (Huh & Scriven Reference Huh and Scriven1971) is regularized in the region of interest. To compute the velocity field, we use that $A+C=Ca\,a$ , which leads to the streamfunction

(2.11) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713} & = & \displaystyle -Ca\,r\sin \unicode[STIX]{x1D719}\nonumber\\ \displaystyle & & \displaystyle +\,ar^{q}\left[\frac{2-q}{2}Ca\cos q\unicode[STIX]{x1D719}+\frac{2-q}{4}\sin q\unicode[STIX]{x1D719}+\frac{q}{2}Ca\cos (q-2)\unicode[STIX]{x1D719}+\frac{q}{4}\sin (q-2)\unicode[STIX]{x1D719}\right].\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Figure 7.

Plots of FEM profiles of the liquid interface (black curved line) at different $Ca$ . The (red) solid straight lines are best fits to the intermediate region $h=ax^{q}$ , where $q$ is defined by (2.10), and the (red) dashed straight lines are those to the cusp exponent $3/2$ . The prefactor $a$ (in units of $l_{c}$ ) is $0.76$ and $0.49$ , respectively. The best fits are made to approximate the shape of the free surface; the slip length is $\unicode[STIX]{x1D706}=10^{-4}$ .

In figure 7 we compare (2.3) and (2.10) to numerical simulations of the full problem for two different values of $Ca$ and $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ . For each $Ca$ , (2.3) is fitted in a region $1\gg h\gg R$ , where the prefactor $a$ is used as a fitting parameter (see figure 8). The fits are shown as solid straight lines, while corresponding fits using the asymptotic value $q=3/2$ are shown as the dashed straight lines. The quality of the fits increases rapidly with $Ca$ , and for $Ca=2.51$ the fit is over three decades in $x$ . Figure 7 also demonstrates that, while $q$ comes quite close to $3/2$ (which is the value used by Jacqmin (Reference Jacqmin2002)), the value (2.10) given by theory still provides an improved fit. This demonstrates that $a$ is determined as a result of the matching between the solution of Benney & Timson (Reference Benney and Timson1980) and an outer solution. By contrast, Ngan & Dussan V. (Reference Ngan and Dussan V.1984) rejected (2.3) on the grounds that $a$ was not determined as part of a local solution.

Figure 8.

The prefactor $a$ in (2.3) as a function of $Ca$ . The circles represent values of $a$ obtained by fitting (2.3) to numerical simulations as shown in figure 7. The solid line is an empirical fit $a=0.45+0.82\exp (-1.31Ca)$ , suggesting that $a$ will rapidly approach a finite value as $Ca\rightarrow \infty$ .

2.2 Region II: cusp tip and the contact line

We begin by analysing the case $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ . In the case of perfect slip, this would be the same solution as that of no wall, with a line of symmetry at $y=0$ . For that case, Jeong & Moffatt (Reference Jeong and Moffatt1992) have found a local similarity solution of the form

(2.12) $$\begin{eqnarray}\displaystyle h=\sqrt{2Rx}+ax^{3/2}, & & \displaystyle\end{eqnarray}$$

where $a$ is a constant and $R$ is the radius of curvature at the tip. The profile (2.12) is the generic form of the singularity of a smooth curve as it is about to make a self-intersection (Eggers & Suramlishvili Reference Eggers and Suramlishvili2017), so we expect it to be valid generically, independent of the flow geometry. In the case of a finite slip length, we now propose on phenomenological grounds that the local cusp solution is

(2.13) $$\begin{eqnarray}\displaystyle h=\sqrt{2Rx}+ax^{q}, & & \displaystyle\end{eqnarray}$$

where $q$ is the exponent (2.10). Just like the original solution of Jeong & Moffatt (Reference Jeong and Moffatt1992), this can be cast in the similarity form

(2.14a-c ) $$\begin{eqnarray}\displaystyle h=R^{q/(2q-1)}H(\unicode[STIX]{x1D709}),\quad \unicode[STIX]{x1D709}=R^{1/(1-2q)}x,\quad H(\unicode[STIX]{x1D709})\approx \sqrt{2\unicode[STIX]{x1D709}}+a\unicode[STIX]{x1D709}^{q}. & & \displaystyle\end{eqnarray}$$

This brings out the fact that the cusp solution possesses a single characteristic length scale, $R$ .

Figure 9.

The composite solution (2.13) (red dashed line) fitted to a FEM simulation of the free surface (black solid line) for two different contact angles; $\unicode[STIX]{x1D706}=10^{-4}$ . The fitting parameters are $R/\unicode[STIX]{x1D706}\approx 3.4\times 10^{-4}$ and $a=0.52$ (a), and $R/\unicode[STIX]{x1D706}=0.011$ $a=0.66$ (b). The inset on the left shows that the composite solution (2.13) slightly underestimates the rate of increase of the interface near the crossover.

Figure 9(a) shows excellent agreement between our asymptotic description (2.13) and a numerical simulation for $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ . Only when taking the first derivative can a small discrepancy in the crossover region be detected. The adjustable parameters are the amplitude $a$ of the cusp solution and the radius of curvature $R$ . It is clear that for $\unicode[STIX]{x1D703}_{e}\neq \unicode[STIX]{x03C0}/2$ the exponent $1/2$ naturally cannot be valid all the way to the contact line. However, figure 9(b) demonstrates that (2.13) is accurate for a remarkable portion of the interface before eventually failing on a very small scale.

This shows that, contrary to the assumptions underlying the conventional theory of the drying transition, the contact line region becomes all but obliterated. Even at modest $Ca\approx 1$ the crossover to the contact line region only takes place at a scale of $10^{-6}$ , which is two orders of magnitude below the slip length $\unicode[STIX]{x1D706}$ . The smallness of $R$ in the continuum description raises fundamental questions as to how the contact line region should be modelled in a physical description, to which we return in the discussion (§ 4). We were not able to provide a comparison for larger $Ca$ , since we can no longer guarantee that the contact line region is fully resolved. We believe the actual behaviour in the contact line region is not as important as is generally believed, since $R$ is the smallest scale that determines the transition. This is at least qualitatively consistent with the conclusions of Eddi, Winkels & Snoeijer (Reference Eddi, Winkels and Snoeijer2013), who find that wetting effects are unimportant for the early stages of drop spreading. Similarly, in Latka et al. (Reference Latka, Boelens, Nagel and de Pablo2018) it was found that the transition to splashing after drop impact on a solid surface is insensitive to the wetting properties of the substrate. The conventional theory of the drying transition based on the Cox–Voinov formula (Cox Reference Cox1986; Kistler Reference Kistler and Berg1993; Vandre et al. Reference Vandre, Carvalho and Kumar2013) aims to describe how the interface angle increases from $\unicode[STIX]{x1D703}_{e}$ to a value close to $\unicode[STIX]{x03C0}$ , making a convex shape. However, at the turning point, the interface becomes concave, forming the cusp region, described above. This shows that the conventional theory fails to describe the interface in even a qualitative fashion.

2.2.1 The radius of curvature $R$

Returning to the tip region, our main task is to develop a theory for the radius of curvature $R$ . This is important, since $R$ determines the smallest scale of the cusp into which the gas phase is confined. As a result, in analogy to Eggers (Reference Eggers2001), $R$ determines the saddle-node bifurcation if gas is included. In the classical theory of Jeong & Moffatt (Reference Jeong and Moffatt1992) (in the absence of a plate), $R$ depends exponentially on the capillary number. In that paper, the authors report an argument by Hinch which represents the cusp tip by a point force of strength $2\unicode[STIX]{x1D6FE}$ , which comes from the vertical upward pull of each side of the cusp. At the tip of the cusp, the upward velocity generated by the point force (often called a stokeslet) must cancel the downward velocity along the cusp walls. Since the speed generated by a stokeslet depends logarithmically on the distance in two dimensions, this leads to an exponential dependence of the tip scale on the imposed velocity field.

In order to adapt this idea to the present situation, we have to replace the free-space stokeslet with its analogue in the presence of a wall. In order for the tip of the cusp to be stationary, the speed generated at the contact line by a force $\boldsymbol{F}=\unicode[STIX]{x1D6FE}\boldsymbol{e}_{x}$ at a distance $r$ above the contact line must equal $U$ . This is the $(x,x)$ component of the Stokes Green function in the presence of a partial-slip wall, assuming that the force is situated on the wall $y=0$ :

(2.15) $$\begin{eqnarray}\displaystyle U=\frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D702}}G_{xx}^{PS}(r,\unicode[STIX]{x1D706}). & & \displaystyle\end{eqnarray}$$

In analogy with the three-dimensional case (Lauga & Squires Reference Lauga and Squires2005), the Green function can be written as the sum of the free-space Green function, its image and a correction factor $W_{xx}^{PS}(r,\unicode[STIX]{x1D706})$ . Since the force is on the wall, the image is the same as the Green function itself, and we obtain

(2.16) $$\begin{eqnarray}\displaystyle G_{xx}^{PS}(r,\unicode[STIX]{x1D706}) & = & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}}(1-\ln r)+\int _{0}^{\infty }\text{e}^{-u/2}(\ln \hat{r}-1)\,\text{d}u\nonumber\\ \displaystyle & & \displaystyle -\,4\unicode[STIX]{x1D706}^{2}\int _{0}^{\infty }\left(1-\left(1+\frac{1}{2}\right)\text{e}^{-u/2}\right)\frac{1}{\hat{r}^{2}}\,\text{d}u,\quad \hat{r}=\sqrt{x^{2}+(\unicode[STIX]{x1D706}u)^{2}}.\end{eqnarray}$$

The slip length $\unicode[STIX]{x1D706}$ enters as a weighted integral of the no-slip Green function with respect to the slip length.

We have seen that the scale $R$ in fact becomes small compared to $\unicode[STIX]{x1D706}$ , hence we analyse (2.16) in the limit of small $r$ , which results in

(2.17) $$\begin{eqnarray}\displaystyle U=\frac{\unicode[STIX]{x1D6FE}}{2\unicode[STIX]{x03C0}\unicode[STIX]{x1D702}}[-\text{ln}\,r/\unicode[STIX]{x1D706}+1+(\unicode[STIX]{x1D6FE}_{E}+1)]+O(r/\unicode[STIX]{x1D706}), & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}_{E}$ is Euler’s constant.

Assuming that the radius of curvature $R=Cr$ of the tip scales like $r$ , we obtain

(2.18) $$\begin{eqnarray}\displaystyle R=\text{e}^{2+\unicode[STIX]{x1D6FE}_{E}}C\unicode[STIX]{x1D706}\exp (-2\unicode[STIX]{x03C0}Ca), & & \displaystyle\end{eqnarray}$$

where $C$ is an empirical constant. As a result of the point force now being at a distance $R$ from the fluid, the $r^{-1}$ stress singularity encountered by Huh & Scriven (Reference Huh and Scriven1971) is now converted into an integrable $r^{-1/2}$ singularity (Jeong & Moffatt Reference Jeong and Moffatt1992; Moffatt Reference Moffatt1993). One can also check that the scaling of (2.18) with $Ca$ is what is needed to make the local dissipation finite, independent of the fluid viscosity. To calculate its exact numerical value of $R$ , a complete theory of the flow in the tip region would be necessary. In figure 10 we find excellent agreement with (2.18), fitting the constant to find $C=4.29$ . Finally, figure 11 shows the dependence for different $\unicode[STIX]{x1D703}_{e}$ . We find that (2.18) is still valid, at least for high enough $Ca$ . This confirms our earlier conclusion that the form of the cusp region is independent of $\unicode[STIX]{x1D703}_{e}$ ; however, the value of $C$ depends significantly on  $\unicode[STIX]{x1D703}_{e}$ .

Figure 10.

Plot of $R/\unicode[STIX]{x1D706}$ for $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ and different slip lengths. The radius $R$ is found by fitting (2.13) to FEM simulations. All data collapse according to (2.18), with $C=4.29$ .

Figure 11.

Plot of $R/\unicode[STIX]{x1D706}$ for different $\unicode[STIX]{x1D703}_{e}$ . The radius $R$ is found by fitting (2.13) to FEM simulations, just like for $\unicode[STIX]{x1D703}_{e}=90^{\circ }$ . As $Ca$ increases, all curves eventually collapse onto a line with exponent $-2\unicode[STIX]{x03C0}$ . The value of $C$ in (2.18) is $C=3.22$ for $\unicode[STIX]{x1D703}_{e}=80^{\circ }$ , $C=4.29$ for $\unicode[STIX]{x1D703}_{e}=90^{\circ }$ and $C=146$ for $\unicode[STIX]{x1D703}_{e}=100^{\circ }$ .

2.3 Region III: the bath

To complete the description of the profile for $M=0$ , we consider how the profile levels off towards a flat bath. The flow far from the interface should be well described by a flow in a rectangular corner, driven by the moving vertical plate. The other side of the corner is formed by the free surface. This amounts to using the GL approximation (see appendix A) in the limit $\unicode[STIX]{x1D703}\rightarrow \unicode[STIX]{x03C0}/2$ , for which $F\approx -4/(3\unicode[STIX]{x03C0})$ , as found from (A 2). Surface tension is subdominant on a scale much larger than $l_{c}$ , as we will confirm self-consistently; of course, we can set $\unicode[STIX]{x1D706}=0$ as well. Thus from (A 1) we obtain to leading order in $\unicode[STIX]{x1D703}-\unicode[STIX]{x03C0}/2$ in dimensional variables

(2.19) $$\begin{eqnarray}\displaystyle \frac{4U\unicode[STIX]{x1D708}}{\unicode[STIX]{x03C0}gh^{2}}=\unicode[STIX]{x1D703}-\unicode[STIX]{x03C0}/2, & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D708}=\unicode[STIX]{x1D702}/\unicode[STIX]{x1D70C}$ . From (A 4) we obtain to leading order

(2.20) $$\begin{eqnarray}\displaystyle \frac{\text{d}h}{\text{d}x}=\frac{1}{\unicode[STIX]{x1D703}-\unicode[STIX]{x03C0}/2}, & & \displaystyle\end{eqnarray}$$

so that

(2.21) $$\begin{eqnarray}\displaystyle \frac{4U\unicode[STIX]{x1D708}}{\unicode[STIX]{x03C0}g}\frac{\text{d}h}{h^{2}}=\text{d}x. & & \displaystyle\end{eqnarray}$$

Integrating, we find

(2.22) $$\begin{eqnarray}\displaystyle h=a_{1}\frac{1}{\unicode[STIX]{x1D6E5}-x},\quad a_{1}=\frac{4U\unicode[STIX]{x1D708}}{\unicode[STIX]{x03C0}g}. & & \displaystyle\end{eqnarray}$$

This agrees well with numerical simulations for $H_{b}=10^{3}$ (dimensionally a domain width $1000$ times the capillary length); as shown in the inset in figure 12, data points approach the theoretical prediction for large  $Ca$ .

Figure 12.

FEM simulations of the interface about the bath for $Ca=0.75$ , 1.04 and 1.49 (bottom to top) for $H_{b}=10^{3}$ (black lines). The short (red) line in the main panel represents the power law $h\sim (\unicode[STIX]{x1D6E5}-x)^{-1}$ . In the inset we show the coefficient of proportionality $a_{1}$ as found by fitting FEM simulations to $h=a_{1}(\unicode[STIX]{x1D6E5}-x)^{-1}$ and compare to the theoretical prediction (2.22) (red line).

3.1 The bifurcation

We solve the gas flow in the narrow space between the liquid and solid for a given $M=0$ solution. Since the geometry is slender, we can use lubrication theory (Eggers & Fontelos Reference Eggers and Fontelos2015) to calculate the pressure inside the gap. As usual in lubrication theory, the pressure is constant over a cross-section of the gap, and the normal stress the gas exerts on the fluid is dominated by the pressure. The shear stress is subdominant in lubrication theory. If $u_{0}(x)$ is the vertical velocity of the fluid on the interface, the vertical velocity in the gap, using a quadratic approximation, is

(3.1) $$\begin{eqnarray}\displaystyle u_{g}(x,y)=a_{2}y^{2}+a_{1}y+a_{1}\unicode[STIX]{x1D706}_{g}-U,\quad a_{2}=\frac{p_{x}}{2M\unicode[STIX]{x1D702}},\quad a_{1}=\frac{u_{0}+U-a_{2}h^{2}}{h+\unicode[STIX]{x1D706}_{g}}. & & \displaystyle\end{eqnarray}$$

It is easily verified that, for $y=h(x)$ , the velocity is continuous, while for $y=0$ (on the plate), the partial-slip condition (1.6) (third equation) is satisfied. From the condition that the total flux $\int _{0}^{h}u_{g}\,\text{d}y$ through the gap must vanish, we obtain

(3.2) $$\begin{eqnarray}\displaystyle \frac{\text{d}p^{lb}}{\text{d}x}=6M\unicode[STIX]{x1D702}\frac{(u_{0}-U)h+2u_{0}\unicode[STIX]{x1D706}_{g}}{h^{2}(h+4\unicode[STIX]{x1D706}_{g})} & & \displaystyle\end{eqnarray}$$

for the derivative of the lubrication pressure with respect to $x$ . To obtain the pressure, one can integrate (3.2) starting from the bath, where the pressure is that of the ambient gas. The lubrication pressure increases rapidly as the gap becomes narrower, which is the root cause of the bifurcation. The growth of the pressure is mitigated somewhat by the presence of the slip $\unicode[STIX]{x1D706}_{g}$ . The liquid flow is calculated for $M=0$ , which supplies $u_{0}$ , from which the lubrication pressure is calculated via (3.2). This produces a correction to the stress balance on the free surface, which changes the liquid flow, so both liquid and gas flow have to be calculated self-consistently. A scheme of calculating the effect of the gas by lubrication theory was implemented numerically by Liu et al. (Reference Liu, Vandre, Carvalho and Kumar2016) and Sprittles (Reference Sprittles2017). In both papers it is confirmed that, at least when the capillary number is sufficiently high, results based on the lubrication approximation in the gas are almost indistinguishable from numerical simulations where both phases are treated equally.

However, this approach still requires a full numerical simulation of the liquid phase for each value of $M$ . In order to capture the effect of the gas analytically, as in Eggers (Reference Eggers2001) we observe that the cusp is similar to a crack in an infinite two-dimensional fluid domain, with a normal load imposed on it. Exploiting the equivalence between linear elasticity and the Stokes equation in two dimensions, one can use classical results from elasticity (Sun Reference Sun2011) to compute the extra velocity $v_{M}(x)$ coming from the stress on the cusp surface:

(3.3) $$\begin{eqnarray}\displaystyle v_{M}(x)=\frac{1}{\unicode[STIX]{x1D702}}\int _{0}^{\unicode[STIX]{x1D6E5}}p^{lb}(x^{\prime })m(x^{\prime }/x)\,\text{d}x^{\prime }, & & \displaystyle\end{eqnarray}$$

where

(3.4) $$\begin{eqnarray}\displaystyle m(x)=\frac{1}{2\unicode[STIX]{x03C0}}\ln \left|\frac{1+\sqrt{x}}{1-\sqrt{x}}\right|. & & \displaystyle\end{eqnarray}$$

The upper limit $\unicode[STIX]{x1D6E5}$ represents the undisturbed bath, where the ambient pressure has been reached. Note the crucial feature that $v_{M}$ is a non-local function of the load $p^{lb}$ , while in the GL approximation, where (3.2) enters into (A 1) through $F(\unicode[STIX]{x1D703},M)$ , the interface description is affected by $p^{lb}$ only locally.

Integrating (3.3) by parts gives

(3.5) $$\begin{eqnarray}\displaystyle v_{M}(x)=-\frac{x}{\unicode[STIX]{x1D702}}\int _{0}^{\unicode[STIX]{x1D6E5}}p_{x}^{lb}(x^{\prime }){\mathcal{M}}(x^{\prime }/x)\,\text{d}x^{\prime }, & & \displaystyle\end{eqnarray}$$

where

(3.6) $$\begin{eqnarray}\displaystyle {\mathcal{M}}(x)=\int _{0}^{x}m(x^{\prime })\,\text{d}x^{\prime }=\frac{1}{\unicode[STIX]{x03C0}}\left[\sqrt{x}+\frac{x-1}{2}\ln \left|\frac{1+\sqrt{x}}{1-\sqrt{x}}\right|\right]. & & \displaystyle\end{eqnarray}$$

The requirement that the free surface be a streamline of the flow is

(3.7) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}h}{\unicode[STIX]{x2202}x}=\frac{v(x)}{u(x)}\approx \frac{v_{0}(x)+v_{M}(x)}{u_{0}(x)}=\frac{\unicode[STIX]{x2202}h_{0}}{\unicode[STIX]{x2202}x}+\frac{\unicode[STIX]{x2202}h^{c}}{\unicode[STIX]{x2202}x},\quad h=h_{0}+h^{c}(x), & & \displaystyle\end{eqnarray}$$

where we have assumed that $v_{M}$ can be added to the base flow in a perturbative fashion, as explained before. Here $h_{0}(x)$ represents the initial $M=0$ base profile, and $h^{c}(x)$ the perturbation from this base profile for a given $M$ . Since $v_{0}$ and $v_{M}$ have opposite signs, the effect of the air is to push the interface so as to make it steeper, effectively narrowing the gap, decreasing $h$ . According to (3.2), the pressure rapidly increases with decreasing $h$ , amplifying the effect of the air. This positive feedback leads to a bifurcation at a critical value of the capillary number.

As described in the introduction and seen in figure 3, to find both stable and unstable branches of the bifurcation curve, one fixes the depression $\unicode[STIX]{x1D6E5}$ , and searches for the corresponding value of $Ca$ . However, since the base profile, which is available to us numerically only, depends on $Ca$ but is by definition independent of $M$ , it is more convenient to fix $\unicode[STIX]{x1D6E5}$ as well as $Ca$ and search for $M$ , as seen in figure 13. Once the critical value of $M$ has been found for different values of $Ca$ , one can produce the conventional plot of $Ca_{cr}$ as a function of $M$ (cf. figure 14).

Equations (3.5) and (3.7) are solved numerically for $h^{c}(x)$ . The unperturbed profile $h_{0}$ , as well as the base solution $u_{0}(x)$ and $v_{0}(x)$ for the velocity field, evaluated at the free surface, are taken from the numerical simulation for $M=0$ . Starting from $\unicode[STIX]{x1D6E5}(M=0)$ , the depression is increased slowly, and the value of $M$ is sought, using the previous solution as an initial condition. At each new value of $\unicode[STIX]{x1D6E5}$ , the profile $h(x)$ is discretized, and (3.7) is solved using Newton’s method. Since the changes in both the profile and $M$ from the preceding step are very small, only a few iterations are needed for both the new profile and value of $M$ to converge to high accuracy. The process is repeated and the value of $M(\unicode[STIX]{x0394})$ recorded. Successive iterations result in a bifurcation curve as shown in figure 13 as the solid black line, for $Ca=0.65$ and $Ca=1.04$ . Thus, to compare how the perturbation $\unicode[STIX]{x0394}(M)$ and $h_{M}^{c}(x)$ increases compared to the numerical simulations, we first compare bifurcation plots perturbing from the same $M=0$ solution, found from the FEM simulations. This is shown in figure 13.

Figure 13.

Bifurcation curve obtained from theory, using the base flow obtained from numerical simulation at $M=0$ for $Ca=0.65$ and $Ca=1.04$ (solid black line), compared to full FEM simulations (red line). The contact angle is $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ and slip lengths are $\unicode[STIX]{x1D706}_{l}=\unicode[STIX]{x1D706}_{g}=10^{-4}$ . Vertical dashed lines mark the critical values $M_{cr}$ .

Figure 14.

The critical capillary number $Ca_{cr}$ as a function of $M$ according to FEM data (red squares) and theory (blue circles), for $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{g}=10^{-4}$ and $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ . The solid line corresponds to the analytical prediction according to (3.14), and the dashed line is the analytical prediction including slip between gas and the liquid; the vertical dotted line denotes the theoretical limit (3.19) below which no bifurcation occurs.

The plots for $M(\unicode[STIX]{x1D6E5})$ in figure 13 show a saddle-node bifurcation: there is a critical $M_{cr}$ at which the upper branch from this point is stable and the lower branch is unstable. This bifurcation point corresponds to the point where a dynamic drying transition would occur as the value of $M$ is raised. The lower branch is unstable. The red line corresponds to the result of the FEM simulation; it was obtained by extrapolating the bifurcation curve $Ca(\unicode[STIX]{x1D6E5})$ obtained numerically for a range of $M$ values to a curve $M(\unicode[STIX]{x1D6E5})$ at fixed $Ca$ . It is seen that, for the larger value of $Ca$  ( $Ca=1.04$ ), the agreement is almost perfect along the upper branch, and the extrapolated value of $M$ where the bifurcation occurs is very close to the bifurcation point predicted by theory. However, owing to computational issues, previously described, we were not able to capture the lower branch in numerical simulations. We estimate the bifurcation point from a sudden increase of the curvature of the bifurcation curve. Trying to extrapolate our data beyond the bifurcation point, we estimate that our critical values of $M$ might be too low by as much as 10 %. For the smaller $Ca$  ( $Ca=0.65$ ), a lower branch is seen in both the numerical simulation and theory. The price we have to pay is that our asymptotic description is not quite as good, and the value of the critical $M$ is overestimated by approximately $20\,\%$ . Still, the bifurcation curve is predicted convincingly, without adjustable parameters. The comparison between FEM results and theory is summarized in figure 14 for an angle of $\unicode[STIX]{x1D703}_{e}=\unicode[STIX]{x03C0}/2$ and slip lengths $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}_{g}=10^{-4}$ . The FEM data are the same as given in figure 2, with slip lengths chosen to match the simulations of Vandre et al. (Reference Vandre, Carvalho and Kumar2013). As shown in Sprittles (Reference Sprittles2017), the experimental results of figure 2 are well reproduced by numerical simulations, if appropriate physical slip lengths are chosen. Good agreement is found between our theory and simulations in figure 14 as long as $M$ values are sufficiently small for our asymptotic theory to be applicable. For small $M$ one observes a small rise of the critical capillary number over a generally logarithmic behaviour. This is because the presence of slip reduces the amount of air being dragged into the cusp region. To gain a better analytical understanding of these behaviours, we now present a simplified analytical theory of the transition.

3.2 Similarity description

We formulate (3.5) and (3.7) in terms of the phenomenological cusp solution (2.13). In order to be able to write the resulting equation in self-similar variables, we assume that the velocity $u_{0}=-U$ along the cusp is constant. If (2.14) represents the base solution for $M=0$ , the full solution reads in similarity variables:

(3.8a-c ) $$\begin{eqnarray}\displaystyle h=R^{q/(2q-1)}H(\unicode[STIX]{x1D709}),\quad \unicode[STIX]{x1D709}=R^{1/(1-2q)}x,\quad H(\unicode[STIX]{x1D709})=H_{0}(\unicode[STIX]{x1D709})+H^{c}(\unicode[STIX]{x1D709}), & & \displaystyle\end{eqnarray}$$

where $H_{0}(\unicode[STIX]{x1D709})$ is the unperturbed ( $M=0$ ) similarity solution described by (2.14), and $H^{c}(\unicode[STIX]{x1D709})$ is the correction coming from the effect of air in the narrow gap. Using $u_{0}=-U$ , we have

(3.9) $$\begin{eqnarray}\displaystyle \frac{\text{d}p^{lb}}{\text{d}x}=-\frac{12M\unicode[STIX]{x1D702}U(h+\unicode[STIX]{x1D706}_{g})}{h^{2}(h+4\unicode[STIX]{x1D706}_{g})}, & & \displaystyle\end{eqnarray}$$

so using (3.7) and (3.6), the equation for $H^{c}(\unicode[STIX]{x1D709})$ becomes

(3.10) $$\begin{eqnarray}\displaystyle \frac{\text{d}H^{c}}{\text{d}\unicode[STIX]{x1D709}}=-12s\unicode[STIX]{x1D709}\int _{0}^{\infty }\frac{{\mathcal{M}}(\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D709})(H(\unicode[STIX]{x1D70F})+\bar{\unicode[STIX]{x1D706}}_{g})\,\text{d}\unicode[STIX]{x1D70F}}{H^{2}(\unicode[STIX]{x1D70F})(H(\unicode[STIX]{x1D70F})+4\bar{\unicode[STIX]{x1D706}}_{g})}, & & \displaystyle\end{eqnarray}$$

where

(3.11a,b ) $$\begin{eqnarray}\displaystyle s=MR^{3(q-1)/(1-2q)},\quad \bar{\unicode[STIX]{x1D706}}_{g}=R^{q/(1-2q)}\unicode[STIX]{x1D706}_{g}. & & \displaystyle\end{eqnarray}$$

Together with (3.8), this is an equation for the perturbation $H^{c}(\unicode[STIX]{x1D709})$ . From the behaviour of the kernel,

(3.12a,b ) $$\begin{eqnarray}\displaystyle \lim _{x\rightarrow 0}{\mathcal{M}}(x)=\frac{2x^{3/2}}{3\unicode[STIX]{x03C0}},\quad \lim _{x\rightarrow \infty }{\mathcal{M}}(x)=\frac{2\sqrt{x}}{\unicode[STIX]{x03C0}}, & & \displaystyle\end{eqnarray}$$

one can deduce that $(\text{d}H^{c}/\text{d}\unicode[STIX]{x1D709})\unicode[STIX]{x1D709}^{1/2}$ for $\unicode[STIX]{x1D709}\rightarrow 0$ and $\text{d}H^{c}/\text{d}\unicode[STIX]{x1D709}\propto \,\unicode[STIX]{x1D709}^{-1/2}$ for $\unicode[STIX]{x1D709}\rightarrow \infty$ . It follows that $H_{c}(\unicode[STIX]{x1D709})$ behaves asymptotically as

(3.13a,b ) $$\begin{eqnarray}\displaystyle \lim _{\unicode[STIX]{x1D709}\rightarrow 0}H^{c}(\unicode[STIX]{x1D709})=-H_{0}\unicode[STIX]{x1D709}^{3/2},\quad \lim _{\unicode[STIX]{x1D709}\rightarrow \infty }H^{c}(\unicode[STIX]{x1D709})=-H_{i}\unicode[STIX]{x1D709}^{1/2}, & & \displaystyle\end{eqnarray}$$

where $H_{0}$ and $H_{i}$ are constants, which need to be found as part of the solution of (3.10). In particular, the integral in (3.10) is convergent at both the lower and upper boundary, so non-universal effects having to do with either the bath or the contact line region do not play a role asymptotically.

Now we solve (3.10) numerically for a given capillary number, putting $H(\unicode[STIX]{x1D709})=\sqrt{2\unicode[STIX]{x1D709}}+a\unicode[STIX]{x1D709}^{q}+H^{c}(\unicode[STIX]{x1D709})$ . This is essentially the same equation as that solved numerically in Eggers (Reference Eggers2001) using Newton’s method. We find a saddle-node bifurcation at a critical value $s=s_{c}$ , above which there is no more solution, while below $s_{c}$ there are two branches, one stable and one unstable. To identify both branches, we treat $s$ as an unknown in (3.10), holding $H_{i}$ (cf. (3.13)) constant. Clearly, the larger values of $H_{i}$ correspond to the unstable branch, the smaller values to the stable branch; solving (3.10), we find $s$ for each value of $H_{i}$ . The maximum of the curve $s(H_{i})$ corresponds to  $s_{c}$ .

With $s_{c}(q,a,\bar{\unicode[STIX]{x1D706}}_{g})$ in hand, we are able to make a theoretical prediction for the critical capillary number shown in figure 14, without using a FEM simulation of the $M=0$ solution. The dependence of the radius of curvature on $Ca$ is given by (2.18), and so the first equation of (3.11) yields

(3.14) $$\begin{eqnarray}\displaystyle M=s_{c}(q,a,\bar{\unicode[STIX]{x1D706}}_{g})(\text{e}^{2+\unicode[STIX]{x1D6FE}_{E}}C\unicode[STIX]{x1D706})^{3(1-q)/(1-2q)}\exp \left(-\frac{6\unicode[STIX]{x03C0}(q-1)}{2q-1}Ca_{cr}\right). & & \displaystyle\end{eqnarray}$$

This is the inverse of $Ca_{cr}(M)$ , and is shown as the solid line in figure 14. We have used $C=4.29$ (cf. figure 10), $q$ as calculated from (2.10), and $a$ using the fit from figure 8. The rescaled slip length $\bar{\unicode[STIX]{x1D706}}_{g}$ in the gas is found from (3.11). The expression (3.14) combines all theoretical results from this paper, providing a unified asymptotic description of the critical capillary number in terms of all relevant parameters.

The solid line agrees well with both full FEM simulations (red squares) and the abridged theory based on knowledge of the full solution for $M=0$ (blue circles), as long as $Ca_{cr}$ is sufficiently high, meaning that $M$ is smaller than approximately $10^{-4}$ . For moderate values of $M$ , meaning that $Ca$ is small, the rescaled slip length $\bar{\unicode[STIX]{x1D706}}_{g}$ is small and can effectively be put to zero. Thus, apart from the weak dependence of $q$ and $a$ on $Ca$ , $s_{c}$ is a constant, and solving (3.14) for $Ca_{cr}$ leads to the logarithmic behaviour

(3.15) $$\begin{eqnarray}\displaystyle Ca_{cr}=\frac{2q-1}{6\unicode[STIX]{x03C0}(1-q)}\ln \frac{M}{s_{c}}+\frac{2+\unicode[STIX]{x1D6FE}_{E}+\ln (C\unicode[STIX]{x1D706})}{2\unicode[STIX]{x03C0}}, & & \displaystyle\end{eqnarray}$$

which is seen in figure 14 for $M\lesssim 10^{-4}$ . For smaller values of $M$ , on the other hand, $Ca_{cr}$ becomes somewhat larger than predicted by this logarithmic law. The reason is that slip between the gas and the solid wall regularizes the growth of the lubrication pressure, so higher speeds are needed to trigger the bifurcation. However, as seen from (3.10), even in the limit of $\bar{\unicode[STIX]{x1D706}}_{g}\rightarrow \infty$ , the gradient is reduced by a factor of $1/2$ only, compared to $\bar{\unicode[STIX]{x1D706}}_{g}=0$ .

The reason for the insensitivity to $\unicode[STIX]{x1D706}_{g}$ is that we still do not allow slip between the gas and the liquid, so that the gas is still dragged into the gap by the liquid’s motion, with the channel effectively twice as wide, as the gas does not encounter any resistance at the wall. If one were to treat the gas flow near the interface with the same slip law as with the wall, as suggested by kinetic effects (Li Reference Li2016; Sprittles Reference Sprittles2017), (3.10) would turn into

(3.16) $$\begin{eqnarray}\displaystyle \frac{\text{d}H^{c}}{\text{d}\unicode[STIX]{x1D709}}=-12s\unicode[STIX]{x1D709}\int _{0}^{\infty }\frac{{\mathcal{M}}(\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D70F}}{H(\unicode[STIX]{x1D70F})(H(\unicode[STIX]{x1D70F})+6\bar{\unicode[STIX]{x1D706}}_{g})}. & & \displaystyle\end{eqnarray}$$

Clearly, the effect of the gas on the right-hand side now becomes small for large $\bar{\unicode[STIX]{x1D706}}$ , i.e. for small $R$ . Instead of the solid line in figure 14, we now obtain the dashed line, which rises rapidly at $M\approx 10^{-6}$ , leading to a sharp deviation from the weak logarithmic dependence usually associated with the drying transition.

To make this observation more quantitative, we consider the asymptotic limit of $\bar{\unicode[STIX]{x1D706}}_{g}$ very large in (3.16), so that $H(\unicode[STIX]{x1D70F})$ can be neglected in comparison. As a result, we obtain

(3.17) $$\begin{eqnarray}\displaystyle \frac{\text{d}H^{c}}{\text{d}\unicode[STIX]{x1D709}}=-2\bar{s}\unicode[STIX]{x1D709}\int _{0}^{\infty }\frac{{\mathcal{M}}(\unicode[STIX]{x1D70F}/\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D70F}}{H(\unicode[STIX]{x1D70F})}, & & \displaystyle\end{eqnarray}$$

where we have introduced $\bar{s}=s/\bar{\unicode[STIX]{x1D706}}_{g}$ . For large $Ca$ , the cusp exponent (2.10) becomes $q=3/2$ , so the only remaining parameter in (3.17) is $a$ ; we take it to be its asymptotic value $a=0.45$ . Solving the integral equation (3.17) as before, we obtain a critical value $\bar{s}_{c}=0.0272$ above which no solution exists. Combining the definition of $\bar{s}_{c}$ with (3.11), we obtain

(3.18) $$\begin{eqnarray}\displaystyle M=\bar{s}_{c}\unicode[STIX]{x1D706}_{g}R^{(3-2q)/(1-2q)}. & & \displaystyle\end{eqnarray}$$

Note that, for the asymptotic value $q=3/2$ , the exponent vanishes, so we have to include the next order $q=3/2+1/(2\unicode[STIX]{x03C0}Ca)$ to obtain $1/(2\unicode[STIX]{x03C0}Ca)$ for the exponent in (3.18) to leading order as $Ca\rightarrow \infty$ . Together with the formula (2.18) for $R$ , this yields, to leading order as $Ca\rightarrow \infty$ , that the right-hand side of (3.18) reaches a finite value in the limit. This means there exists a critical value

(3.19) $$\begin{eqnarray}\displaystyle M_{c}=0.0272\unicode[STIX]{x1D706}_{g}/e & & \displaystyle\end{eqnarray}$$

of $M$ below which no bifurcation occurs, regardless of how large the capillary number is. For the parameters of figure 14, this is $\log _{10}M_{c}=-6$ (shown as the dotted vertical line), in good agreement with the sharp rise of the dashed line at that point.