Transition in a numerical model of contact line dynamics and forced dewetting

TL;DR

This paper models contact line dynamics and forced dewetting transitions using a numerical approach that aligns with Cox's theory, revealing how contact angle and viscosity influence the critical capillary number.

Contribution

It introduces a numerical method with a discretized contact angle that accurately captures contact line behavior and generalizes the Landau-Levich-Derjaguin transition prediction for various angles and viscosity ratios.

Findings

Cox's theory applies to the numerical model with a gauge function.

Agreement between numerics and theory depends on contact angle and viscosity ratio.

The critical capillary number is generalized for large contact angles and arbitrary viscosity ratios.

Abstract

We investigate the transition to a Landau-Levich-Derjaguin film in forced dewetting using a quadtree adaptive solution to the Navier-Stokes equations with surface tension. We use a discretization of the capillary forces near the receding contact line that yields an equilibrium for a specified contact angle $\theta_\Delta$ called the numerical contact angle. Despite the well-known contact line singularity, dynamic simulations can proceed without any explicit additional numerical procedure. We investigate angles from $15^\circ$ to $110^\circ$ and capillary numbers from $0.00085$ to $0.2$ where the mesh size $\Delta$ is varied in the range of $0.0035$ to $0.06$ of the capillary length $l_c$. To interpret the results, we use Cox's theory which involves a microscopic distance $r_m$ and a microscopic angle $\theta_e$. In the numerical case, the equivalent of $\theta_e$ is the angle $\theta_\Delta$ and we find that Cox's theory also applies. We introduce the scaling factor or gauge function $\phi$ so that $r_m = \Delta/\phi$ and estimate this gauge function by comparing our numerics to Cox's theory. The comparison provides a direct assessment of the agreement of the numerics with Cox's theory and reveals a critical feature of the numerical treatment of contact line dynamics: agreement is poor at small angles while it is better at large angles. This scaling factor is shown to depend only on $\theta_\Delta$ and the viscosity ratio $q$. In the case of small $\theta_e$, we use the prediction by Eggers [Phys. Rev. Lett., vol. 93, pp 094502, 2004] of the critical capillary number for the Landau-Levich-Derjaguin forced dewetting transition. We generalize this prediction to large $\theta_e$ and arbitrary $q$ and express the critical capillary number as a function of $\theta_e$ and $r_m$. An analogy can be drawn between $r_m$ and the numerical slip length.

Tables

Table 1. Table 1: Summary of the simulation Setups. When no value is indicated, the corresponding number is computed from the other numbers using Eq. ( 1 ).

Setup	$N_G={\rho }_1^2{V}_s^3/\left({\mu }_1\left({\rho }_1-{\rho }_2\right)g\right)$	$\mathrm{Re}={\rho }_1{V}_s{l}_c/{\mu }_1$	$1/q={\mu }_1/{\mu }_2$	${\rho }_1/{\rho }_2$
A	$25/64$	-	1	5
B	-	1	1	5
C	-	1	50	5

Equations

$$\text{mbox}Ca={\mu }_1{V}_s/\sigma ,$$

$$\mathrm{Mo}=\frac{({\rho }_1-{\rho }_2)g{\mu }_1^4}{{\rho }_1^2{\sigma }^3}.$$

$$N_G=\frac{{\rho }_1^2{V}_s^3}{{\mu }_1({\rho }_1-{\rho }_2)g}.$$

$$\mathrm{Re}=N_G^{1/2}{\mathrm{Ca}}^{-1/2}.$$

$$\mathrm{Re}=\frac{5}{8}{\mathrm{Ca}}^{-1/2}.$$

$$\rho \left({\partial }_t\text{u}+\text{u}\cdot \nabla \text{u}\right)=-\nabla p+\nabla \cdot \left[\mu \left(\nabla \text{u}+\nabla {\text{u}}^{\top }\right)\right]+\sigma \kappa {\delta }_s\text{n}+\rho \text{g},$$

$$\nabla \cdot \text{u}=0,$$

$${\partial }_t\rho +\text{u}\cdot \nabla \rho =0,$$

$${\partial }_t\chi +\text{u}\cdot \nabla \chi =0,$$

$$n_{x}x+n_{y}y=\alpha ,$$

$$h_i=\sum _{j=j_1}^{j_2}{C}_{i,j}$$

$${\kappa }_{\mathrm{III}}=l_c^{-1}\sqrt{2-2\sin {\theta }_{\mathrm{a}}}$$

$$G[\theta (r)]=G({\theta }_e)-\mathrm{Ca}\ln (r/{\lambda }_c)-\mathrm{Ca}\frac{Q_i}{f({\theta }_e,q)}+o(\mathrm{Ca}),$$

$$G(\theta )={\int }_0^{\theta }\frac{d{\theta }^{\prime }}{f({\theta }^{\prime },q)},$$

$$f(\theta ,q)=\frac{2\sin \theta \{q^2({\theta }^2-{\sin }^2\theta )+2q[\theta (\pi -\theta )+{\sin }^2\theta ]+[(\pi -\theta {)}^2-{\sin }^2\theta ]\}}{q({\theta }^2-{\sin }^2\theta )[(\pi -\theta )+\cos \theta \sin \theta ]+[(\pi -\theta {)}^2-{\sin }^2\theta ](\theta -\cos \theta \sin \theta )}.$$

$$G(\theta )\approx {\theta }^3/9.$$

$$r_m={\lambda }_c\exp [f({\theta }_e,q)/Q_i],,$$

$$G(\theta )=G({\theta }_e)-\mathrm{Ca}\ln (r/r_m),$$

$${v|}_{x=0}-V_s=\lambda {\frac{\partial v}{\partial x}|}_{x=0},$$

$${\eta }^{{}^{\prime \prime \prime }}=\frac{3\mathrm{Ca}}{{\eta }^2+3\lambda \eta },$$

$${\eta }^{\prime }(\zeta )\sim {\theta }_e-\frac{3\mathrm{Ca}}{{\theta }_e^2}\ln \left(\frac{\zeta }{r_m}\right),$$

$$r_m=\frac{3\lambda }{\mathrm{e}{\theta }_e}+\mathcal{O}(\mathrm{Ca}).$$

$${\eta }^{\prime }(\zeta {)}^3\sim {\theta }_e^3-9\mathrm{Ca}\ln (\zeta /r_m).$$

$$r_m=\lambda /\varphi ,$$

$$\varphi ({\theta }_e)=\frac{\mathrm{e}{\theta }_e}{3}.$$

$$r_m=\Delta /{\varphi }_{\mathrm{num}},$$

$$\lambda =\varphi \Delta /{\varphi }_{\mathrm{num}},$$

$$G({\theta }_e)-G[\theta (r)]=\mathrm{Ca}\ln (r/\Delta )+\mathrm{Ca}\ln \varphi .$$

$${\mathrm{Ca}}_{cr}=\frac{G({\theta }_e)}{\ln ({\zeta }_I/r_m)},$$

$$ϵ=1/\ln ({\zeta }_I/r_m)={\mathrm{Ca}}_{\mathrm{cr}}/G({\theta }_e).$$

Full text

††journal: Journal of Computational Physics

1 Introduction

Wetting of solids by liquids, in which a liquid displaces another fluid on a solid substrate, is an ubiquitous phenomenon with applications ranging from coating [2] and tear films on the cornea [3] to micro-layer formation in wall boiling [4, 5] and CO2 sequestration [6]. However, despite the abundance of applications, the precise mechanism of wetting is only partially understood. From the numerical modeling point of view, difficulties arise due to the highly multiscale nature of the problem (length scales extending from the macroscopic to the molecular sizes). Another major challenge in numerical simulations is the so-called contact line singularity that arises when a continuum description of moving contact lines is used in combination with a no-slip boundary condition at the liquid-solid interface. Because of this singularity, the continuum description is untenable below a certain scale. Thus a transition to a different, nonsingular physics must occur as the scale is reduced. The most obvious such transition is the appearance of molecular effects at nanometer scales. However a variety of other “microscopic” contact line physics, some of which would “kick-in” at scales much larger than the nanometer, have been considered in the literature. It is difficult to be exhaustive but these involve precursor film models [7], diffuse-interface models [8, 9, 10, 11] and the related issue of evaporation [12], interface formation models [13, 14], and surface roughness [15]. The reader may find references to other mechanisms in review papers [16, 17, 18, 19]. Slip of the contact line is of particular interest as a possible physical mechanism to allow motion of the contact line on the microscopic scale, mostly because it conveniently does not require to change the Navier–Stokes equations, see e.g. [20, 21, 22, 23, 24, 25, 26, 27, 28]. However, numerical simulations involving slip-length modeling are unfeasible in most physical problems since the true slip should be related to molecular interactions between the liquid and the solid substrate [16], which based on experimental measurements lies in the nanometer range [29, 30]. Thus a regularization of the numerics based on the slip length leads to computationally inconvenient large ratios of scales. On the other hand, if no regularization is performed or if the slip length is dependent on grid spacing, the moving contact line solutions become themselves dependent on grid spacing (see e.g. [27, 31, 32]).

This paper pursues several goals. One is to describe the method used in the numerical framework, Gerris [33] to impose the contact angle and the related “banded” advection of the Volume-Of-Fluid method. (Although the contact angle in a version of Gerris has been described in [34, 35], the exact method used in the mainline distribution of the Gerris code has not been described before.) Another goal is to attempt to extend numerically the theory of Eggers and his coworkers [36, 37, 38] for the dewetting transition to the case of finite microscopic angles. One consequence of this analysis applied to the numerical case is a precise description of the behavior of any numerical method in the vicinity of the contact line. Another consequence is to aid models such as those in [27, 39] or [40] that are used to perform grid-independent models.

