S Afkhami

J Buongiorno

A Guion

Stéphane Popinet

Y Saade

R Scardovelli

Stéphane Zaleski
email: zaleski@dalembert.upmc.fr

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

Keywords: Dynamic contact line, Contact angle, Contact line stress singularity, Slip boundary condition, Landau-Levich-Derjaguin film, Forced dewetting, Wetting failure, Cox-Voinov model, Volume-Of-Fluid (VOF), Gerris, Slip length, Navier slip, Partial slip

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

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 [START_REF] Darhuber | Selective dip-coating of chemically micropatterned surfaces[END_REF] and tear films on the cornea [START_REF] Li | Tear film dynamics with evaporation, wetting, and time-dependent flux boundary condition on an eye-shaped domain[END_REF] to micro-layer formation in wall boiling [START_REF] Kim | Detection of liquid-vapor-solid triple contact line in two-phase heat transfer phenomena using high-speed infrared thermometry[END_REF][START_REF] Guion | Dynamics of the liquid microlayer underneath a vapor bubble growing at a heated wall[END_REF] and CO2 sequestration [START_REF] Kimbrel | Experimental characterization of nonwetting phase trapping and implications for geologic CO2 sequestration[END_REF]. 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 [START_REF] Mahady | A volume of fluid method for simulating fluid/fluid interfaces in contact with solid boundaries[END_REF], diffuse-interface models [START_REF] Ding | Wetting condition in diffuse interface simulations of contact line motion[END_REF][START_REF] Yue | Wall energy relaxation in the Cahn-Hilliard model for moving contact lines[END_REF][START_REF] Seppecher | Moving contact lines in the Cahn-Hilliard theory[END_REF][START_REF] Pismen | Disjoining potential and spreading of thin liquid layers in the diffuse-interface model coupled to hydrodynamics[END_REF] and the related issue of evaporation [START_REF] Pomeau | Contact line moving on a solid[END_REF], interface formation models [START_REF] Shikhmurzaev | Moving contact lines in liquid/liquid/solid systems[END_REF][START_REF] Sprittles | Finite element simulation of dynamic wetting flows as an interface formation process[END_REF], and surface roughness [START_REF] Hocking | A moving fluid interface on a rough surface[END_REF]. The reader may find references to other mechanisms in review papers [START_REF] Blake | The physics of moving wetting lines[END_REF][START_REF] Bonn | Wetting and spreading[END_REF][START_REF] Snoeijer | Moving contact lines: Scales, regimes, and dynamical transitions[END_REF][START_REF] Sui | Numerical simulations of flows with moving contact lines[END_REF]. 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. [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF][START_REF] Greenspan | On the motion of a small viscous droplet that wets a surface[END_REF][START_REF] Dussan | On the spreading of liquids on solid surfaces: Static and dynamic contact lines[END_REF][START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF][START_REF] Renardy | Numerical simulation of moving contact line problems using a volume-of-fluid method[END_REF][START_REF] Spelt | A level-set approach for simulations of flows with multiple moving contact lines with hysteresis[END_REF][START_REF] Devauchelle | Forced dewetting on porous media[END_REF][START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF][START_REF] Manservisi | A variational approach to the contact angle dynamics of spreading droplets[END_REF]. 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 [START_REF] Blake | The physics of moving wetting lines[END_REF], which based on experimental measurements lies in the nanometer range [START_REF] Joseph | Direct measurement of the apparent slip length[END_REF][START_REF] Lauga | Microfluidics: The No-Slip Boundary Condition[END_REF]. 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. [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF][START_REF] Moriarty | Effective slip in numerical calculations of moving-contact-line problems[END_REF][START_REF] Weinstein | Scale dependence of contact line computations[END_REF]).

This paper pursues several goals. One is to describe the method used in the numerical framework, Gerris [START_REF] Popinet | The Gerris flow solver[END_REF] 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 [START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF][START_REF] Afkhami | Height functions for applying contact angles to 3D VOF simulations[END_REF], 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 [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF][START_REF] Eggers | Contact line motion for partially wetting fluids[END_REF][START_REF] Chan | Theory of the forced wetting transition[END_REF] 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 [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF][START_REF] Legendre | Comparison between numerical models for the simulation of moving contact lines[END_REF] or [START_REF] Abu-Al-Saud | Level-set method for accurate modeling of two-phase immiscible flow with moving contact lines[END_REF] 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 [START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF][START_REF] Afkhami | Height functions for applying contact angles to 3D VOF simulations[END_REF] 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, 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 [START_REF] Landau | Dragging of a liquid by a moving plate[END_REF][START_REF] Derjaguin | On the thickness of a layer of liquid remaining on the walls of vessels after their emptying, and the theory of the application of photoemulsion after coating on the cine film[END_REF]. 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 [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF][START_REF] Eggers | Contact line motion for partially wetting fluids[END_REF][START_REF] Chan | Theory of the forced wetting transition[END_REF], provided a hydrodynamic prediction, based on the lubrication approximation theory, of the critical capillary number Ca cr . Cox [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF] and Voinov [START_REF] Voinov | Hydrodynamics of wetting[END_REF] described, on the other hand, how the singularity drove a peculiar curved form of the fluid wedge at small 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 [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF][START_REF] Legendre | Comparison between numerical models for the simulation of moving contact lines[END_REF]. 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. [START_REF] Kim | Detection of liquid-vapor-solid triple contact line in two-phase heat transfer phenomena using high-speed infrared thermometry[END_REF], we report the results of our simulations including the method for determining the critical capillary number Ca cr . In Sec. 5, we introduce the main new theoretical concept of this paper, the gauge function φ. 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 [START_REF] Blake | Kinetics of liquid/liquid displacement[END_REF]. 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.

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 ≤ x, y ≤ 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 µ i and ρ i , respectively. The capillary number is defined as

Ca = µ 1 V s /σ,
where µ i is the viscosity of fluid i and σ the surface tension. We set L ≈ 9 l c , where l c is the capillary length l c = σ/[(ρ 1 -ρ 2 )g] with g the gravitational acceleration. The Reynolds number is then defined based on the capillary length as Re = ρ 1 V s l c /µ 1 . Thus Re = Ca Mo -1/2 where the Morton number is

Mo = (ρ 1 -ρ 2 ) gµ 4 1 ρ 2 1 σ 3 .
We define a "gravity-wave-damping" number as

N G = ρ 2 1 V 3 s µ 1 (ρ 1 -ρ 2 ) g .
Indeed, the wavelength of gravity waves traveling at the same speed as the withdrawing tape is

L gw = ρ 1 V 2 s /[(ρ 1 -ρ 2 )g] and it can be connected to N G through N G = ρ 1 L gw V s /µ 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 Re = N

1/2 G Ca -1/2 . (1)
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 µ 1 /µ 2 = 1 and the density ratio ρ 1 /ρ 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 ρ 1 = 5, ρ 2 = 1, g = 16, µ 1 = 1, V s = 1.) Thus from Eq. ( 1) the Reynolds number based on l c varies as

Re = 5 8 Ca -1/2 . ( 2
)
As a result the Reynolds number based on l c increases as the capillary number decreases. Varying the Reynolds number between 0 (Stokes approximation) and Re = 3 has no effect on the results, however increasing Re beyond this value introduces significant inertial effects and interface oscillations that modify the conclusions of our investigations. At small 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 Re = 1. In a final set of simulations (Setup C), we keep ρ 1 /ρ 2 = 5 and Re = 1 but let µ 1 /µ 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 ≈ 3.1 l c . A no-slip boundary condition is prescribed at the substrate (x = 0). Symmetry boundary conditions are imposed on the

Setup N G = ρ 2 1 V 3 s /(µ 1 (ρ 1 -ρ 2 ) g) Re = ρ 1 V s l c /µ 1 1/q = µ 1 /µ 2 ρ 1 /ρ 2 A 25/64 - 1 5 B - 1 1 5 C - 1 50 5
Table 1: Summary of the simulation Setups. When no value is indicated, the corresponding number is computed from the other numbers using Eq. (1).

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 Ca tends to zero. How simulations with a contact angle are performed is described in the next section.

Numerical model

We use Gerris [START_REF] Popinet | The Gerris flow solver[END_REF][START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible Euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF] to numerically solve the Navier-Stokes, continuity and density equations,

ρ (∂ t u + u • ∇u) = -∇p + ∇ • µ ∇u + ∇u + σκδ s n + ρg, (3)
∇ • u = 0, (4)
∂ t ρ + u • ∇ρ = 0, (5)
respectively. Here, u is the velocity field, p the pressure,

ρ = ρ 1 χ + ρ 2 (1 -χ), µ = µ 1 χ + µ 2 (1 -χ)
, κ is the interface curvature, n the normal to the interface (pointing from fluid 1 to fluid 2), δ s the delta function centered at the interface, ρg = -ρg ŷ the body force due to gravity and ŷ the unit vector in the ydirection, and χ (= 1 in fluid 1 and 0 in fluid 2) the characteristic function, where δ s n = ∇χ. Note that Eq. ( 5) is equivalent to

∂ t χ + u • ∇χ = 0, (6)
which is solved using the Volume-Of-Fluid (VOF) interface capturing method [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF][START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF][START_REF] Popinet | Numerical models of surface tension[END_REF]. The Continuous Surface Force (CSF) method is used for the surface tension force with curvature computed using the Height-Function method [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | Numerical models of surface tension[END_REF]. Viscous forces are implemented using a partially implicit method described in [START_REF] Lagrée | The granular column collapse as a continuum: validity of a two-dimensional Navier-Stokes model with a µ(I)rheology[END_REF].

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,

B A B' A' i-1/2,j i+1/2,j u u i,j i,j i+1,j
Figure 2: Advection of the interface along the x-direction: standard "Lagrangian Explicit" or "CIAM" method with the cell centered at grid point i, j being advected and expanded/compressed by the flow, in this case it is compressed since ∂u/∂x < 0; before advection (left) and after advection (right).

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 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. [START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]) or the "mixed Youngs-centered" (MYC) method (see e.g. [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]), 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 [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]) from the knowledge of n and C i,j . Thus the equation of the segment is written

n x x + n y y = α, (7)
where the scalar α characterizes the position of the interface. The knowledge of n and α 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 [START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]) 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 [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Tryggvason | Direct Numerical Simulations of Gas-Liquid Multiphase Flows[END_REF]. 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 i+1,j i,j i,j

Figure 3: Advection of the interface along the x-direction with the cell subdivided in m = 4 equal bands: with "banding" the effect of shearing, that is a ∂u/∂y derivative, is taken into account more accurately than in the standard method; before advection (left) and after advection (right).

the "banded" advection [START_REF] Popinet | Rotation of a straight interface[END_REF]. In this test case, a straight interface is advected by a pure shear flow, at a shear rate of 1, in a 1 × 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

h i = j2 j=j1 C i,j
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,j1 = 1) and the top cell is empty (C i,j2 = 0), and there is a single interface in the column, the height h i approaches the exact interface height to second order [START_REF] Bornia | On the properties and limitations of the height function method in two-dimensional cartesian geometry[END_REF]. 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 [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF][START_REF] Popinet | Numerical models of surface tension[END_REF]).

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 [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF]. 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 θ ∆ as n x = -sin(θ ∆ ), n y = cos(θ ∆ ). The value of α 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(θ ∆ )∆ is used, where ∆ 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. [START_REF] Mahady | A volume of fluid method for simulating fluid/fluid interfaces in contact with solid boundaries[END_REF].

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 [START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF], 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 κ∆ small enough. The critical value of κ∆ was estimated independently by one of the authors through tests on a large number of random circles showing that the minimum value of κ∆ 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 ∆ 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 ∆ 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.

(a) (b) (c) (d)

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 τ → ∞, where the nondimensional time is τ = 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, ∆/l c , and the imposed contact angle, θ ∆ . 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 ∆/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 < ∆/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 ∆ l c L, where is the thickness of the film at the transition, of the order = 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 ∆ 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 ∞ -h 0 )/l c , where h ∞ = h CL (τ → ∞), as a function of the nondimensional mesh size, ∆/l c , for various capillary numbers and for various contact angles, 50 • ≤ θ ∆ ≤ 90 • . The results show the grid sizes, ∆/l c , where a stationary meniscus forms and no film is deposited on the substrate, for the range of considered θ ∆ .

More importantly, the results show that the computed height of the contact line is a function of ∆/l c and that for small enough θ ∆ no steady state menisci can be attained when ∆/l c ≤ 0.014. This lower limit of θ ∆ , for which steady state contact lines can be achieved, gets larger as 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 Ca is increased or θ ∆ 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 Ca cr , we propose a procedure as follows. Fig. 9(a) shows contact line heights as a function of time for various values of Ca, when θ ∆ = 90 • and ∆/l c = 0.007. As shown, for sufficiently small Ca, a stationary meniscus can be reached while for large 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τ |, does not go to zero as a function of τ . This critical capillary number Ca cr is depicted in red in Figs. 9(a) and (b). Using the procedure above, we can therefore determine 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, τ , for two mesh sizes, ∆/l c , when varying the wall contact angle θ ∆ , for a fixed Ca. As shown, the transition from a stationary meniscus not only depends on θ ∆ , but also on ∆/l c . As illustrated, the transition occurs at a larger θ ∆ for smaller ∆/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 Ca = 0.043, θ ∆ = 90 • , and ∆/l c = 0.014 (for this set of parameters, Ca cr = 0.52).

θ ∆ =50 o θ ∆ =90 o (h ∞ -h 0 )/l c ∆/l c 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
0.01 0.02 0.03 0.04 0.05 0.06 0.07