In this paper we do not attempt to implement a sophisticated dynamic contact angle model, but instead take an existing, simple numerical method already documented in [34, 35] for static cases and apply it “as is” for dynamic simulations. This is a kind of “implicit modeling” approach similar to the implicit subgrid scale modeling frequent in Large Eddy Simulations of turbulence. This approach has the merit of simplicity, and it then remains to assess how this “numerical boundary condition” affects the flow.

We study a specific, complex physical problem: the dewetting transition. In a number of applications, the interface is forced to move along a solid in a manner that can result either in a receding contact line or in the formation of a thin film on the solid. One example of such a flow is the withdrawing-tape experiment whose geometry is illustrated on Fig. 1. A solid substrate is withdrawn on the left from a viscous liquid pool of quiescent liquid. The interface may either sustain a stationary state meniscus, if below a critical capillary number, ${\rm Ca}_{cr}$, or continue to move up the substrate until depositing a thin film to arbitrary heights. The latter is called a Landau–Levich–Derjaguin (LLD) film [41, 42]. This transition can be understood in terms of the imbalance between the surface tension, gravity and viscous forces that leads to the vanishing of the contact line. The analysis of the transition process on a partially wetting substrate is however complicated due to the singularity of the moving contact line. On one hand, Eggers and his coworkers, in a series of papers [36, 37, 38], provided a hydrodynamic prediction, based on the lubrication approximation theory, of the critical capillary number $\mbox{Ca}_{cr}$. Cox [23] and Voinov [43] described, on the other hand, how the singularity drove a peculiar curved form of the fluid wedge at small ${\rm Ca}$. We use these theories to predict the numerically observed transition. We shall show that 1) the microscopic length $r_{m}$ is entirely described by a gauge function depending only on the equilibrium contact angle imposed numerically and 2) that the knowledge of this numerical gauge function can be used to mimic the effect of actual subgrid scale microscopic physics. This second point is related to the notion of grid-independent simulations in [27, 39].

The rest of this paper is organized as follows. In Sec. 2, we describe the physical problem and explain our choice of parameters. In Sec. 3, we describe the computational setup and the numerical method for the implementation of the contact angle. In Sec. 4, we report the results of our simulations including the method for determining the critical capillary number ${\rm Ca}_{cr}$. In Sec. 5, we introduce the main new theoretical concept of this paper, the gauge function $\phi$. We discuss the general theoretical features of the transition and study it in two cases, the small angle, small viscosity ratio case, free-surface case, where Eggers’s results can be applied directly and the more general case where the use of matched asymptotics and Airy functions is replaced by heuristic arguments. In Sec. 6, we discuss the connection with theories where the microscopic angle varies with the capillary number as in the Molecular Kinetic Theory [44]. In Sec. 7, we apply the concepts developed above to describe an improved procedure for obtaining grid-independent computations of problems with dynamic contact lines. Finally, we draw our conclusions in Sec. 8.

2 Problem setup: forced dewetting

We consider a solid plate being withdrawn from a liquid reservoir with a constant velocity $V_{s}>0$. The computational domain is $0\leq x,y\leq L$, with fluid 1 occupying $y<h_{0}$ and fluid 2 occupying $y>h_{0}$ at $t=0$ (see Fig. 1). The viscosity and density of fluid $i=1,2$ are $\mu_{i}$ and $\rho_{i}$, respectively. The capillary number is defined as

[TABLE]

where $\mu_{i}$ is the viscosity of fluid $i$ and $\sigma$ the surface tension. We set $L\approx 9\,l_{c}$, where $l_{c}$ is the capillary length $l_{c}=\sqrt{\sigma/[(\rho_{1}-\rho_{2})g]}$ with $g$ the gravitational acceleration. The Reynolds number is then defined based on the capillary length as $\mbox{Re}=\rho_{1}V_{s}l_{c}/\mu_{1}$. Thus ${\rm Re}\,={\rm Ca}\,{\rm Mo}\,^{-1/2}$ where the Morton number is

[TABLE]

We define a “gravity-wave-damping” number as

[TABLE]

Indeed, the wavelength of gravity waves traveling at the same speed as the withdrawing tape is $L_{gw}=\rho_{1}V_{s}^{2}/[(\rho_{1}-\rho_{2})g]$ and it can be connected to $N_{G}$ through $N_{G}=\rho_{1}L_{gw}V_{s}/\mu_{1}$. Thus $N_{G}$ is also the Reynolds number based on the wavelength of gravity waves. It is related to the capillary length Reynolds number by

[TABLE]

We use several setups for the simulations, presented in Tab. 1. In Setups A and B, for relatively more efficient computations, we set the ratios of physical properties to moderate values with the viscosity ratio $\mu_{1}/\mu_{2}=1$ and the density ratio $\rho_{1}/\rho_{2}=5$. The other parameters are chosen in the following way. In Setup A, the number $N_{G}$ is arbitrarily set to $N_{G}=25/64$. (It is the value corresponding to the arbitrary choices of $\rho_{1}=5,\rho_{2}=1,g=16,\mu_{1}=1,V_{s}=1$.) Thus from Eq. (1) the Reynolds number based on $l_{c}$ varies as

[TABLE]

As a result the Reynolds number based on $l_{c}$ increases as the capillary number decreases. Varying the Reynolds number between [math] (Stokes approximation) and ${\rm Re}\,=3$ has no effect on the results, however increasing ${\rm Re}\,$ beyond this value introduces significant inertial effects and interface oscillations that modify the conclusions of our investigations. At small ${\rm Ca}$, we therefore switch to Setup B, where $N_{G}$ is free to vary and the Reynolds number based on $l_{c}$ is fixed to ${\rm Re}\,=1$. In a final set of simulations (Setup C), we keep $\rho_{1}/\rho_{2}=5$ and ${\rm Re}\,=1$ but let $\mu_{1}/\mu_{2}=50$. This allows to bring the simulations somewhat closer to air/water conditions without encountering the numerical problems arising with very large density ratios.

We begin the simulations by considering a flat interface between the two fluids initially at the height $h_{0}\approx 3.1\,l_{c}$. A no-slip boundary condition is prescribed at the substrate ($x=0$). Symmetry boundary conditions are imposed on the right ($x=L$), top ($y=L$), and bottom ($y=0$) boundaries of the domain. We note that we have checked that the results are insensitive to the computational domain size.

In the neighborhood of the contact line, no specific choice of parameters is required except the equilibrium or static contact angle that is specified in the numerical model. It is also expected that the numerical model leads to a solution varying continuously with the withdrawing tape velocity, so that the contact angle tends to the static contact angle as ${\rm Ca}$ tends to zero. How simulations with a contact angle are performed is described in the next section.

3 Numerical model

We use Gerris [33, 45, 46] to numerically solve the Navier–Stokes, continuity and density equations,

[TABLE]

[TABLE]

[TABLE]

respectively. Here, u is the velocity field, $p$ the pressure, $\rho=\rho_{1}\chi+\rho_{2}(1-\chi)$, $\mu=\mu_{1}\chi+\mu_{2}(1-\chi)$, $\kappa$ is the interface curvature, n the normal to the interface (pointing from fluid 1 to fluid 2), $\delta_{s}$ the delta function centered at the interface, $\rho\textbf{g}=-\rho g{\hat{y}}$ the body force due to gravity and ${\hat{y}}$ the unit vector in the $y$-direction, and $\chi$ ($=1$ in fluid 1 and [math] in fluid 2) the characteristic function, where $\delta_{s}\textbf{n}=\nabla\chi$. Note that Eq. (5) is equivalent to

[TABLE]

which is solved using the Volume-Of-Fluid (VOF) interface capturing method [46, 47, 34, 48]. The Continuous Surface Force (CSF) method is used for the surface tension force with curvature computed using the Height-Function method [46, 48]. Viscous forces are implemented using a partially implicit method described in [49].

It is useful to describe the procedure used near the contact line in the main distribution of the Gerris code, which we used in the computations reported below. First, without giving the full details that can be found in the references, let us outline the procedure used away from the contact line. There, the VOF discretization of Eq. (6) consists in the definition of a variable $C_{i,j}$ on each grid point $i,j$ that is equal to the volume fraction of the reference fluid, fluid “1”, in the cell. The VOF method proceeds in two steps, first the reconstruction of the interface followed by its advection. In the first part of the reconstruction step, the interface normal ${\bf n}=(n_{x},n_{y})$ in cell $i,j$ is determined from the values $C_{i,j}$ in neighboring cells, using the “Height-Function” method described below (see e.g. [34, 46]) or the “mixed Youngs-centered” (MYC) method (see e.g. [50]), if the former method fails. In the second part of the reconstruction step, the position of a linear segment representing the interface in the cell is determined using elementary geometry (see [50]) from the knowledge of ${\bf n}$ and $C_{i,j}$. Thus the equation of the segment is written

[TABLE]

where the scalar $\alpha$ characterizes the position of the interface. The knowledge of ${\bf n}$ and $\alpha$ is then used in turn in the second, advection step, where the interface is displaced by the fluid velocity field. On Fig. 2, the standard “Lagrangian-Explicit” (see [50]) advection step is represented. It is useful to describe the advection process in some detail. The collocated velocities $u_{i,j},v_{i,j}$, defined on cell centers, are used to compute an auxiliary set of velocities $u_{i+1/2,j}$ and $v_{i,j+1/2}$ on different cell face centers using a projection method. The determination of the motion of the piece of interface shown on Fig. 2 is identical to the “Lagrangian Explicit” or “CIAM” method [46, 50]. However a recent implementation in Gerris uses a “banded” advection. In this approach, the cell is subdivided in $m$ equal bands, the default being $m=4$ as on Fig. 3, and the advection is performed separately in each band. After that, the bands in the cells are aggregated to produce the final volume fraction. This ensures a better representation of shearing or rotating velocity fields, while volume conservation is enforced by the requirement that the average of the horizontal velocity in the bands, for example, on the right side of the cell in Fig. 3 (left panel) to be equal to the face center velocity $u_{i+1/2,j}$.