θ ∆ =54 o θ ∆ =90 o (h ∞ -h 0 )/l c ∆/l c (a) (b) 0.4 0.6 0.8 1 1.2 1.4
0.01 0.02 0.03 0.04 0.05 0.06 0.07 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 Ca = 0.048, θ ∆ = 60 • , and ∆/l c = 0.014 (for this set of parameters, Ca cr = 0.024). The figure shows a typical evolution of the interface and the transition to the film. Figs. 12(d)-( 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.

θ ∆ =66 o θ ∆ =90 o (h ∞ -h 0 )/l c ∆/
θ ∆ =78 o θ ∆ =90 o (h ∞ -h 0 )/l c ∆/l c (c) (d)
θ ∆ =80 o θ ∆ =50 o (h CL -h 0 )/l c τ (a) (b)
We now turn our attention to computing the critical capillary number, Ca cr , of the dewetting transition, as a function of ∆, θ ∆ , and q. We carry out a series of simulations when incrementing the capillary number, holding ∆ and θ ∆ fixed, until we find 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 ∆/l c ≈ 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 Ca cr as a function of ∆/l c for a range of θ ∆ 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 θ ∆ . The error bars represent the increment in the capillary number in the procedure of finding 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 Ca cr ; in general, lowering the viscosity of fluid 2 results in an increase in 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 Ca cr and give the scaling of it with ∆/l c .

Hydrodynamic theories of the dynamic contact line and the dewetting transition

Generalities about the asymptotic description

To interpret our numerical results, we will use a theoretical framework that extends the work in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF][START_REF] Chan | Theory of the forced wetting transition[END_REF]. The theory in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF][START_REF] Chan | Theory of the forced wetting transition[END_REF] 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 [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF]. 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 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 λ. We make no assumptions about region I except that it is limited in length scales to a microscopic scale λ that remains much smaller than all other relevant length scales.

θ ∆ =110 o θ ∆ =90 o θ ∆ =70 o θ ∆ =60 o θ ∆ =50 o θ ∆ =40 o θ ∆ =30 o θ ∆ =20 o
θ ∆ =110 o θ ∆ =90 o θ ∆ =70 o θ ∆ =60 o θ ∆ =50 o θ ∆ =40 o θ ∆ =30 o θ ∆ =20 o Ca cr ∆/l c (b) 0.001 0.01 0.1 0.01 0.1 θ ∆ =110 o θ ∆ =90 o θ ∆ =40 o θ ∆ =30 o θ ∆ =20 o Ca cr ∆/l c (c)
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 [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF]. This region starts at lengths λ and ends at lengths 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 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 [START_REF] Landau | Fluid Mechanics[END_REF]. We define the "apparent contact angle" θ a as the angle seen in meniscus variables that is on scales such that y ∼ l c . The meniscus solution quadrature gives the curvature

κ III = l -1 c 2 -2 sin θ a (8)
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 θ 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 θ 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 θ a = θ e . This condition closes the problem and gives in particular the height of the contact line h cl = l c √ 2 -2 sin θ e . As the capillary number increases, the curvature and the bending increase in region II. The apparent contact angle decreases below θ 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 [START_REF] Derjaguin | On the thickness of a layer of liquid remaining on the walls of vessels after their emptying, and the theory of the application of photoemulsion after coating on the cine film[END_REF] and Landau and Levich [START_REF] Landau | Dragging of a liquid by a moving plate[END_REF] 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 [START_REF] Landau | Dragging of a liquid by a moving plate[END_REF] have shown that in region IV and at small capillary number, the film thickness behaves as h ∼ 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 [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF][START_REF] Chan | Theory of the forced wetting transition[END_REF], 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 is small enough. This parameter is not simply Ca as the Landau-Levich film scaling would seem to imply, but ∼ Ca/G(θ e ) where G(θ e ) is the function introduced in [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF] and defined below. It can also be shown that at the transition, an equivalent estimate of the small parameter is ∼ 1/ ln(l c /r m ). This value of is not exceedingly small, which has consequences which will be spelled out below.

The theory of Cox

We thus focus on the analysis of region II. In this region, the wedge solution of Huh & Scriven is assumed [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF]. It is then possible to use it as done by Cox [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF] 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

G[θ(r)] = G(θ e ) -Ca ln(r/λ c ) -Ca Q i f (θ e , q) + o(Ca), (9)
where λ c is a characteristic scale for microscopic effects, θ e is the equilibrium angle in agreement with our assumption above for the Ca → 0 limit, f and G are defined below, and q = µ 2 /µ 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 Ca, higher-order terms exist whose form is however unknown to the authors. Finally the functions f and G are given by

G(θ) = θ 0 dθ f (θ , q) , (10)
and

f (θ, q) = 2 sin θ{q 2 (θ 2 -sin 2 θ) + 2q[θ(π -θ) + sin 2 θ] + [(π -θ) 2 -sin 2 θ]} q(θ 2 -sin 2 θ)[(π -θ) + cos θ sin θ] + [(π -θ) 2 -sin 2 θ](θ -cos θ sin θ) .
In some specific asymptotic limits, special forms of G(θ) can be used. For θ 1 and q = 0 as in thin liquid film in air situations, it can be shown that

G(θ) ≈ θ 3 /9. (11)
Also, Sheng and Zhou [START_REF] Sheng | Immiscible-fluid displacement: Contact-line dynamics and the velocity-dependent capillary pressure[END_REF] show that G(θ) -G(θ e ) ≈ (cos θ e -cos θ)/5.63 when q = 1 and | cos θ| < 0.6. In this work, we account for G(θ) using Eq. ( 10) directly.

A more convenient notation is to introduce a "microscopic length scale"

r m = λ c exp[f (θ e , q)/Q i ], , (12)
as in Voinov [START_REF] Voinov | Hydrodynamics of wetting[END_REF][START_REF] Voinov | Wetting: Inverse dynamic problem and equations for microscopic parameters[END_REF] so that in the general case

G(θ) = G(θ e ) -Ca ln(r/r m ), (13)
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 θ e and r m . Eq. ( 9) is very attractive because it provides a universal description giving the dependence of θ(r) on 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 λ c ∼ ∆.

The microscopic length r m can be obtained in a number of special cases, or by dimensional analysis. Eggers [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF] has assumed several slip length models. For the Navier-slip model, and using our notations above,

v| x=0 -V s = λ ∂v ∂x x=0 , ( 14
)
where λ is the slip length. Using lubrication theory valid for small slopes or angles, one obtains the equation [START_REF] Duffy | A third-order differential equation arising in thin-film flows and relevant to tanner's law[END_REF]

η = 3Ca η 2 + 3λη , (15)
where we introduced a variable ζ that measures distance away from the contact line along the direction of the solid surface. We also note that η(ζ) is the local thickness of the film and we use the notation η = dη/dζ, η = d 2 η/dζ 2 , etc.. A solution of Eq. ( 15) can be obtained [START_REF] Eggers | Existence of receding and advancing contact lines[END_REF] in the form

η (ζ) ∼ θ e - 3Ca θ 2 e ln ζ r m , (16)
where r m = 3λ e θ e + O(Ca).