We next present a test case to clearly illustrate the improvement when using the “banded” advection [51]. In this test case, a straight interface is advected by a pure shear flow, at a shear rate of $1$, in a $1\times 1$ computational domain. The exact solution is simply a rotation of the interface around the center of the domain. We note that both the interface and the velocity field are described exactly by a second-order method. Fig. 4 illustrates the rotation of the straight interface under the shear flow. The green segments are the VOF reconstructed interfaces obtained with $m=1$ and the red segments are when using a “banded” advection method with $m=4$ bands. Fig. 5 illustrates the evolution of the error with time. For $t=1$, the interface is at $45$ degrees and the errors in fluxes cancel out exactly.

In order to compute capillary forces, we use the Height-Function method, in which the local height of the interface is computed from summing over a column of cells with

[TABLE]

where $j_{1}$ is the index of the bottom cell in the column and $j_{2}$ the index of the top cell. When the bottom cell is full ($C_{i,j_{1}}=1$) and the top cell is empty ($C_{i,j_{2}}=0$), and there is a single interface in the column, the height $h_{i}$ approaches the exact interface height to second order [52]. Using finite differences of the local Height-Function then provides the curvature, as well as the interface normals, used to compute the surface tension force by the CSF method (see [46, 48]).

The numerical method given in this paper has often been tested only for regular grids. For example tests of the curvature estimation by the height function method were performed only on circles or spheres where the curvature is uniform. A uniform curvature is naturally associated in these tests to a uniform grid. However, there are some indirect tests, where near a singularity, a range of scales of curvature appear, such as the pinching thread test in [46]. The contact line flow is another such case where the singular flow near the dynamic contact line forces a wide range of curvature. Tests of more complex flows as performed in this paper are also useful to assess the accuracy of capillary force modeling on non-uniform quadtree grids.

Near the contact line, we consider a cell $i,j$ containing the contact line C as shown on Fig. 6(a). The solid-fluid boundary is located exactly on the lower boundary of this cell. (Other locations for the solid boundary have not been explored by the authors.) For such a computational cell, the normal vector is recovered trivially from the specified contact angle $\theta_{\Delta}$ as $n_{x}=-\sin(\theta_{\Delta})$, $n_{y}=\cos(\theta_{\Delta})$. The value of $\alpha$ in Eq. (7) is then obtained using elementary geometry. The interface can then be linearly extended into the solid cell $i,j-1$ as shown by the dashed line on Fig. 6(b). It is assumed there is no other contact line in the immediate neighborhood. The computation of the normal in cell $i+1,j$ is not immediately possible but this difficulty is easily circumvented by assigning “ghost” $C_{i,j}$ values to the cells in the first solid layer $j-1$.

Then as shown on Fig. 6(c), the heights $h_{i}$ and $h_{i+1}$ can be reconstructed in columns $i$ and $i+1$. To construct the height $h_{i-1}$ in the column beyond the contact line a first-order extrapolation of the form $h_{i-1}=h_{i}-\tan(\theta_{\Delta})\Delta$ is used, where $\Delta$ is the local grid size.

The heights over the three columns $i-1,i,i+1$ can be used to compute the curvature in the cell containing the contact line C. This is a good approximation when both the interface slope and its curvature are small enough. Alternatively, it is possible to fit a parabolic approximation of the interface through the two heights $h_{i}$ and $h_{i+1}$ and the contact point C of Fig. 6(a) computed from Eq. (7).

For vertical interfaces one uses “widths” $w_{j}$ instead, as shown on Fig. 6(d). Finally when neither heights nor widths are available, several other strategies are used to compute the curvature, as outlined by Popinet [46], using either a mixture of heights and widths (the so-called mixed-heights method) or if the mixed-heights method fails, a polynomial fit to the mid points of the segments in each cell. These alternative strategies can be adapted to the vicinity of the contact line, provided extrapolations of the volume fraction $C_{i,j}$ and slope in cells $i-1,j-1$ and $i,j-1$ are used. It is important to remark that the Height-Function method for curvature (without mixed heights) always provides a curvature for $\kappa\Delta$ small enough. The critical value of $\kappa\Delta$ was estimated independently by one of the authors through tests on a large number of random circles showing that the minimum value of $\kappa\Delta$ at which the Height-Function always works is about $0.06$.

Once the interface positions and the curvature are computed, there is no special difficulty in computing the velocity field using the standard methods. No special provision is made for the discontinuity of velocities or the divergence of viscous stresses and pressures, which are computed as elsewhere in the domain using finite volumes and finite differences.

We note however that the Gerris code uses the staggered (face) velocities to advect interface pieces, and therefore, in the contact line cell C on Fig. 6, the tangential (face) velocities used for advecting the interface in that cell will not be equal to the solid velocity, since they are defined half-a-cell above the solid boundary which is at $j-1/2$ on Fig. 6. Intuitively, this allows a kind of effective slip. However the amount of slip is reduced by the banding method described above and shown in Fig. 3.

It is interesting to note that at the scale $\Delta$ of the cell $i,j$ the interface and the velocity field are represented in a coarsely averaged manner that is very far from capturing the flow reversal expected inside the fluid 1 wedge (although on scales larger than $\Delta$ this flow is well captured). Whether it is possible or desirable to have a more sophisticated discretization approach near the contact line is beyond the scope of this paper.

4 Results: transition to film formation in forced dewetting

We focus on the problem of a partially wetting substrate withdrawn from a liquid reservoir. We find two parameter ranges: first, the stationary regime, for which the contact line motion along the substrate can evolve to a steady state, as $\tau\to\infty$, where the nondimensional time is $\tau=V_{s}t/l_{c}$; second, the unsteady regime for which a steady state solution cannot be found and the contact line height continues to increase, covering the substrate by a film. All the results presented in this section are for Setup A, unless stated otherwise. We begin by presenting various scenarios characterized by different nondimensional grid sizes, $\Delta/l_{c}$, and the imposed contact angle, $\theta_{\rm{}_{\Delta}}$. Fig. 7 shows the instantaneous contact line height, $h_{CL}$, from the reference height, $h_{0}$, nondimensionalized by $l_{c}$. Fig. 7 shows that depending on $\Delta/l_{c}$, different equilibrium configurations can be obtained. Also, it shows that when the contact angle is decreased, the contact line is raised to a new equilibrium height for the large grid sizes while a stationary meniscus cannot be achieved for the smallest grid size. Here we take $0.007<\Delta/l_{c}<0.057$. We note that for the capillary numbers that we consider (from $0.001$ to $0.1$), it is very difficult to have a larger range of grid sizes since one needs to satisfy $\Delta\ll\ell\ll l_{c}\ll L$, where $\ell$ is the thickness of the film at the transition, of the order $\ell={\rm Ca}^{2/3}l_{c}$. We however note that the Gerris code uses quadtree grids (or octree in 3D), that allow to refine the grid where necessary. This is a very useful feature for dynamical contact line problems as it allows to use a very small grid size $\Delta$ in the immediate neighborhood of the contact line and a larger grid size elsewhere.

Next, we elaborate on the stationary state results in detail. Figs. 8(a)-(d) show the nondimensional stationary contact line height $(h_{\infty}-h_{0})/l_{c}$, where $h_{\infty}=h_{CL}(\tau\to\infty)$, as a function of the nondimensional mesh size, $\Delta/l_{c}$, for various capillary numbers and for various contact angles, $50^{\circ}\leq\theta_{\rm{}_{\Delta}}\leq 90^{\circ}$. The results show the grid sizes, $\Delta/l_{c}$, where a stationary meniscus forms and no film is deposited on the substrate, for the range of considered $\theta_{\rm{}_{\Delta}}$.

More importantly, the results show that the computed height of the contact line is a function of $\Delta/l_{c}$ and that for small enough $\theta_{\rm\Delta}$ no steady state menisci can be attained when $\Delta/l_{c}\leq 0.014$. This lower limit of $\theta_{\rm\Delta}$, for which steady state contact lines can be achieved, gets larger as ${\rm Ca}$ is increased. The results clearly depend on the chosen value for the smaller mesh size near the contact line, an effect that is expected and will be explained below. This dependence of the results with mesh refinement becomes more marked as ${\rm Ca}$ is increased or $\theta_{\rm\Delta}$ is decreased.

Next we present the results, for which the contact line continues to move upward and a liquid film is then deposited on the substrate. We can understand the onset of film deposition, i.e. the forced dewetting transition, as when the balance between the surface tension and viscous forces close to the contact line region can no longer hold, resulting in wetting failure. (We note that gravity is only involved, asymptotically, in the outermost region.) At the transition, however, we will need to allow the computations to run for a very long time in order to determine when the contact line motion along the wall cannot reach a stationary state. For efficient and accurate determination of the numerical values of the transition ${\rm Ca}_{cr}$, we propose a procedure as follows. Fig. 9(a) shows contact line heights as a function of time for various values of ${\rm Ca}$, when $\theta_{\rm\Delta}=90^{\circ}$ and $\Delta/l_{c}=0.007$. As shown, for sufficiently small ${\rm Ca}$, a stationary meniscus can be reached while for large ${\rm Ca}$, the contact line height keeps increasing. We then use the information in Fig. 9(a) to obtain Fig. 9(b). We then pick the transition capillary number, for which the relative velocity of the contact line, $|d(h_{CL}/l_{c})/d\tau|$, does not go to zero as a function of $\tau$. This critical capillary number ${\rm Ca}_{cr}$ is depicted in red in Figs. 9(a) and (b). Using the procedure above, we can therefore determine ${\rm Ca}_{cr}$ with a very good precision.

Figs. 10(a)-(b) show the contact line height, $(h_{CL}-h_{0})/l_{c}$, as a function of nondimensional time, $\tau$, for two mesh sizes, $\Delta/l_{c}$, when varying the wall contact angle $\theta_{\rm\Delta}$, for a fixed ${\rm Ca}$. As shown, the transition from a stationary meniscus not only depends on $\theta_{\rm\Delta}$, but also on $\Delta/l_{c}$. As illustrated, the transition occurs at a larger $\theta_{\rm\Delta}$ for smaller $\Delta/l_{c}$. This observation begs a further exploration of how the critical capillary number depends on the contact angle and the grid size. We study these effects in what follows.

To shed more light on the transition mechanism, we analyze the flow, for both when a stationary meniscus forms and when a steady state contact line cannot be attained. Fig. 11(a) provides an example of a stationary meniscus for ${\rm Ca}=0.043$, $\theta_{\rm\Delta}=90^{\circ}$, and $\Delta/l_{c}=0.014$ (for this set of parameters, $\mbox{Ca}_{cr}=0.52$). The inset shows the magnified flow field and the pressure distribution. Fig. 11(b) shows a magnified view of the contact line region and the computational mesh. The fine structure of the flow field and the pressure distribution in the contact line region are illustrated. As shown, large gradients of velocity and pressure necessitate a high mesh resolution around the contact line region. As illustrated, the interface is highly curved close to the contact line, leading to an intensified pressure gradient around that region, while the pressure gradient remains weak outside the vicinity of the contact line, leading to gentle bending of the interface away from the contact line.

Fig. 12 shows an example of when the contact line cannot attain a steady state, leading to the formation of a film deposited on the substrate, for ${\rm Ca}=0.048$, $\theta_{\rm\Delta}=60^{\circ}$, and $\Delta/l_{c}=0.014$ (for this set of parameters, $\mbox{Ca}_{cr}=0.024$). The figure shows a typical evolution of the interface and the transition to the film. Figs. 12(d)-(f) also show the sagging of the interface behind the contact line after the film formation. The insets of Figs. 12(a)-(c) show the magnified flow field and the pressure distribution. Figs. 12(a)-(c) also show a further magnified view of the contact line region and the flow streamlines. The fine structure of the flow field and the pressure distribution in the contact line region are illustrated. The results reveal the strong pressure gradients close to the contact line region; when capillary forces can no longer balance this strong pressure gradient, the consequence is the wetting failure.

We now turn our attention to computing the critical capillary number, ${\rm Ca}_{cr}$, of the dewetting transition, as a function of $\Delta$, $\theta_{\Delta}$, and $q$. We carry out a series of simulations when incrementing the capillary number, holding $\Delta$ and $\theta_{\Delta}$ fixed, until we find ${\rm Ca}_{cr}$. We repeat the procedure for Setups A, B, and C in Tab. 1. We note that the smallest mesh size we choose here is $\Delta/l_{c}\approx 0.004$; such an unprecedented fine grid is only made possible owing to the adaptive mesh refinement capability of the Gerris code, that allows focusing the grid refinement at around the interface and in particular the contact line region while the grid size is free to adapt elsewhere. Fig. 13 presents ${\rm Ca}_{cr}$ as a function of $\Delta/l_{c}$ for a range of $\theta_{\rm\Delta}$ for Setups A, B, and C. The symbols are direct simulation results and the solid lines are drawn to guide the eye through a set of points corresponding to the same $\theta_{\rm\Delta}$. The error bars represent the increment in the capillary number in the procedure of finding ${\rm Ca}_{cr}$. These results illustrate the logarithmic dependence of the dewetting transition on the mesh size for all the ranges of the contact angles that are considered here for various Setups. Moreover, as expected, for large capillary numbers, Setup A and B behave more or less the same. In addition, we show the effect of the viscosity ratio on computed ${\rm Ca}_{cr}$; in general, lowering the viscosity of fluid 2 results in an increase in ${\rm Ca}_{cr}$ compared to when the viscosity of both fluids is matched, with this effect being more prominent at large contact angles. As shown, a decrease in contact angle promotes the dewetting transition. In the next section, we will develop an improved understanding of the onset of this dewetting transition. We study the effects of the contact angle and the grid size on $\mbox{Ca}_{cr}$ and give the scaling of it with $\Delta/l_{c}$.

5 Hydrodynamic theories of the dynamic contact line and the dewetting transition

5.1 Generalities about the asymptotic description

To interpret our numerical results, we will use a theoretical framework that extends the work in [36, 38]. The theory in [36, 38] is valid for small capillary number and equilibrium contact angles, large viscosity and density ratios, and several specific slip length models. Since our numerical approach is not based on a slip length model, we cannot apply this theory directly. However, as we shall see, it is straightforward to replace the slip length analysis, by an asymptotic matching, with the small length scale behavior, using the theory of Cox [23]. In the general physical case, it is also assumed that any microscopic physics with a sufficiently small length scale can be represented in Cox’s theory.

In the small ${\rm Ca}$ limit, the problem is analyzed by asymptotic matching of three or four different regions:

1. A region I near the contact line, where microscopic effects dominate in the physical reality and numerical effects dominate in our simulations. The size of that region is the microscopic scale $\lambda$. We make no assumptions about region I except that it is limited in length scales to a microscopic scale $\lambda$ that remains much smaller than all other relevant length scales.

2. A region II at length scales large enough that continuum mechanics and the no-slip condition on the wall hold, but close enough to the contact line that a special logarithmically-scaling solution holds, as discussed by Cox in [23]. This region starts at lengths $\ell\gg\lambda$ and ends at lengths $\ell\ll{\rm Ca}^{1/3}l_{c}$. In that region, capillary forces are balanced by viscous drag forces. The interface thus bends, and as shown below, it has a curvature proportional to ${\rm Ca}/r$ where $r$ is the distance to the contact line.

3. An outer region III, where viscous effects are negligible and surface tension balances gravity. This region scales like the capillary length $l_{c}$. The solution in region III is the famous static meniscus solution [53]. We define the “apparent contact angle” $\theta_{\rm a}$ as the angle seen in meniscus variables that is on scales such that $y\sim l_{c}$. The meniscus solution quadrature gives the curvature

[TABLE]

The general properties of the solution can be seen as follows at small capillary numbers. At the overlap of region I and region II the angle is the microscopic angle noted $\theta_{e}$. It is often assumed that it is the equilibrium contact angle, but this is true only if the microscopic physics in region I are close to equilibrium. It is quite possible to assume that $\theta_{e}$ differs from the equilibrium angle in a manner determined by experiments, or by microscopic theories of the physics beyond the scope of this paper. For small capillary numbers, the bending in region II is sufficiently small that the angle changes little from region I to region III and thus $\theta_{\rm a}=\theta_{e}$. This condition closes the problem and gives in particular the height of the contact line $h_{cl}=l_{c}\sqrt{2-2\sin\theta_{e}}$. As the capillary number increases, the curvature and the bending increase in region II. The apparent contact angle decreases below $\theta_{e}$ while the height of the contact line increases.

It can thus be hypothesized that the critical condition for dewetting can be expressed as the condition of vanishing apparent contact angle. For a vertical plate, this means that the inner solution in region II matches into the outer solution in the meniscus region III with a vertical slope or yet that at the inflection point, the slope vanishes. Thus a region IV, overlapping regions II and III is found in which the slope is small. This hypothesis is implicit in the works of Derjaguin [42] and Landau and Levich [41] since their assumption of a lubrication theory implies a near-vertical slope in the matching region, the existence of region IV and a zero apparent contact angle. Moreover Landau and Levich [41] have shown that in region IV and at small capillary number, the film thickness behaves as $h\sim{\rm Ca}^{2/3}l_{c}$, the famed Landau-Levich film thickness.

This hypothesis of transition at zero apparent contact angle has been considered by Eggers in the small angle case [36, 38], but has not yet been investigated for large angles or when the free-surface model has to be replaced by an interfacial model at non-vanishing viscosity ratio $q$.

The matched asymptotics allowing to connect the regions defined above are valid upon the hypothesis that a parameter $\epsilon$ is small enough. This parameter is not simply ${\rm Ca}$ as the Landau-Levich film scaling would seem to imply, but $\epsilon\sim{\rm Ca}/G(\theta_{e})$ where $G(\theta_{e})$ is the function introduced in [23] and defined below. It can also be shown that at the transition, an equivalent estimate of the small parameter is $\epsilon\sim 1/\ln(l_{c}/r_{m})$. This value of $\epsilon$ is not exceedingly small, which has consequences which will be spelled out below.

5.2 The theory of Cox

We thus focus on the analysis of region II. In this region, the wedge solution of Huh & Scriven is assumed [20]. It is then possible to use it as done by Cox [23] to obtain the variation of pressure in the wedge and thus by Laplace’s law the variation of curvature. After integration, one may obtain the variation of slope. Cox’s solution in his “intermediate region” identified with our region II, from his Eqs. (7.13) and (7.18) is

[TABLE]

where $\lambda_{c}$ is a characteristic scale for microscopic effects, $\theta_{e}$ is the equilibrium angle in agreement with our assumption above for the ${\rm Ca}\to 0$ limit, $f$ and $G$ are defined below, and $q=\mu_{2}/\mu_{1}$ is the viscosity ratio. $Q_{i}$ is an integration constant that is obtained by matching with region I and thus depends on region I characteristics. Since we actually write the first order of an expansion in small ${\rm Ca}$, higher-order terms exist whose form is however unknown to the authors.

Finally the functions $f$ and $G$ are given by

[TABLE]

and

[TABLE]

In some specific asymptotic limits, special forms of $G(\theta)$ can be used. For $\theta\ll 1$ and $q=0$ as in thin liquid film in air situations, it can be shown that