(Notice that equation ( 33) of [START_REF] Eggers | Existence of receding and advancing contact lines[END_REF] misses the above factor of e because of an incorrect derivation from equation ( 32) of [START_REF] Eggers | Existence of receding and advancing contact lines[END_REF]). This can be compared to the solution obtained by Cox in [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF]. From Eq. ( 11) and Eq. ( 13), using the fact that for small angles r ∼ ζ = h CL -y and θ ∼ η (ζ), one obtains

η (ζ) 3 ∼ θ 3 e -9Ca ln(ζ/r m ). ( 18
)
For Ca θ 3 e , the above yields again Eq. ( 16). Thus the length scale r m can be written as

r m = λ/φ, ( 19
)
where φ is a gauge function depending on θ e equal, at first order, to

φ(θ e ) = e θ e 3 . (20)
Different slip-length models give different gauge functions φ. 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 φ having a dependence on θ 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

r m = ∆/φ num , (21)
where the gauge function φ num may depend on q, θ 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 φ on the viscosity ratio q implicit. Comparing expressions ( 19) and ( 21) leads to the relation

λ = φ∆/φ num , (22)
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

G(θ e ) -G[θ(r)] = Ca ln(r/∆) + Ca ln φ. ( 23
)
This form is verified numerically in the next Section.

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 • ≤ θ ∆ ≤ 110 • and 0.004 ≤ ∆/l c ≤ 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, θ ∆ , 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. ζ/l c = (h CL -y)/l c , (see Fig. 1 for the illustration of the distance ζ). Figs. 14 to 22 depict the curvature plots (black symbols) along with the angle (red symbols) θ a = θ(ζ/l c ) that the interface makes with the substrate as a function of ζ/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 θ ∆ = 110 • ; for example, for ∆ = 0.004, curvature is tending to zero around -0.75l c . However, the maximum curvature does not reach the asymptotic value of √ 2l -1 c when the asymptotic parameter = Ca/G(θ e ) is not small enough. As shown, for small angles, where the smallness of is guaranteed, the maximum curvature approaches √ 2l -1 c . 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 [START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF] on Figs. 14 to 22. We plot the computed (black symbols) G(θ e = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆), using the computed interface slope, θ(r), at the distance r = x 2 i + y 2 i , 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(θ e = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. ( 23) using the best fitted value of φ (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(θ e = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) should be unity. In fact, there is a deviation from unity which is largest at small angles (15 • and 30 • ) 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 φ in the range where Cox's theory is expected to be asymptotically valid ∆ r r max . The expected upper limit r max = Ca 1/3 l c will be derived below. We consider the results for 0.004 ≤ ∆/l c ≤ 0.016 and show that the fitted values of φ are similar and thus approximately independent of ∆ when Ca is maintained at criticality. We extend the discussion on these results next. 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 θ ∆ = 110 • , 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 θ ∆ ≤ 30 • , even the smallest grid size appears to be just sufficient, illustrating the challenging computations at smaller angles.

General case

In the arbitrary angle case or when q does not vanish, one cannot use lubrication theory as done in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF]. 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 (η = 0 when η = 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 • ≤ θ ∆ ≤ 110 • and 0.004 ≤ ∆/l c ≤ 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 [START_REF] Snoeijer | Free-surface flows with large slope: Beyond lubrication theory[END_REF]. A mapping from the equation valid for large slopes in [START_REF] Snoeijer | Free-surface flows with large slope: Beyond lubrication theory[END_REF] to the lubrication approximation for small slopes has been performed in [START_REF] Qin | Asymptotic theory of fluid entrainment in dip coating[END_REF], allowing to reuse the result in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF]; 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 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 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 ζ ∼ ζ I with ζ I r m . The bending occurs over a large distance η 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 η I (paralell to the wall) and ζ I (perpendicular to the wall). The meaning of all the geometrical quantities s, r, η, ζ, θ and θ e is summarized on Fig. 24. As a result of the thinness, we can approximate r ∼ ζ in that region. Writing the small slope assumption θ = 0 at

and Ca cr is small if 1/ ln(ζ I /r m ) is small. The small parameter in the expansion is

= 1/ ln(ζ I /r m ) = Ca cr /G(θ e ). (25)
This determines the scale of Ca and we now turn to determining the thickness η 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 η I ∼ Ca 2/3 l c solutions near the inflection point that match the meniscus, there are also solutions with an asymptotically flat interface, η → η ∞ as ζ → ∞. These solutions have η ∞ ∼ 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 η = Constant in the generalized lubrication equation or in the large slope equation in [START_REF] Snoeijer | Free-surface flows with large slope: Beyond lubrication theory[END_REF]. 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 Ca is small, it is possible to deduce that the curvature is still small in most of region II thus ds/dr ∼ 1 where s is the curvilinear abscissa along the interface. Then

η I ∼ ζ I rm sin G -1 G(θ e ) -G(θ e ) ln r r m dr, (26)
where we have used the fact that r I ∼ ζ I , dη/ds = sin θ(s) and Ca Ca cr . After some work it can be shown that this integral behaves as

η I ∼ 1/3 G(θ e ) 1/3 exp 1 r m , (27)
while by definition

ζ I ∼ exp 1 r m . (28)
We thus have η I ∼ Ca 1/3 ζ I as expected for the inflection point in the thin film dewetting transition. In region IV, the length scales are thus η I across the film and ζ I parallel to the film. The curvature thus scales as

η ∼ η I /ζ 2 I ∼ 1/3 G(θ e ) 1/3 exp - 1 r -1 m . (29)
To match region IV and region III, one should have the curvature η of same order as the curvature l -1 c of the static meniscus. Writing η ∼ l -1 c one obtains

1/3 G(θ e ) 1/3 exp - 1 l c r m ∼ 1. ( 30
)
Identifying as above

C(q)Ca 1/3 cr l c r m exp - G(θ e ) Ca cr = 1, (31)
where we have introduced a constant C(q) = 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 [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF], as shown in Appendix A, to obtain

C(0) = 3 1/3 2 -1/3 πeAi 2 (s max )κ ∞ l c . ( 32
)
Or equivalently (see A.8)

3 1/3 2 -1/3 Ca 1/3 cr πeAi 2 (s max )r m κ ∞ exp - G(θ e ) Ca cr = 1. (33)
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 µ 1 V s /η(y) while the stress in the outer fluid scales as the much smaller µ 2 V s /y, realizing an approximation of the free surface condition. We note that expression ( 33) is equivalent at small angles θ e to the one obtained in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF]. It is also equivalent to the expression obtained in [START_REF] Qin | Asymptotic theory of fluid entrainment in dip coating[END_REF] at all angles apart for the fact that we make no hypothesis on the behavior of r m at large angles, unlike [START_REF] Qin | Asymptotic theory of fluid entrainment in dip coating[END_REF] and [START_REF] Snoeijer | Free-surface flows with large slope: Beyond lubrication theory[END_REF]. 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 [START_REF] Huh | Hydrodynamic model of steady movement of a solid/liquid/fluid contact line[END_REF], but for different boundary conditions.