[TABLE]

Also, Sheng and Zhou [54] show that $G(\theta)-G(\theta_{e})\approx(\cos{\theta_{e}}-\cos{\theta})/5.63$ when $q=1$ and $|\cos{\theta}|<0.6$. In this work, we account for $G(\theta)$ using Eq. (10) directly.

A more convenient notation is to introduce a “microscopic length scale”

[TABLE]

as in Voinov [43, 55] so that in the general case

[TABLE]

The microscopic length scale $r_{m}$ encompasses the properties of the microscopic region. Thus the effect of the corresponding small scales can be summarized with two parameters $\theta_{e}$ and $r_{m}$.

Eq. (9) is very attractive because it provides a universal description giving the dependence of $\theta(r)$ on ${\rm Ca}$, to the leading order, without any specification of the microscopic physics or any necessity to calculate the details of the macroscopic flow in the outer region. In this paper, we also claim that it is valid when the concept of “microscopic physics” is replaced by “numerical scheme”; the numerical scheme introducing deviations from the continuum equations at scales of order $\lambda_{c}\sim\Delta$.

The microscopic length $r_{m}$ can be obtained in a number of special cases, or by dimensional analysis. Eggers [36] has assumed several slip length models. For the Navier-slip model, and using our notations above,

[TABLE]

where $\lambda$ is the slip length. Using lubrication theory valid for small slopes or angles, one obtains the equation [56]

[TABLE]

where we introduced a variable $\zeta$ that measures distance away from the contact line along the direction of the solid surface. We also note that $\eta(\zeta)$ is the local thickness of the film and we use the notation $\eta^{\prime}=d\eta/d\zeta$, $\eta^{\prime\prime}=d^{2}\eta/d\zeta^{2}$, etc.. A solution of Eq. (15) can be obtained [57] in the form

[TABLE]

where

[TABLE]

(Notice that equation (33) of [57] misses the above factor of $\mathrm{e}$ because of an incorrect derivation from equation (32) of [57]). This can be compared to the solution obtained by Cox in [23]. From Eq. (11) and Eq. (13), using the fact that for small angles $r\sim\zeta=h_{CL}-y$ and $\theta\sim\eta^{\prime}(\zeta)$, one obtains

[TABLE]

For ${\rm Ca}\ll\theta^{3}_{e}$, the above yields again Eq. (16). Thus the length scale $r_{m}$ can be written as

[TABLE]

where $\phi$ is a gauge function depending on $\theta_{e}$ equal, at first order, to

[TABLE]

Different slip-length models give different gauge functions $\phi$. Moreover, alternate models of the microscopic physics beyond slip length, such as Van der Waals forces or phase field, would still result in an expression of the form Eq. (13), with $\phi$ having a dependence on $\theta_{e}$ and other microscopic model parameters to be determined.

In the numerical case, the equivalent expression for the microscopic length scale in region II is

[TABLE]

where the gauge function $\phi_{\rm num}$ may depend on $q$, $\theta_{e}$ and the details of the numerical implementation such as the number of shear bands $n$, any specific numerical slip etc. In what follows, for simplicity, we shall drop the “num” subscript and leave the dependence of $\phi$ on the viscosity ratio $q$ implicit. Comparing expressions (19) and (21) leads to the relation

[TABLE]

which shows that the numerical model can be interpreted as having an effective slip. Using the definition of the numerical gauge function Eq. (21), Cox’s solution (Eq. (9)) can be rewritten in the form

[TABLE]

This form is verified numerically in the next Section.

5.3 Numerical verification of Cox’s theory and matched asymptotics

We next present the numerical results of computed curvature and interface shapes, for all the three Setups A, B, and C, and compare the results with the theoretical predictions presented in previous section. We consider $15^{\circ}\leq\theta_{\rm\Delta}\leq 110^{\circ}$ and $0.004\leq\Delta/l_{c}\leq 0.016$. In what follows, we show that in general, the curvature is negative and higher close to the contact line and it decreases away from the contact to zero, where it then turns and becomes positive to assume the curvature of the meniscus and eventually approaches the zero curvature of the flat film far away from the wall. This means that the curvature effects are more prominent close to the contact line, and that there always exists an inflection point. We show that interface bending becomes significant as we increase the contact angle and/or decrease the mesh size. The divergence of the curvature near the contact line is as expected and discussed in the previous section. We also present the interface apparent slope along with the curvature to enlighten the connection of the interface apparent zero slope with the presence of an inflection point. Indeed, we notice that the angle is near zero at inflection point when the curvature vanishes, as assumed in Sec. 5.1. The interface slope results also confirm that we approach the imposed numerical contact angle, $\theta_{\rm\Delta}$, as we approach the contact line. We nondimensionalize the curvature by $l_{c}$ and present the results as a function of the vertical height from the contact line, i.e. $\zeta/l_{c}=(h_{CL}-y)/l_{c}$, (see Fig. 1 for the illustration of the distance $\zeta$).

Figs. 14 to 22 depict the curvature plots (black symbols) along with the angle (red symbols) $\theta_{a}=\theta(\zeta/l_{c})$ that the interface makes with the substrate as a function of $\zeta/l_{c}$. We find that the maximum curvature is obtained at the contact line and that this maximum curvature increases with decreasing the mesh size. As shown, there is a substantial bending of the interface for $\theta_{\rm\Delta}=110^{\circ}$; for example, for $\Delta=0.004$, curvature is tending to zero around $-0.75l_{c}$. However, the maximum curvature does not reach the asymptotic value of $\sqrt{2}l_{c}^{-1}$ when the asymptotic parameter $\epsilon={\rm Ca}/G(\theta_{e})$ is not small enough. As shown, for small angles, where the smallness of $\epsilon$ is guaranteed, the maximum curvature approaches $\sqrt{2}l_{c}^{-1}$. Figs. 14 to Figs. 22 also show that the angle is indeed near zero at inflection point when the curvature vanishes as assumed in Sec. 5.1.

We also compare the computed interface slope close to the contact line region with the prediction of the Cox’s theory by plotting relation (23) on Figs. 14 to 22. We plot the computed (black symbols) $G(\theta_{e}=\theta_{\Delta})-G[\theta(r)]$ versus ${\rm Ca}\ln(r/\Delta)$, using the computed interface slope, $\theta(r)$, at the distance $r=\sqrt{x_{i}^{2}+y_{i}^{2}}$, where $x_{i}$ and $y_{i}$ are the $x$ and $y$ coordinates of the cell center of an interfacial cell relative to the contact line coordinates. Figs. 14 to 22 also show the comparison of the computed (black symbols) $G(\theta_{e}=\theta_{\Delta})-G[\theta(r)]$ versus ${\rm Ca}\ln(r/\Delta)$ with the prediction of Eq. (23) using the best fitted value of $\phi$ (solid lines). However, we also find three to four aberrant grid points at the smallest values of $r$ corresponding probably to the limited accuracy of our numerical method at the grid scale near the contact line. Moreover, if the relation (23) from the Cox-Voinov theory were exact, the slope of $G(\theta_{e}=\theta_{\Delta})-G[\theta(r)]$ versus ${\rm Ca}\ln(r/\Delta)$ should be unity. In fact, there is a deviation from unity which is largest at small angles ($15^{\circ}$ and $30^{\circ}$) and smallest at large angles. This is particularly significant since at small angles and small $q$ the theory has a simple derivation from lubrication theory and the assumption of a parabolic flow in the thin liquid wedge. A possible explanation is that this theoretical flow is not well approximated by the numerical method. Indeed, for a small angle, there is a large region where there are very few grid points across the thin liquid wedge, with the boundary conditions on the solid and the free surface being imperfectly approximated by the finite volume method and the interpolations used for viscous stresses. This type of fit thus appears to be a good test of the accuracy of the method in the vicinity of the contact line.

The fitting also yields a best value for $\phi$ in the range where Cox’s theory is expected to be asymptotically valid $\Delta\ll r\ll r_{\rm max}$. The expected upper limit $r_{\rm max}={\rm Ca}^{1/3}l_{c}$ will be derived below. We consider the results for $0.004\leq\Delta/l_{c}\leq 0.016$ and show that the fitted values of $\phi$ are similar and thus approximately independent of $\Delta$ when ${\rm Ca}$ is maintained at criticality. We extend the discussion on these results next.

Figs. 23 illustrate the interface profiles corresponding to examples in Figs. 14 to 22. The results show the interface shapes, close to the contact line where $r\ll r_{\rm max}$, for various cases, demonstrating how the interface profiles vary at the transition. The background mesh depicts the adaptive mesh refinement that is utilized for the simulations. The results show the monotonically increasing film thickness that connects to the flat region of the interface, corresponding to the zero apparent interface angle, that is in turn connected to the static liquid reservoir meniscus (only the film near the contact line is shown here). For $\theta_{\rm\Delta}=110^{\circ}$, there is a significant bending of the interface, more markedly close to the contact line. This bending of the interface close to the contact line is better captured at higher mesh resolutions. As shown, for $\theta_{\rm\Delta}\leq 30^{\circ}$, even the smallest grid size appears to be just sufficient, illustrating the challenging computations at smaller angles.

5.4 General case

In the arbitrary angle case or when $q$ does not vanish, one cannot use lubrication theory as done in [36]. However we can proceed by assuming a vanishing apparent contact angle at the transition. For a vertical plate, this means that at the inflection point the slope vanishes ($\eta^{\prime}=0$ when $\eta^{\prime\prime}=0$). We have indeed confirmed this observation by our numerical results. As shown in Figs. 14 to 22 (b,d,f), for the range of contact angles and mesh sizes considered here, $15^{\circ}\leq\theta_{\rm\Delta}\leq 110^{\circ}$ and $0.004\leq\Delta/l_{c}\leq 0.016$, the vanishing slope of the interface coincides with the inflection point at the transition, confirming the criteria for the dewetting transition.

Thus there is a region IV of small slope overlapping with region II and III in which the lubrication approximation is valid. However because the slope is not small everywhere in region II the lubrication approximation does not apply there and asymptotic matching is less obvious than in the literature. One possibility would be to use a generalized differential equation also valid for large slopes as done in [58]. A mapping from the equation valid for large slopes in [58] to the lubrication approximation for small slopes has been performed in [59], allowing to reuse the result in [36]; however, for several reasons explained below, we prefer to perform directly a matching from Cox’s theory valid from region II to region IV.

In what follows, we first perform a heuristic derivation of the central relation of our theory, that gives the critical value of the capillary number ${\rm Ca}_{cr}$ implicitly as function of the other parameters of the problem. We then compare these results to those of lubrication theory obtained for small angles, and to our numerics for small angles.

The heuristic description starts by determining the scale of ${\rm Ca}$ and the thickness of the liquid in region III near the inflexion point. Using the Cox-Voinov analysis, the slope is given by Eq. (13) and we write the condition that it vanishes in region IV. Let the inflection point occur at $\zeta\sim\zeta_{I}$ with $\zeta_{I}\gg r_{m}$. The bending occurs over a large distance $\eta_{I}$ and the “final phase” of the bending where the slope and the curvature are small occupies almost all of region II. Thus region II “looks like” a very thin wedge seen at scale $\eta_{I}$ (paralell to the wall) and $\zeta_{I}$ (perpendicular to the wall). The meaning of all the geometrical quantities $s,r,\eta,\zeta,\theta$ and $\theta_{e}$ is summarized on Fig. 24.

As a result of the thinness, we can approximate $r\sim\zeta$ in that region. Writing the small slope assumption $\theta=0$ at $r=\zeta_{I}$ in Eq. (13) gives

[TABLE]

and ${\rm Ca}_{cr}$ is small if $1/\ln(\zeta_{I}/r_{m})$ is small. The small parameter in the expansion is

[TABLE]

This determines the scale of ${\rm Ca}$ and we now turn to determining the thickness $\eta_{I}$ of the liquid in region III near the inflection point. This is not as trivial as it may seem since in addition to the $\eta_{I}\sim{\rm Ca}^{2/3}l_{c}$ solutions near the inflection point that match the meniscus, there are also solutions with an asymptotically flat interface, $\eta\rightarrow\eta_{\infty}$ as $\zeta\rightarrow\infty$. These solutions have $\eta_{\infty}\sim{\rm Ca}^{1/2}l_{c}$ and involve a balance of gravity and viscosity instead of a balance of surface tension and viscosity. A solution may be found trivially by letting $\eta=$ Constant in the generalized lubrication equation or in the large slope equation in [58]. It is thus a fixed point solution. The following serves to show that one matches to the correct inflection point solution, not to the fixed point solution. Recalling that ${\rm Ca}$ is small, it is possible to deduce that the curvature is still small in most of region II thus $ds/dr\sim 1$ where $s$ is the curvilinear abscissa along the interface. Then

[TABLE]

where we have used the fact that $r_{I}\sim\zeta_{I}$, $d\eta/ds=\sin\theta(s)$ and ${\rm Ca}\simeq{\rm Ca}_{\rm cr}$. After some work it can be shown that this integral behaves as

[TABLE]

while by definition

[TABLE]

We thus have $\eta_{I}\sim{\rm Ca}^{1/3}\zeta_{I}$ as expected for the inflection point in the thin film dewetting transition. In region IV, the length scales are thus $\eta_{I}$ across the film and $\zeta_{I}$ parallel to the film. The curvature thus scales as

[TABLE]

To match region IV and region III, one should have the curvature $\eta^{\prime\prime}$ of same order as the curvature $l_{c}^{-1}$ of the static meniscus. Writing $\eta^{\prime\prime}\sim l_{c}^{-1}$ one obtains

[TABLE]

Identifying $\epsilon$ as above

[TABLE]

where we have introduced a constant $C(q)={\cal O}(1)$.

In the free surface case when $q=0$, we may match region III with region IV in a manner similar to the one in [36], as shown in A, to obtain

[TABLE]

Or equivalently (see 52)

[TABLE]

It is possible that $C(q)$ does not differ very much from $C(0)$ if in region IV the outer fluid stress is negligible. Indeed in the thin layer the stress scales as $\mu_{1}V_{s}/\eta(y)$ while the stress in the outer fluid scales as the much smaller $\mu_{2}V_{s}/y$, realizing an approximation of the free surface condition.

We note that expression (33) is equivalent at small angles $\theta_{e}$ to the one obtained in [36]. It is also equivalent to the expression obtained in [59] at all angles apart for the fact that we make no hypothesis on the behavior of $r_{m}$ at large angles, unlike [59] and [58]. Moreover, as we wrote above, we prefer not to obtain our expression through a generalized lubrication equation, but rather through asymptotic matching. We believe that the effect of some models of the microscopic physics, such as the slip length model, in region I, cannot be found for large angles by asymptotic analysis but would require a full solution of the Stokes equation analogous to that of Huh and Scriven [20], but for different boundary conditions.

We also note that at the transition the condition (30) yields $\eta_{I}\sim{\rm Ca}^{2/3}l_{c}$, and that we used essentially the same matching argument as Derjaguin [42] and Landau and Levich [41]. This can also be compared to the prediction of Voinov [55] which is similar but with a constant $C(0)$ differing from that obtained by asymptotic matching.

Considering now the numerical case, we substitute the numerical expression (21) for the microscopic length scale $r_{m}$ in the critical matching expression (33) to obtain the central relation of our theory

[TABLE]

where we identified $\theta_{\Delta}$ with $\theta_{e}$.

It is interesting to investigate explicit expressions for ${\rm Ca}_{\rm cr}$ obtained by an asymptotic analysis of Eq. (34). To perform this analysis, we let $\delta=1/\ln({l_{c}/\Delta})$. This small parameter is of the same order as $\epsilon$ above but unlike $\epsilon$ is explicitly defined using the parameters of the problem. To give an order of magnitude for $\delta$ and $\epsilon$, if we use, as in the above simulation, a grid size $\Delta/l_{c}=1/2^{8}$, we obtain $\delta\approx 0.18$. Thus our small parameters is small but not exceedingly so. It is useful to introduce the auxiliary parameter ${\hat{\mu}}$ defined as ${\hat{\mu}}=-1/\ln[C(q)\phi(\theta_{\Delta},q)G(\theta_{\Delta},q)]$. The parameter ${\hat{\mu}}$ will be small for transcendentally small $\theta_{\Delta}$. Otherwise if $\hat{\mu}$ is considered order 1, solutions of Eq. (34) expand as

[TABLE]

At small angles, with the expected scaling $\phi(\theta_{\Delta})\sim\theta_{\Delta}$, one gets $\hat{\mu}\sim-1/(4\ln\theta_{\Delta})$. Moreover the higher order terms in Eq. (35) are small if ${\hat{\mu}}\gg\delta$ which is verified even for the smallest $\theta_{\Delta}$ and typical values of $l_{c}/\Delta$ in our calculations. For Setup B, with $\theta_{\Delta}=110^{\circ},\Delta/l_{c}=1/2^{8},q=1$, and using the best fit estimate $\phi\simeq 3$, the first order of expansion of Eq. (35) gives ${\rm Ca}_{\rm cr}\simeq G\delta\simeq 0.0558$, using all the terms in expansion (35) gives ${\rm Ca}_{\rm cr}\simeq 0.0626$, while a full iterative solution of Eq. (34) yields ${\rm Ca}_{\rm cr}\simeq 0.0637$. For Setup C, in the case $\theta_{\Delta}=110^{\circ},\Delta/l_{c}=1/2^{8},q=0.02$, and using the best fit estimate $\phi\simeq 3.5$, the first order gives ${\rm Ca}_{\rm cr}\simeq G\delta\simeq 0.127$, using all the terms in expansion (35) gives ${\rm Ca}_{\rm cr}\simeq 0.169$, while a full iterative solution of Eq. (34) yields ${\rm Ca}_{\rm cr}\simeq 0.132$. Thus in a wide range of parameters, the approximation

[TABLE]

is correct within 15% accuracy.

5.5 Comparison with numerics

In our simulations, we measure ${\rm Ca}_{cr}$ and specify $\Delta$ and $\theta_{\Delta}$. For the two special cases described above, we compare the value of ${\rm Ca}_{\rm cr}$ from full simulations to the solutions of the central relation, Eq. (34). The simulations values are in remarkable agreement with the above solutions. For setup B, with $\theta_{\Delta}=110^{\circ},\Delta/l_{c}=1/2^{8},q=1$, we find numerically $Ca_{cr}=0.0555$ compared to the theoretical value $0.0637$. For setup C, in the case $\theta_{\Delta}=110^{\circ},\Delta/l_{c}=1/2^{8},q=0.02$, we find numerically $Ca_{cr}=0.12125$ compared to the theoretical value $0.132$.

Alternately, to perform the comparison, Eq. (34) can be inverted to provide an estimate of $\tilde{\phi}=C(q)/C(0)\phi$ given the numerical estimate of ${\rm Ca}_{cr}$. The estimate is

[TABLE]

If our theory is valid, the right hand side should not depend on $\Delta$. Moreover if lubrication theory can be used in region IV, as is possibly the case even for $q>0$, the above expression should yield the gauge function $\phi$ defined in relation with Cox’s analysis and estimated numerically above.