We also note that at the transition the condition (30) yields η I ∼ Ca 2/3 l c , and that we used essentially the same matching argument as Derjaguin [START_REF] Derjaguin | On the thickness of a layer of liquid remaining on the walls of vessels after their emptying, and the theory of the application of photoemulsion after coating on the cine film[END_REF] and Landau and Levich [START_REF] Landau | Dragging of a liquid by a moving plate[END_REF]. This can also be compared to the prediction of Voinov [START_REF] Voinov | Wetting: Inverse dynamic problem and equations for microscopic parameters[END_REF] 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 [START_REF] Greenspan | On the motion of a small viscous droplet that wets a surface[END_REF] for the microscopic length scale r m in the critical matching expression [START_REF] Popinet | The Gerris flow solver[END_REF] to obtain the central relation of our theory

C(q)φ(θ ∆ , q)Ca 1/3 cr l c ∆ exp - G(θ ∆ ) Ca cr = 1, (34)
where we identified θ ∆ with θ e . It is interesting to investigate explicit expressions for Ca cr obtained by an asymptotic analysis of Eq. [START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF]. To perform this analysis, we let δ = 1/ ln(l c /∆). This small parameter is of the same order as above but unlike is explicitly defined using the parameters of the problem. To give an order of magnitude for δ and , if we use, as in the above simulation, a grid size ∆/l c = 1/2 8 , we obtain δ ≈ 0.18. Thus our small parameters is small but not exceedingly so. It is useful to introduce the auxiliary parameter μ defined as μ = -1/ ln[C(q)φ(θ ∆ , q)G(θ ∆ , q)]. The parameter μ will be small for transcendentally small θ ∆ . Otherwise if μ is considered order 1, solutions of Eq. ( 34) expand as Ca cr = δG(θ ∆ , q) 1 -δ ln δ -δ μ-1 + O(δ 2 ln δ) + O(δ 2 μ-2 ) .

At small angles, with the expected scaling φ(θ ∆ ) ∼ θ ∆ , one gets μ ∼ -1/(4 ln θ ∆ ). Moreover the higher order terms in Eq. ( 35) are small if μ δ which is verified even for the smallest θ ∆ and typical values of l c /∆ in our calculations. For Setup B, with θ ∆ = 110 • , ∆/l c = 1/2 8 , q = 1, and using the best fit estimate φ 3, the first order of expansion of Eq. [START_REF] Afkhami | Height functions for applying contact angles to 3D VOF simulations[END_REF] gives Ca cr Gδ 0.0558, using all the terms in expansion [START_REF] Afkhami | Height functions for applying contact angles to 3D VOF simulations[END_REF] gives Ca cr 0.0626, while a full iterative solution of Eq. ( 34) yields Ca cr 0.0637. For Setup C, in the case θ ∆ = 110 • , ∆/l c = 1/2 8 , q = 0.02, and using the best fit estimate φ 3.5, the first order gives Ca cr Gδ 0.127, using all the terms in expansion [START_REF] Afkhami | Height functions for applying contact angles to 3D VOF simulations[END_REF] gives Ca cr 0.169, while a full iterative solution of Eq. ( 34) yields Ca cr 0.132. Thus in a wide range of parameters, the approximation

Ca cr = G(θ ∆ , q) ln(l c /∆) , (36)
is correct within 15% accuracy.

Comparison with numerics

In our simulations, we measure Ca cr and specify ∆ and θ ∆ . For the two special cases described above, we compare the value of Ca cr from full simulations to the solutions of the central relation, Eq. [START_REF] Afkhami | Height functions for applying contact angles to 2D VOF simulations[END_REF]. The simulations values are in remarkable agreement with the above solutions. For setup B, with θ ∆ = 110 • , ∆/l c = 1/2 8 , q = 1, we find numerically Ca cr = 0.0555 compared to the 1) The microscopic angle is discontinuous, that is there are two angles, the advancing angle θ a and the receding angle θ r , so that θ r < θ e < θ a and that for small Ca, θ in ∼ θ a for an advancing contact line and θ in ∼ θ r for a receding contact line. The conclusion is that the asymptotic analysis is unchanged with θ e replaced by either θ r or θ a .

θ ∆ =110 o θ ∆ =90 o θ ∆ =70 o θ ∆ =60 o θ ∆ =50 o θ ∆ =40 o θ ∆ =30 o θ ∆ =20 o φ ∆/l c 0 0.5 1
θ ∆ =110 o θ ∆ =90 o θ ∆ =70 o θ ∆ =60 o θ ∆ =50 o θ ∆ =40 o θ ∆ =30 o θ ∆ =20 o θ ∆ =15 o φ ∆/l c 0 0.5 1 1.5 2 2.5 3 3.5
θ ∆ =110 o θ ∆ =90 o θ ∆ =40 o θ ∆ =30 o θ ∆ =20 o θ ∆ =15 o φ ∆/l c
θ ∆ =110 o θ ∆ =90 o θ ∆ =30 o φ ∆/l c (e) (f)
2) The microscopic angle is continuous so that at first order

θ in ∼ θ e + Caf Ca , (40)
with f Ca a constant. In this case, and the small Ca assumption, Eq. ( 39) can be rewritten at first order in Ca in the form of Eq. ( 9)

G[θ(r)] = G(θ e ) -Ca ln(r/λ) -Ca Q i f (θ e , q) + o(Ca) (41)
the change in microscopic angle with Ca being absorbed into the new integration constant Q i . We note that the relation ( 40) is in agreement with the Molecular Kinetics Theory (MKT) of Blake and Haynes [START_REF] Blake | Kinetics of liquid/liquid displacement[END_REF], 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 Ca down to 10 -6 , although the dewetting transition itself was studied only for Ca > 10 -4 . For very small 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 = ∆/V s as if the interface was pinned with spatial periodicity equal to the grid size ∆. However, no such irregularity is observed for Ca = 10 -4 . Thus we may infer that at least in the range of parameters used in this study 10 -3 Ca 0.1, the numerical solution is a continuous function of the contact line position.

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 [START_REF] Legendre | Comparison between numerical models for the simulation of moving contact lines[END_REF] 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 [START_REF] Legendre | Comparison between numerical models for the simulation of moving contact lines[END_REF] is to specify numerically the angle

G(θ ∆ ) = G(θ in ) -Ca ln(∆/r m ) (42)
for some θ 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 θ ∆ is the numerically imposed contact angle while θ in is an experimentally measured or determined quantity whose meaning is the contact angle at the microscopic scale r m . Since r m and ∆ 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 θ ∆ is the interface angle at distance ∆, 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)

G(θ) = G(θ in ) -Ca ln rφ(θ ∆ ) r m . (43)
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 θ ∆ is the angle at scale ∆ by a factor ln φ(θ ∆ ) = -f (θ ∆ , q)/Q i . The approach in the previous publication by two of the authors [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF] is identical to the above one except for an approximation where