In Fig. 25, we plot the values of the RHS of Eq. (37) for Setups A, B and C, along with the computed values from the best fit of $\phi$ presented in Sec. 5.3 (see Figs. 14 to 22). In the near-free-surface Setup C ($q=1/50$), a very clean estimate of the gauge function $\phi$ is obtained at small angles and we find $\phi(\theta_{\Delta})\simeq\theta_{\Delta}$. Numerical results from the direct comparison with Eq. (23) are also in reasonable agreement with the predictions at small angles. Eq. (22) therefore suggests an effective slip length $\lambda\simeq\Delta$. Thus our numerical model may be viewed as having an effective slip of the order of a grid cell.

In Setups A and B, the plot is more scattered indicating deviations from our theory, which is not surprising given the not-so-small value of the asymptotic parameter. However, a rough proportionality to the angle $\theta_{\Delta}$ is observed at small angles. Moreover, for all the Setups, large contact angle data do not follow the linear proportionality to the angle $\theta_{\Delta}$, with the viscosity ratio having a more prominent effect on this deviation at large angles. To shed more light on the behavior of $\phi$ as a function of the mesh size, as well as exemplifying the errors in the computation of $\phi$, in Fig. 26, we plot $\phi(\theta_{\Delta})$ versus $\Delta/l_{c}$. The results show some small variations; this is however expected because the error in estimating $\phi$ also depends on $\Delta$, while decreasing the mesh size does not reduce the transition capillary by much.

6 Dynamic inner contact angle

A possible reinterpretation of Cox’s theory is to consider that the inner region angle is not the equilibrium angle but an angle $\theta_{in}$ depending on the contact line velocity so that

[TABLE]

where dependence on the fluids and solid material properties is implicit. Then Eq. (9) must be rewritten as

[TABLE]

and we can proceed as before. We note however that there are two possibilities.

1. The microscopic angle is discontinuous, that is there are two angles, the advancing angle $\theta_{a}$ and the receding angle $\theta_{r}$, so that $\theta_{r}<\theta_{e}<\theta_{a}$ and that for small ${\rm Ca}$, $\theta_{in}\sim\theta_{a}$ for an advancing contact line and $\theta_{in}\sim\theta_{r}$ for a receding contact line. The conclusion is that the asymptotic analysis is unchanged with $\theta_{e}$ replaced by either $\theta_{r}$ or $\theta_{a}$.

2. The microscopic angle is continuous so that at first order

[TABLE]

with $f^{\prime}_{{\rm Ca}}$ a constant. In this case, and the small ${\rm Ca}$ assumption, Eq. (39) can be rewritten at first order in ${\rm Ca}$ in the form of Eq. (9)

[TABLE]

the change in microscopic angle with ${\rm Ca}$ being absorbed into the new integration constant $Q^{\prime}_{i}$. We note that the relation (40) is in agreement with the Molecular Kinetics Theory (MKT) of Blake and Haynes [44], which states that the contact line is displaced by small random molecular jumps at the surface of the solid due to thermal fluctuations. The average size of these small jumps is determined by the intermolecular interactions between the liquid and solid. The frequency of these random molecular jumps and the average distance of each jump can either be determined empirically or by direct comparison to molecular dynamics simulations. This model postulates that there is an out-of-balance surface tension force, of non-hydrodynamic origin, as a result of the contact line moving on a solid surface. The model then relates the contact line speed to this driving force, resulting in a dynamic contact angle.

The idea of advancing and receding angles is related to the idea of physical roughness of the surface. This may connect to the idea of a discontinuity in the global numerical solution itself. We have conducted tests of the numerical model for advancing menisci for ${\rm Ca}$ down to $10^{-6}$, although the dewetting transition itself was studied only for ${\rm Ca}>10^{-4}$. For very small ${\rm Ca}$, the contact line evolves in an irregular manner, with intermittent spikes in velocity akin to the motion on a rough surface. These spikes have a time periodicity of $T=\Delta/V_{s}$ as if the interface was pinned with spatial periodicity equal to the grid size $\Delta$. However, no such irregularity is observed for ${\rm Ca}=10^{-4}$. Thus we may infer that at least in the range of parameters used in this study $10^{-3}\lesssim{\rm Ca}\lesssim 0.1$, the numerical solution is a continuous function of the contact line position.

7 Consequences for dynamic contact line computations in general

How to conduct realistic simulations of flows including a moving contact line in practice is a difficult problem that has elicited significant research efforts (See [39] for a systematic comparison of some models.) Here we only discuss how our Cox type analysis of numerical models connects with other work on mesh independent simulations. In the computation of a number of problems involving dynamic contact lines, such as droplet spreading, drop impact or drop sliding on inclined plates, the result typically varies with grid resolution because of hydrodynamic effects. One already-suggested approach [39] is to specify numerically the angle

[TABLE]

for some $\theta_{in}$ that may be function of the capillary number as in Eq. (38), and for some microscopic length $r_{m}$. We remind the reader that $\theta_{\Delta}$ is the numerically imposed contact angle while $\theta_{in}$ is an experimentally measured or determined quantity whose meaning is the contact angle at the microscopic scale $r_{m}$. Since $r_{m}$ and $\Delta$ are not necessary equal and not even of the same order of magnitude, we distinguish these two quantities. The loose rationale for expression (42) is that if $\theta_{\Delta}$ is the interface angle at distance $\Delta$, then the Cox theory, i.e. Eq. (13), reduces to Eq. (42). A more systematic approach amounts to note that if Eq. (42) is used then in region II of our asymptotic analysis, we have from Eqs. (42) and (21)

[TABLE]

The above equation predicts in region II an angle independent of the grid size. If the simulation is sufficiently well resolved and grid independent in the other regions of the domain, then via asymptotic matching of the outer scales with region II, the simulation should also be mesh-independent globally. This most systematic approach differs from the naive interpretation that $\theta_{\Delta}$ is the angle at scale $\Delta$ by a factor $\ln\phi(\theta_{\Delta})=-f(\theta_{\Delta},q)/Q_{i}$.

The approach in the previous publication by two of the authors [27] is identical to the above one except for an approximation where

[TABLE]

valid only in the range defined by Sheng and Zhou [54], that is $q=1$ and $|\cos{\theta}|<0.6$. The approach embodied in Eq. (42) is more general, and can be used to define a procedure to perform mesh-independent simulations of realistic problems. (However, predictive simulations of moving contact lines are not possible if a physical model at the microscopic scale, beyond the Navier-Stokes equations, is not available.)

In such a procedure, one of the main issues is whether contact line physics are obtained by experiment or by a reduction to a microscopic theory such as molecular dynamics or phase field.

To develop a strategy for realistic simulations where experiments would be used, the first step would be to perform a series of simulations at varying $\Delta$ and for the conditions of the experiments. Obviously, simulations and experiments are not performed near the critical ${\rm Ca}$ but below it. Comparing the angles observed in the region where the theory is still valid (region II) would cross-validate the numerics, the Cox’s theory, and the experimental reality. This would in turn fix the parameters $\theta_{in}$ and $r_{m}$ that could be used in simulations. This approach is of course made difficult by the fact that there is no evidence so far that a single pair $\theta_{in},r_{m}$ could predict experiments over a range of different flow configurations. Actually the authors of [9] state that “In the literature we have not found a single pair of experiments in different geometries using precisely the same materials.”. This highlights the difficulty of a predictive simulation approach based on experiments. Moreover, it should be noted that the authors of references [60] and [61] claim precisely the opposite, that no inner angle $\theta_{in}$, even a dynamical one depending on the capillary number, can predict the whole range of experiments they have performed or simulated.

Another approach could be based on mesoscale or nanoscale physics, using molecular dynamics or phase field simulations to predict interface shapes, and again attempting to fit the predicted shapes to those predicted by Navier-Stokes simulations performed using Eq. (42) and a pair $\theta_{in},r_{m}$. The models would have to be simulated up to scales much larger than the nanoscale or the mesoscale so the Cox asymptotics become observable. This could be difficult since in particular the phase field models involve an intrinsic length scale $l_{d}$ (in the notations of ref. [9]) that may be large enough to preclude the appearance of a “Cox region” such as region II inserted between the mesoscale model scale $l_{d}$ and the larger scales of the simulation such as the capillary length $l_{c}$.

All the above considerations must be subject to the proviso that the Cox region II is really observable. There are several circumstances when that cannot be the case. For both withdrawing or advancing contact angles, above a critical ${\rm Ca}$ number, a transition occurs to a liquid film or air film solutions, marking the disappearance of region II. These transitions fix an upper limit on ${\rm Ca}$ for moving contact lines to simply exist. A lower limit for our approach to be applicable is obtained at small capillary numbers. Indeed, the Cox solution is valid only if region II exists, which may be very narrow for small ${\rm Ca}$. Indeed region II transitions directly to region III at small ${\rm Ca}$ for viscous forces of the same order than the hydrostatic pressure gradient $\mu U/\ell^{2}\sim\rho g$ which is equivalent to $\ell~{}\sim l_{c}{\rm Ca}^{-1/2}$ where the capillary length appears. For ${\rm Ca}=10^{-6}$ and the classic capillary length $l_{c}=2.110^{-3}m$ we find $\eta_{I}\sim 2.1\mu$m. This distance may become so small that simulations at this scale are not practical and it is preferable to set the grid size at larger scales at which the contact angle is simply the apparent contact angle or the angle of the meniscus solution.

8 Conclusions

We focus on the problem of forced dewetting transition of a partially wetting plate withdrawn from a reservoir using direct numerical simulations.

We compare numerical solutions in the vicinity of the contact line with the theory of Cox and Voinov and point out the existence of a gauge function $\phi$ that corresponds to one of the integration constants in these theories. We find that Cox-Voinov theory is not well verified at small angles pointing to necessary improvements in the numerical treatment of the dynamic contact line. These improvements may be related to the treatment of boundary conditions or to the treatment of viscosity in mixed cells. They should be the focus of future studies on this topic.

Using the gauge function $\phi$, we provide a numerically-validated approximate correlation for the critical capillary number, ${\rm Ca}_{cr}$, at which the dewetting transition occurs, as a function of the mesh size, $\Delta$, and the numerically-imposed equilibrium contact angle, $\theta_{\Delta}$.

Using the asymptotic hydrodynamic theory of the vicinity of the contact line and matching it to the static theory of menisci, we generalize the hydrodynamic theory of the dewetting transition. We use it to derive an equation for the effect of $\Delta$ and $\theta_{\Delta}$ on the larger-scale regions of the simulation. The critical capillary number is then predicted by an implicit equation for ${\rm Ca}_{cr}$. This equation contains the unknown gauge function $\phi$ that characterizes the contact line dynamics and is akin to a coefficient that determines the amount of slip. This gauge function is specific to the numerical model used. Our numerical simulations allow to quantify this “numerical slip” and confirm that it varies linearly with the grid size $\Delta$.

Of particular interest is the work of Snoeijer in [58] which generalizes the lubrication theory to large angles. The numerical verification of the lubrication equation developed in [58] is topic of our future work.

We suggest a manner in which simulations can be made convergent upon grid refinement, despite the singularity at the contact line. This involves adjusting the numerical contact angle as a function of the grid size. This adaptation of the contact angle involves the numerical gauge function $\phi$ and improves in several ways over the model proposed in [27]. Indeed, it is valid for arbitrary angles and viscosity ratios. However, the microscopic parameters and even more generally the microscopic physics are not known, except for the case treated in this paper of a postulated Navier-slip model. Thus the applicability of the grid-independent model may be limited.

The perspectives of this work are a systematic determination of the gauge function $\phi$ for a range of numerical or physical contact line models used in practice and the verification of the procedure for grid-independent computations in a number of flows. Such grid independent simulations should be performed in conjunction with experiments on contact line dynamics.

Acknowledgements

This research was supported by Electricité de France (EdF), Contract No. 8610-5910129228 (“Advances in Subgrid Models for Subcooled Flow Boiling in Pressurized Water Reactors”) and NSF grant Nos. DMS-1320037 and CBET-1604351. S. A. gratefully acknowledges the support from CNRS for his visit to Institut Jean Le Rond d’Alembert during the preparation of this paper. S. Z. gratefully acknowledges the support from CNRS and MOST (Taiwan) for supporting his visit to NTU during the preparation and writing of this paper.

The authors thank CINES and its team for the grant of computer time and technical assistance on the OCCIGEN supercomputer in the framework of the GENCI allocation x20162b7325. The authors also thank D. Fuster, G. Galliero, C. Josserand, D. Legendre, I. Lunati, A. Malan, H. Tchelepi, M. Abu-Al-Saud and A.-B. Wang for fruitful discussions. S. Z. gratefully acknowledges helpful discussions with Y. Pomeau on contact line modeling, and thanks Jens Eggers for his helpful correspondence about the derivation of equation (37).

Appendix A Matching procedure in the small angle, free-surface case

In the small $\theta_{e}$ case, lubrication theory can be used following ref. [36]. We already expressed the first step of that theory when we obtained the slope in the vicinity of the contact line. It was seen that lubrication theory can be used when $\theta_{e}\ll Ca^{1/3}$. The range of validity of the theory is estimated in ref. [57] and for the advancing case to be ${\rm Ca}/\theta_{e}^{3}<0.05$. This is equivalent for small angles to ${\rm Ca}/G(\theta_{e})\sim\delta<0.5$ which is not too small but larger than the value in the numerical case above.

In region IV, the lubrication equation is well known to be the Airy equation [56]

[TABLE]

and we assume the scaling

[TABLE]

so that the thickness of the film is ${\rm Ca}^{2/3}$ much smaller than its horizontal extent ${\rm Ca}^{1/3}$. Eq. (45) can be solved using Airy functions. The analysis given in [56] is reproduced in B. One finds that for $\zeta\ll{\rm Ca}^{1/3}l_{c}$

[TABLE]

where $\beta$ is a parameter characterizing the Airy function solution of Eq. (45). Matching Eq. (47) with the Cox’s solution Eq. (18) gives

[TABLE]

The matching is performed for $\zeta$ small in region IV variables but large in region II variables, $r_{m}\ll\zeta\ll{\rm Ca}^{1/3}l_{c}$. This is consistent with an upper bound $r_{\rm max}={\rm Ca}^{1/3}l_{c}$ on the validity of the Cox’s solution Eq. (13). At the upper end of region IV, for ${\rm Ca}^{1/3}l_{c}\ll\zeta\ll l_{c}$, one can find the curvature as shown in B

[TABLE]

Eq. (49) is obtained from the Airy function solution, which is parameterized by both $\beta$ and $s_{1}$. The determination of $s_{1}$ is more subtle [36, 38]. It is seen that expression (49) predicts a range of possible curvatures depending on the value of $s_{1}$. The smallest possible curvature obtains when the Airy function assumes its global maximum ${\rm Ai}(s_{\rm max})\simeq 0.53$ for $s_{1}=s_{\rm max}\simeq-1.0$. Thus for a given ${\rm Ca}$ and $\beta$, a minimum curvature is given by

[TABLE]

Indeed if the curvature determined at the lower limit of region III is larger, that is $\kappa^{\rm III}>\kappa_{\rm IV,min}$, it is always possible to match with a solution parameterized by some $s_{1}$ such that ${\rm Ai}(s_{1})<{\rm Ai}(s_{\rm max})$ and ${\rm Ai}(s_{\rm max})$ is small enough. However, if curvature $\kappa^{\rm III}<\kappa_{\rm IV,min}$, the matching is impossible. Thus the critical ${\rm Ca}$ is given by $\kappa^{\rm III,max}=\kappa_{\rm IV,min}$

[TABLE]

We now eliminate $\beta$ between Eq. (48) and Eq. (51) and use the notation $\kappa^{\rm III,max}=\kappa_{\infty}$,

[TABLE]

In the slip-length model case, with $r_{m}$ given by Eq. (19) and $\phi$ given by Eq. (20), we obtain

[TABLE]

This is in agreement with Eggers’s result in [36]. Eggers considers the case of a plate inclined with a small slope $\theta_{p}$ above the horizontal, while we consider a vertical plate. This does not change the nature of the asymptotics and as already pointed out in ref. [36], the theory transposes as well to the case of a vertical plate (although the full numerical solution is then more difficult to obtain). In the small plate angle case $\theta_{p}$ is equal to the curvature $\kappa_{\infty}$ of the meniscus solution, while in our case the curvature is given by Eq. (8). There is agreement between Eq. (53) and Eq. (9) of [36] if we take note of the use of dimensionless variables in [36], while we use dimensional variables, and substitute $\kappa_{\infty}$ for $\theta$ (using ref. [36] notation for $\theta_{p}$).

Appendix B Analysis of the Airy equation

We first outline the Airy function solution of Eq. (45) given in [56]. With the transformation $\eta=(3\mbox{Ca})^{1/3}H$, Eq. 45 becomes

[TABLE]

This equation can be turned into Airy’s equation for a new variable $z(s)$, that is $z^{{}^{\prime\prime}}=sz$, upon the substitutions

[TABLE]

and

[TABLE]

whose solution is

[TABLE]

which implies

[TABLE]

where $\alpha$ and $\beta$ are two arbitrary constants, while $\zeta(s)$ is an antiderivative of $-2^{1/3}z^{-2}(s)$. After integrating explicitly for $\zeta(s)$, Duffy and Wilson [56] conclude that for $\beta>0$ and any given $\alpha$, there is a single branch that behaves as desired for the contact line problem, that it grows monotonically from $H=0$ for $\zeta=0$ and for $s\to\infty$ and has $H\to\infty$ for $s\to s_{1}$ and $\zeta\to\infty$ where $s_{1}$ is a root of

[TABLE]

That single branch cannot go to infinity inside $s\in[s_{1},\infty)$, which implies that $s_{1}$ is the largest among the countably infinite set of roots of Eq. (57). Thus the solution can be characterized either 1) by the arbitrary pair $\alpha,\beta$ ($\beta>0$) or 2) by $\beta$ ($\beta>0$) and $s_{1}$ chosen arbitrarily with $\alpha$ given by $\alpha=-\beta\mbox{Bi}(s_{1})/\mbox{Ai}(s_{1})$ provided there is no larger root of Eq. (57). In what follows, we use the second characterization. It can be verified that

[TABLE]

is an antiderivative of $-2^{1/3}z^{-2}(s)$ and that it satisfies $\zeta=0$ for $s\to\infty$. Thus the solution given parametrically by Eqs. (56) and (58) starts at $H=0$ for $\zeta=0$ and grows monotonically with $H\to\infty$ for $\zeta\to\infty$. Expanding the solution for $s\gg 1$, Duffy and Wilson [56] find that

[TABLE]

where $c=\ln[\pi/(2^{2/3}\beta^{2})]$. Differentiating, and returning to the original variable $\eta$, leads to

[TABLE]

For $\zeta\gg 1$, Eq. (56) reduces to

[TABLE]

Differentiating twice, one obtains the second derivative. As discussed in Sec. 5, it is possible to equate the curvature with the second derivative since 1) ${\rm Ai}^{\prime}(s_{1})=0$, and 2) the curvature is of the order of $l_{c}^{-1}$ and the matching is performed at distances small in the region III variables thus over distances $\zeta\ll l_{c}$. Over such distances, the curvature of the parabola (Eq. (61)) is close to the curvature at its apex. Reverting to the original variable $\eta$, we obtain the curvature $\kappa_{\infty}\simeq\eta^{\prime\prime}(\zeta)$ in the asymptotic range ${\rm Ca}^{-1/3}l_{c}\ll\zeta\ll l_{c}$ as

[TABLE]