G(θ) -G(θ in ) [cos θ -cos θ in ]/5.63 (44)
valid only in the range defined by Sheng and Zhou [START_REF] Sheng | Immiscible-fluid displacement: Contact-line dynamics and the velocity-dependent capillary pressure[END_REF], that is q = 1 and | cos θ| < 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 ∆ and for the conditions of the experiments. Obviously, simulations and experiments are not performed near the critical Ca but below it. Comparing the angles observed in the region where the theory is still valid (region II) would crossvalidate the numerics, the Cox's theory, and the experimental reality. This would in turn fix the parameters θ 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 θ in , r m could predict experiments over a range of different flow configurations. Actually the authors of [START_REF] Yue | Wall energy relaxation in the Cahn-Hilliard model for moving contact lines[END_REF] 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 [START_REF] Blake | Experimental evidence of nonlocal hydrodynamic influence on the dynamic contact angle[END_REF] and [START_REF] Wilson | Nonlocal hydrodynamic influence on the dynamic contact angle: Slip models versus experiment[END_REF] claim precisely the opposite, that no inner angle θ 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 θ 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. [START_REF] Yue | Wall energy relaxation in the Cahn-Hilliard model for moving contact lines[END_REF]) 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 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 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 Ca. Indeed region II transitions directly to region III at small Ca for viscous forces of the same order than the hydrostatic pressure gradient µU/ 2 ∼ ρg which is equivalent to ∼ l c Ca -1/2 where the capillary length appears. For Ca = 10 -6 and the classic capillary length l c = 2.110 -3 m we find η I ∼ 2.1µ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.

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 φ 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 φ, we provide a numerically-validated approximate correlation for the critical capillary number, Ca cr , at which the dewetting transition occurs, as a function of the mesh size, ∆, and the numerically-imposed equilibrium contact angle, θ ∆ .

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 ∆ and θ ∆ on the larger-scale regions of the simulation. The critical capillary number is then predicted by an implicit equation for Ca cr . This equation contains the unknown gauge function φ 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 ∆. Of particular interest is the work of Snoeijer in [START_REF] Snoeijer | Free-surface flows with large slope: Beyond lubrication theory[END_REF] which generalizes the lubrication theory to large angles. The numerical verification of the lubrication equation developed in [START_REF] Snoeijer | Free-surface flows with large slope: Beyond lubrication theory[END_REF] 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 φ and improves in several ways over the model proposed in [START_REF] Afkhami | A mesh-dependent model for applying dynamic contact angles to VOF simulations[END_REF]. 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 φ 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. value of s 1 . The smallest possible curvature obtains when the Airy function assumes its global maximum Ai(s max ) 0.53 for s 1 = s max -1.0. Thus for a given Ca and β, a minimum curvature is given by κ IV,min = (3Ca) 1/3 2 1/6 β πAi(s max ) 2 .

(A.6) Indeed if the curvature determined at the lower limit of region III is larger, that is κ III > κ IV,min , it is always possible to match with a solution parameterized by some s 1 such that Ai(s 1 ) < Ai(s max ) and Ai(s max ) is small enough. However, if curvature κ III < κ IV,min , the matching is impossible. Thus the critical Ca is given by κ III,max = κ IV,min (3Ca cr ) 1/3 2 1/6 β πAi(s max ) In the slip-length model case, with r m given by Eq. ( 19) and φ given by Eq. ( 20), we obtain θ e 18 1/3 πAi 2 (s max )

Ca 1/3 cr λ exp - θ 3 e 9Ca cr = κ ∞ (A.9)
This is in agreement with Eggers's result in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF]. Eggers considers the case of a plate inclined with a small slope θ 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. [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF], 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 θ p is equal to the curvature κ ∞ of the meniscus solution, while in our case the curvature is given by Eq. [START_REF] Ding | Wetting condition in diffuse interface simulations of contact line motion[END_REF]. There is agreement between Eq. (A.9) and Eq. ( 9) of [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF] if we take note of the use of dimensionless variables in [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF], while we use dimensional variables, and substitute κ ∞ for θ (using ref. [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF] notation for θ p ).

Appendix B. Analysis of the Airy equation

We first outline the Airy function solution of Eq. (A.1) given in [START_REF] Duffy | A third-order differential equation arising in thin-film flows and relevant to tanner's law[END_REF]. With the transformation η = (3Ca) 1/3 H, Eq. A.1 becomes

H (ζ) = 1 H 2 (ζ)
. 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) Ai (s 1 ) = 0, and 2) the curvature is of the order of l -1 c and the matching is performed at distances small in the region III variables thus over distances ζ l c . Over such

Fluid 2 ,Figure 1 :

 21 Figure 1: Schematic illustrating the contact angles θe and θ(r), corresponding to length scales λ and r, respectively. h CL represents the contact line height and h 0 is the unperturbed initial interface profile. ζ = h CL -y is the vertical distance of the interface at y from the contact line. η(y) is the horizontal distance of the interface from the solid substrate.

Figure 4 :

 4 Figure 4: Evolution of the VOF interface in a shear flow at a shear rate of 1. (a) t = 0, (b) t = 1, (c) t = 2, and (d) t = 5; for m = 1 ( ) and m = 4 ( ). Note that in (d), the VOF reconstructions for m = 1 ( ) in the center of the domain become identical to the black lines of the grid and hidden by them.

Figure 5 :

 5 Figure 5: Volume fraction error as a function of time; for m = 1 ( ) and m = 4 ( ).

Figure 6 :

 6 Figure 6: Reconstruction of the interface and height functions in the vicinity of the contact line. (a) Reconstruction in cell i, j. (b) Linear extension to cell i, j -1. (c) Standard Height-Function method in columns i and i + 1 together with extrapolation of the height to cell i -1. (d) Same as (c) but with widths in horizontal segments j -1, j and j + 1.

Figure 7 :

 7 Figure 7: Contact line height as a function of time for θ ∆ = 90 • ( ) and θ ∆ = 60 • ( ) at various ∆/lc. For ∆/lc = 0.014 and θ ∆ = 60 • , the contact line elevation continues increasing as τ increases. In this figure, Ca = 0.03.

  0.02 0.03 0.04 0.05 0.06 0.07

  0.02 0.03 0.04 0.05 0.06 0.07

Figure 8 :Figure 9 :

 89 Figure 8: Nondimensional stationary contact line height, (h∞ -h 0 )/lc, as a function of the nondimensional mesh size, ∆/lc, for (a) Ca = 0.01, (b) Ca = 0.02, (c) Ca = 0.03, and (d) Ca = 0.04, for various contact angles, θ ∆ ; the contact angle difference between each set is 2 • . For (Ca,θ ∆ )=(0.02,54 • ), (0.03,66 • ), and (0.04,78 • ), no steady state menisci can be attained when ∆/lc ≤ 0.014.

Figure 10 :

 10 Figure 10: Contact line height as function of time for 50 • ≤ θ ∆ ≤ 80 • and (a) ∆/lc = 0.028 and (b) ∆/lc = 0.014; the increment in θ ∆ is 2 • . Ca = 0.03

Figure 11 :

 11 Figure 11: (a) A stationary meniscus forms when Ca < Cacr as τ → ∞. The inset shows the magnified flow field and the pressure distribution. Ca = 0.043, θ ∆ = 90 • , and ∆/lc = 0.014. (b) 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. The pressure colors show the maximum (dark red) and minimum (dark blue) of the pressure distribution.

  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.

Figure 12 :

 12 Figure 12: Time evolution of the interface for Ca > Cacr at τ = 4 (a), 4.8 (b), 5.9 (c), 6.8 (d), 7.9 (e), 8.7 (f). The insets show the magnified flow field and the pressure distribution. (a)-(c) Right panels show 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. Ca = 0.048, θ ∆ = 60 • , and ∆/lc = 0.014. The pressure colors show the maximum (dark red) and minimum (dark blue) of the pressure distribution.

Figure 13 :

 13 Figure 13: Cacr as a function of ∆/lc for a range of θ ∆ for Setups A (a), B (b), and C (c). Symbols present the numerical results and the solid lines are drawn to guide the eye. The error bars represent the increment in the capillary number in the procedure of finding Cacr.

Figure 14 :

 14 Figure 14: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 0.64 (Setup A), (c) 0.85 (Setup B), and (e) 0.77 (Setup C) (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. θ ∆ = 15 • and ∆/lc = 0.0035 (Setup A) and 0.004 (Setup B and C), at Ca = 0.0009, τ = 1.52 (Setup A), Ca = 0.0009, τ = 1.2 (Setup B), and Ca = 0.00095, τ = 1.52 (Setup C).

Figure 15 :

 15 Figure 15: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 2.15 (Setup B) and (c) 1.51 (Setup C) (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. θ ∆ = 30 • and ∆/lc = 0.004 (Setup B and C), at Ca = 0.004, τ = 5.37 (Setup B) and Ca = 0.0049, τ = 5.95 (Setup C).

Figure 16 :

 16 Figure 16: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 2.71 (Setup A), (c) 2.01 (Setup B), and (e) 1.47 (Setup C) (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. θ ∆ = 30 • and ∆/lc = 0.007 (Setup A) and 0.008 (Setup B and C), at Ca = 0.004, τ = 1.52 (Setup A), Ca = 0.005, τ = 3.18 (Setup B), and Ca = 0.0065, τ = 3.04 (Setup C).

Figure 17 :

 17 Figure 17: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 3.08 (Setup A), (c) 2.42 (Setup B), and (e) 1.52 (Setup C) (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. θ ∆ = 40 • and ∆/lc = 0.007 (Setup A) and 0.008 (Setup B and C), at Ca = 0.008, τ = 5.55 (Setup A), Ca = 0.009, τ = 5.52 (Setup B), and Ca = 0.013, τ = 6.27 (Setup C).

Figure 18 :

 18 Figure 18: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 3.14 (Setup A), (c) 2.42 (Setup B), (e) 2.78 (Setup B)(black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. θ ∆ = 60 • and ∆/lc = 0.014 (a-b) (Setup A), 0.016 (c-d) (Setup B), and 0.008 (e-f) (Setup B) at Ca = 0.024, τ = 8.85 (a-b) (Setup A), Ca = 0.024, τ = 10.48 (c-d) (Setup B), and Ca = 0.022, τ = 11.12 (e-f) (Setup B).

Figure 19 :

 19 Figure 19: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 2.91, (c) 3.05, (e) 2.72 (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. Setup B, θ ∆ = 90 • and ∆/lc = 0.004 (a-b), 0.008 (c-d), and 0.016 (e-f) at Ca = 0.042, τ = 17.41 (a-b), Ca = 0.047, τ = 22.76 (c-d), and Ca = 0.055, τ = 22.28 (e-f).

Figure 20 :

 20 Figure 20: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 3.42, (c) 3.79, (e) 3.67 (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. Setup C, θ ∆ = 90 • and ∆/lc = 0.004 (a-b), 0.008 (c-d), and 0.016 (e-f) at Ca = 0.082, τ = 15.79 (a-b), Ca = 0.09, τ = 22.5 (c-d), and Ca = 0.1, τ = 39.52 (e-f).

Figure 21 :

 21 Figure 21: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 2.94, (c) 3.22, (e) 3.37 (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. Setup B, θ ∆ = 110 • and ∆/lc = 0.004 (a-b), 0.008 (c-d), and 0.016 (e-f) at Ca = 0.058, τ = 21.67 (a-b), Ca = 0.064, τ = 32.88 (c-d), and Ca = 0.074, τ = 29.92 (e-f).

Figure 22 :

 22 Figure 22: Left panel: Comparison of the computed (black symbols) G(θe = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) with the best φ value of (a) 5.05, (c) 6.37, (e) 6.08 (black solid line). The fit is performed three to four grid points away from the contact line to minimize the inaccuracies amplified in our numerical method at the grid scale near the contact line. Right panel: Nondimensional curvature (black symbols) and the slope that the interface makes with the substrate (red symbols), as a function of the nondimensional vertical distance of the interface from the contact line, ζ/lc. Setup C, θ ∆ = 110 • and ∆/lc = 0.004 (a-b), 0.008 (c-d), and 0.016 (e-f) at Ca = 0.125, τ = 30.05 (a-b), Ca = 0.135, τ = 47.76 (c-d), and Ca = 0.155, τ = 35.43 (e-f).

Figs. 23

 23 Figs. 14 to 22 depict the curvature plots (black symbols) along with the angle (red symbols) θ a = θ(ζ/l c ) that the interface makes with the substrate as a function of ζ/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 θ ∆ = 110 • ; for example, for ∆ = 0.004, curvature is tending to zero around -0.75l c . However, the maximum curvature does not reach the asymptotic value of √ 2l -1 c when the asymptotic parameter = Ca/G(θ e ) is not small enough. As shown, for small angles, where the smallness of is guaranteed, the maximum curvature approaches √ 2l -1 c . 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[START_REF] Cox | The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow[END_REF] on Figs. 14 to 22. We plot the computed (black symbols) G(θ e = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆), using the computed interface slope, θ(r), at the distance r = x 2 i + y 2 i , 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(θ e = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) with the prediction of Eq. (23) using the best fitted value of φ (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(θ e = θ ∆ ) -G[θ(r)] versus Ca ln(r/∆) should be unity. In fact, there is a deviation from unity which is largest at small angles (15 • and 30 • ) 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 φ in the range where Cox's theory is expected to be asymptotically valid ∆ r r max . The expected upper limit r max = Ca 1/3 l c will be derived below. We consider the results for 0.004 ≤ ∆/l c ≤ 0.016 and show that the fitted values of φ are similar and thus approximately independent of ∆ when 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 r 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

Figure 23 :

 23 Figure 23: Interface profiles and the background adaptive mesh refinement for (a) Setup A, θ ∆ = 15 • , ∆/lc = 0.003, Ca = 0.0009, τ = 1.52, (b) Setup C, θ ∆ = 30 • , ∆/lc = 0.004, Ca = 0.004, τ = 5.37, (c) Setup A, θ ∆ = 40 • , ∆/lc = 0.007, Ca = 0.008, τ = 5.55, (d) Setup B, θ ∆ = 60 • , ∆/lc = 0.016, Ca = 0.022, τ = 11.12, (e) Setup B, θ ∆ = 90 • , ∆/lc = 0.004, Ca = 0.042, τ = 17.41, (f) Setup C, θ ∆ = 90 • , ∆/lc = 0.004, Ca = 0.082, τ = 15.79, (g)Setup B, θ ∆ = 110 • , ∆/lc = 0.004, Ca = 0.058, τ = 21.67, (h) Setup C, θ ∆ = 110 • , ∆/lc = 0.004, Ca = 0.125, τ = 30.05, and (i) Setup C, θ ∆ = 110 • , ∆/lc = 0.008, Ca = 0.135, τ = 47.76.

Figure 24 :

 24 Figure 24: The geometry of the wedge. The case of a large θe angle at the microscopic scale is illustrated. The various geometrical quantities s, r, η, ζ, θ and θe are discussed in the text.The figure emphasizes the important fact that our theory is also valid for θe > π/2 but that despite this fact we have r/s ζ/s 1 at large r.

Figure 26 :

 26 Figure 26: φ(θ ∆ ) as a function of ∆/lc for (a) Setup A, (c) Setup B, and (d) Setup C, and exemplifying the errors in the computation of φ for (b) Setup A, (d) Setup B, and (f) Setup C.

2 = 1 /3 2 - 1

 2121 κ III,max . (A.[START_REF] Mahady | A volume of fluid method for simulating fluid/fluid interfaces in contact with solid boundaries[END_REF] We now eliminate β between Eq. (A.4) and Eq. (A.7) and use the notation κ III,max = κ ∞ , 3

(B. 1 ) 3 (

 13 This equation can be turned into Airy's equation for a new variable z(s), that is z = sz, upon the substitutionsH(ζ) = z -2 (s),anddζ ds = -2 1/3 z -2 (s), whose solution is z(s) = αAi(s) + βBi(s), s) + βBi(s)] 2 , (B.3)where α and β are two arbitrary constants, while ζ(s) is an antiderivative of -2 1/3 z -2 (s). After integrating explicitly for ζ(s), Duffy and Wilson[START_REF] Duffy | A third-order differential equation arising in thin-film flows and relevant to tanner's law[END_REF] conclude that for β > 0 and any given α, there is a single branch that behaves as desired for the contact line problem, that it grows monotonically from H = 0 for ζ = 0 and for s → ∞ and has H → ∞ for s → s 1 and ζ → ∞ where s 1 is a root of αAi(s 1 ) + βBi(s 1 ) = 0. (B.4) That single branch cannot go to infinity inside s ∈ [s 1 , ∞), which implies that s 1 is the largest among the countably infinite set of roots of Eq. (B.4). Thus the solution can be characterized either 1) by the arbitrary pair α, β (β > 0) or 2) by β (β > 0) and s 1 chosen arbitrarily with α given by α = -βBi(s 1 )/Ai(s 1 ) provided there is no larger root of Eq. (B.4). In what follows, we use the second characterization. It can be verified that ζ = 2 1/3 π Ai(s) β[αAi(s) + βBi(s)] , (B.5) is an antiderivative of -2 1/3 z -2 (s) and that it satisfies ζ = 0 for s → ∞. Thus the solution given parametrically by Eqs. (B.3) and (B.5) starts at H = 0 for ζ = 0 and grows monotonically with H → ∞ for ζ → ∞. Expanding the solution for s 1, Duffy and Wilson [56] find that H(ζ) ∼ ζ[-3(ln ζ + c)] 1/3 , (B.6) where c = ln[π/(2 2/3 β 2 )]. Differentiating, and returning to the original variable η, leads to η (ζ) ∼ {9Ca ln[π/(2 2/3 β 2 eζ)]} 1/3 + O | ln(β 2 ζ)| -5/

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

theoretical value 0.0637. For setup C, in the case θ ∆ = 110 • , ∆/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 φ = C(q)/C(0)φ given the numerical estimate of Ca cr . The estimate is φ(θ ∆ , q) = πeAi 2 (s max )

If our theory is valid, the right hand side should not depend on ∆. 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 φ 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 φ 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 φ is obtained at small angles and we find φ(θ ∆ ) θ ∆ . 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 λ ∆. 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 θ ∆ is observed at small angles. Moreover, for all the Setups, large contact angle data do not follow the linear proportionality to the angle θ ∆ , with the viscosity ratio having a more prominent effect on this deviation at large angles. To shed more light on the behavior of φ as a function of the mesh size, as well as exemplifying the errors in the computation of φ, in Fig. 26, we plot φ(θ ∆ ) versus ∆/l c . The results show some small variations; this is however expected because the error in estimating φ also depends on ∆, while decreasing the mesh size does not reduce the transition capillary by much.

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 θ in depending on the contact line velocity so that

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

and we can proceed as before. We note however that there are two possibilities. Appendix A. Matching procedure in the small angle, free-surface case

In the small θ e case, lubrication theory can be used following ref. [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF]. 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 θ e Ca 1/3 . The range of validity of the theory is estimated in ref. [START_REF] Eggers | Existence of receding and advancing contact lines[END_REF] and for the advancing case to be Ca/θ 3 e < 0.05. This is equivalent for small angles to Ca/G(θ e ) ∼ δ < 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 [START_REF] Duffy | A third-order differential equation arising in thin-film flows and relevant to tanner's law[END_REF] η = 3Ca/η 2 , (

and we assume the scaling

so that the thickness of the film is Ca 2/3 much smaller than its horizontal extent Ca 1/3 . Eq. (A.1) can be solved using Airy functions. The analysis given in [START_REF] Duffy | A third-order differential equation arising in thin-film flows and relevant to tanner's law[END_REF] is reproduced in Appendix B. One finds that for ζ

where β is a parameter characterizing the Airy function solution of Eq. (A.1). Matching Eq. (A.3) with the Cox's solution Eq. ( 18) gives

The matching is performed for ζ small in region IV variables but large in region II variables, r m ζ Ca 1/3 l c . This is consistent with an upper bound r max = Ca 1/3 l c on the validity of the Cox's solution Eq. ( 13). At the upper end of region IV, for Ca 1/3 l c ζ l c , one can find the curvature as shown in Appendix B

Eq. (A.5) is obtained from the Airy function solution, which is parameterized by both β and s 1 . The determination of s 1 is more subtle [START_REF] Eggers | Hydrodynamic theory of forced dewetting[END_REF][START_REF] Chan | Theory of the forced wetting transition[END_REF]. It is seen that expression (A.5) predicts a range of possible curvatures depending on the distances, the curvature of the parabola (Eq. (B.8)) is close to the curvature at its apex. Reverting to the original variable η, we obtain the curvature κ ∞ η (ζ) in the asymptotic range Ca -1/3 l c ζ l c as κ ∞ ∼ (3Ca) 1/3 2 1/6 β πAi(s 1 )

2 .

(B.9)