Cox–Voinov theory with slip

Abstract Most of our understanding of moving contact lines relies on the limit of small capillary number ${Ca}$. This means the contact line speed is small compared to the capillary speed $\gamma /\eta$, where $\gamma$ is the surface tension and $\eta$ the viscosity, so that the interface is only weakly curved. The majority of recent analytical work has assumed in addition that the angle between the free surface and the substrate is also small, so that lubrication theory can be used. Here, we calculate the shape of the interface near a slip surface for arbitrary angles, and for two phases of arbitrary viscosities, thereby removing a key restriction in being able to apply small capillary number theory. Comparing with full numerical simulations of the viscous flow equation, we show that the resulting theory provides an accurate description up to $Ca \approx 0.1$ in the dip coating geometry, and a major improvement over theories proposed previously.


Introduction
The theory of the moving contact line at small capillary numbers was founded by Voinov (1976) and generalized to arbitrary viscosity ratios M by Cox (1986). The problem is that if the no-slip boundary condition were to apply down to arbitrarily small scales (Huh & Scriven 1971;Bonn et al. 2009;, a contact line would not be able to move. Therefore, one needs to invoke a small length scale on which the conventional equations for fluid motion are relaxed. The simplest, and often physically realistic, such choice is the introduction of a Navier slip length λ (Lauga, Brenner & Stone 2008), over which a fluid may slip past a solid interface. Thus we will always assume that, as the wall is approached, the tangential velocity relative to the wall is λ times the shear rate. The slip length is generally of the order of a few nanometres, but increases somewhat for hydrophobic surfaces (Barrat & Bocquet 1999;Cottin-Bizonne et al. 2005). A slip length is used very widely in contact line problems (Bonn et al. 2009;Vandre, Carvalho & Kumar 2014;Sprittles 2015) in order to regularize the local flow and to thus allow contact line motion, but is usually not important elsewhere in the flow. The inner length scale λ is contrasted with an outer 'macroscopic' length scale R, for example the radius of a spreading drop or the capillary length scale in the problem. Cox (1986) clarified the structure of low capillary number problems in terms of the ratio = 0 λ/R between the two length scales; 0 is a numerical factor to be determined. From a general analysis, Cox (1986) obtained (see also Sibley, Nold & Kalliadasis 2015) g(θ eq ) = g(θ app ) + Ca ln , (1.1) where for a single phase (a liquid volume, neglecting the effect of a surrounding gas). Here, θ eq is the (equilibrium) contact angle on the microscopic scale, and θ app is the apparent contact angle. In the case of a small spreading drop, this is the angle a spherical cap, fitted to the drop, makes with the solid. In Cox's analysis, the dimensionless capillary number is assumed small, where U is the speed of the contact line. Here, η is the fluid viscosity, and γ is the surface tension. For steady flow, it is often more convenient to instead think of the contact line being stationary and the wall to be moving at speed U, because we can then think of the interface as time independent.
In the present paper, we take the approach proposed by Snoeijer (2006) and Chan et al. (2013), who describe the interface shape as an evolution equation for the local slope θ . For simplicity, we neglect external forcing (such as gravity), which only becomes important on a macroscopic scale. The requirement that the interface (described by its thickness h(s)), is assumed steady, while the substrate is moving at speed U to the right, then leads to the 'generalized lubrication' (GL) equation. One obtains (Snoeijer 2006) where s is the arclength along the interface, and M = η g /η l the ratio of the viscosity of the outer phase (for example a gas), and the viscosity of the liquid (cf. figure 1). In the case of a single (liquid) phase only, .
(1.5)  In contrast to the approach of asymptotic matching between micro and macroscales adopted by Cox (1986) or Hocking & Rivers (1982), the GL equation (1.4), for small Ca, gives a full representation of the interfacial profile, given that boundary conditions at the contact line and the far field are properly imposed, for arbitrary contact angles and viscosity ratios. On the large scale, boundary conditions may be the requirement for the interface to assume a drop shape, or to asymptote to a flat interface, as we will require below in our numerical tests of the equation.
The unknown length scale appearing in Cox's (1986) theory is now no longer required, but should in principle come out of the solution of (1.4). However, the slip effect is not properly implemented in (1.4) and as a result, the equation does not have solutions corresponding to a moving contact line. To account for this, Chan et al. (2013), motivated by the form of the lubrication equation, introduced the modified GL equation to yield where c is a constant to be chosen such that the solution matches properly with the slip region. In Chan et al. (2013), in the absence of any theory for how c should be chosen, this was done assuming that c = 3, which is true only for small angles and for a single phase (Hocking 1983). Here, we calculate the dependence of c on the microscopic contact angle θ eq and on the viscosity ratio M. Once c is known, one can integrate (1.7) from the contact line to whichever macroscopic configuration is required for the problem. Let us start by sketching the structure of the analysis to compute c. For small Ca, it is sufficient to consider a linear perturbation of the free surface shape around a wedge with microscopic (equilibrium, for simplicity and without loss of generality) angle θ eq . Therefore we set θ = θ eq − ϕ(s), and, to linear order in ϕ, all calculations can be performed assuming a wedge geometry h = s sin θ eq for the flow (cf. figure 1). Linearizing (1.7) in ϕ and integrating twice, we have We are interested in solutions which only grow logarithmically, corresponding to vanishing curvature at infinity, and so C = 0. From the boundary condition ϕ(0) = 0 it follows that C 1 = 3Ca F(θ eq , M)(1 − ln cλ)/ sin 2 θ eq , and so to leading order as s/λ → ∞ we have ϕ(s) = 3Ca F(θ eq , M) sin 2 θ eq ln se sin θ eq cλ . (1.8) Here, we calculate c for arbitrary angles θ eq and arbitrary viscosity ratios M, based on earlier work by Hocking (1977), who calculates the stress on a slip wall, assuming a straight interface h = s sin θ eq . Using the fact that in the absence of inertia (Stokes dynamics), the total force on any fluid volume vanishes, we can convert the wall force to a force on the interface, which leads to bending of the interface, allowing us to determine c. In addition, one may also allow for an arbitrary ratio of slip lengths λ 1 /λ 2 at the two fluid-solid interfaces, but no explicit results are available for this case. So in the interest of simplicity, we will always assume λ 1 = λ 2 ≡ λ.

A single phase
We start the analysis by considering the case of a single phase (the liquid, as shown in figure 1), which enables us to show some of the calculations in greater detail. The idea is to consider a force balance over the volume shown in figure 1, within which the Stokes equation ∇ · σ = 0 is satisfied; σ is the stress tensor. This is motivated by Hocking's (1977) calculation of the force exerted by a moving contact line on a solid boundary. Using Gauss' theorem, and only considering the force in the x-direction, we obtain S n · σ · e x = 0, (2.1) where n is the outward normal, and S = S i + S c + S w . The volume is the slice of radius s inside the fluid, where S i is the fluid-gas interface, S c a circular arc of radius s inside the fluid and S w the wall, see figure 1. The force on V coming from the interface is Here, we have used the stress boundary condition (Landau & Lifshitz 1984) n · σ = −γ nκ, where κ = dϕ/ds is the interface curvature and n · e x − sin θ eq . The integral over S w , representing the total force w on the wall between the origin and x = s, has been calculated by Hocking (1977) for Stokes flow in a corner with slip at the wall, and a free-slip condition at the interface. The contribution v of the force on S c can be inferred from the far-field limit of the flow, calculated by Huh & Scriven (1971). Thus (2.1) gives which we compare to (1.8) to find c.
To find v and w, we need to consider the flow in the wedge-shaped fluid domain (cf. figure 1) of opening angle θ = θ eq . In Hocking (1977), the flow is solved subject to the slip boundary condition (with u the horizontal component of the velocity) on the wall. It is convenient to count the angle φ from the interface; the velocity field is given in terms of the streamfunction The boundary conditions on both boundaries are zero normal velocity, and vanishing shear stress σ rφ = 0 on the interface. Thus for φ = θ eq (the wall) we have In the limit r λ we can neglect the effects of slip, and we should fall back on the similarity solution ψ = rf (θ eq ) found by Huh & Scriven (1971). With the ansatz the boundary conditions are The coefficients are (D = θ eq − sin θ eq cos θ eq ), a 1 = − Uθ eq cos θ eq D , a 2 = 0, a 3 = 0, a 4 = − U sin θ eq D , and so the stresses become Using (2.5a,b) and n · σ · e x = σ rr cos(θ eq − φ) + σ rφ sin(θ eq − φ), where F(θ eq ) is the angle dependence from the GL model given previously in (1.5). As an aside, F(θ eq ) can easily be calculated from σ φφ , evaluated at the interface φ = 0, and using the stress boundary condition γ κ = −σ φφ | φ=0 . On the other hand, using h = s sin θ eq and integrating (1.4) once, we have comparing with σ φφ | φ=0 yields F(θ eq ). Finally, Hocking (1977) has calculated the force on the wall where h 1 (θ eq ) is the solution of an integral equation, which in general has to be solved numerically. Using (2.3), it follows that, in the general case, which is precisely of the expected form (1.8). Comparing the two, we find which is our main result in the case of a single phase. The function c(θ eq ) as obtained from (2.7) is shown in figure 2. To compute h 1 (θ eq ), we have solved the integral equation (3.6) of Hocking (1977) numerically using a relatively simple iterative procedure. The red cross is the exact value h 1 (π/2) = 4(γ E − ln 2)/π, giving c(π/2) = exp(ln 2 − γ E ) ≈ 1.12. (2.8) For small values of θ eq , Hocking's (1977) asymptotic argument leads to h 1 = −k lnk, wherek = −3F(θ eq )/ sin θ eq . However, in order to obtain a consistent expansion,k had to be expanded to the next order beyond (3.18) of (Hocking 1977), which results in c = 3 − (9/10)θ 2 eq . (2.9) For θ eq = 0, one of course recovers c = 3, as known from lubrication theory. On the other hand for θ eq → π, using Hocking's (1977) leading-order expression for h 1 , we obtain c = (π − θ eq ) exp −π/(π − θ eq ) . (2.10) The limits (2.9) and (2.10) are shown in figure 2 as the red dashed lines, and are seen to agree very well with our numerical solution. Interestingly, c becomes very small for angles close to π. For example, c(0.9π) = 1.04 · 10 −4 , showing that in general c cannot be inferred from a comparison with lubrication theory.

Two phases
We now generalize to the case of two phases (arbitrary viscosity ratio M), for which the constant c has never been calculated, not even for small angles. We explain the general structure, but explicit results are available for θ eq = π/2 only. θ eq /π c(θ eq ,0) FIGURE 2. The c-factor as a function of θ eq /π for M = 0, calculated numerically by evaluating the integral expressions for h 1 as given by Hocking (1977). The (red) dashed lines are the asymptotic expansions (2.9) and (2.10) for small θ eq and as θ eq → π, respectively. The (red) cross marks the analytical value c(π/2) = exp(ln 2 − γ E ) (see (2.8)).
In contrast to the previous subsection, we now have the integral (2.1) over the surface of two volumes V 1 and V 2 ; V 1 is the volume over the 'liquid' phase as before (labelled 1), V 2 is the corresponding slice of the same radius s over phase 2 (the 'gas'). Then (2.1), written for each phase, becomes where v and w have the same meanings as before, but for each phase separately. Using the stress condition n · (σ 1 − σ 2 ) = −γ nκ, this yields The solution for a no-slip flow, for general M, has been given by Huh & Scriven (1971). At the interface between the two fluids we now have continuity of the tangential velocity (the normal velocity vanishes), and of shear stress. If the streamfunctions in either phase are ψ 1/2 = rf 1/2 (θ eq ), the boundary conditions are The results are a little too complicated to write out here, but as in (2.6), we can calculate the general form of F(θ, M) in (1.4), by computing the normal stress on the interface from either phase, giving Comparison with (1.4) then shows that note the sign error in Chan et al. (2013). Furthermore, it follows from the Huh-Scriven solution, after a lengthy but elementary calculation, that The shear force on the wall now has contributions from both phases, and it follows from Hocking (1977) that Taken together with (2.12), this yields the interface deformation and comparing with (1.8) (2.14) now involving two numerical constants h 1 and h 2 . This completes our calculation, but it remains to calculate the constants h 1 and h 2 , which are functions of θ eq and M. According to Hocking (1977) this can be done explicitly for a right angle θ eq = π/2, for which where h a = (4/π)(γ E − ln 2) and h b = −1.539. The geometrical factor is

Discussion
To illustrate and validate our results, we consider a plate being pushed steadily into a container of fluid 2, with an outer fluid 1, as shown in figure 4(a). This is the classical dip-coating configuration, used industrially to coat a solid with liquid, and studied in many experiments (Blake & Ruschak 1979 one has to add a term − cos(θ ) to the modified GL equation (1.7) (Chan et al. 2013). One of the main aims is to calculate the apparent contact angle as well as the critical capillary number at which fluid 1 is entrained into fluid 2. As shown in Kamal et al. (2018), for small values of M (for example when fluid 1 is air), the critical capillary number is so high that it cannot be calculated within an expansion for small Ca. However, here we focus on the regime of smaller Ca up to 0.1, and compute the apparent contact angle, to assess the improvement of our predictions over Chan et al. (2013), where the results it was assumed that c = 3. We compare full numerical simulations of the dip-coating problem, using the finite elements method (FEM) described in detail in Kamal et al. (2018), and based on a framework benchmarked in Sprittles (2015), with results of the GL equation for Ca up to 0.1. We show an example with M = 0 and one with M = 1, with a contact angle of π/2 (figure 4(b,c), respectively). On the left, the solid line is the interface profile found from FEM, the red dashed line comes from integrating (1.7) using the new expression of c in (2.7), as obtained by the present theory; very good agreement is found.
To demonstrate the importance of using the correct value of c in order to obtain this agreement, we also plot the profile obtained using the value c = 3, proposed previously (Snoeijer 2006;Chan et al. 2013). However, this value is only appropriate in the limit of small contact angles and no outer fluid (M = 0). As a result, the dot-dashed blue, which is the profile thus obtained from the 'old' GL equation, differs significantly from the 'exact' result, obtained from direct numerical simulation, including poorly predicting the (macroscopic) height far away from the contact line.
To focus on the solutions very close to the contact line, we also plot the relative deviation between the FEM simulation and solutions obtained from integrating (1.7) on the right of figure 4. For the solid line we use the correct value of c, the dashed line represents the value c = 3. For the smaller capillary number, using the c-value as calculated in the present paper, the deviation is negligible, and even for a capillary number Ca = 0.1, the relative error remains small; when the capillary number becomes of order unity the GL equation fails (Kamal et al. 2018), since it is an expansion for small Ca only. On the other hand for c = 3 the relative error is significant in both cases as one approaches the contact line.  Having benchmarked the new theory against FEM calculations, in order to assess the improvement of the current theory over the previous GL theory using c = 3, we compare the apparent contact angle θ app in both cases, for a wide range of parameters M and θ eq . For the dip-coating geometry shown in figure 4, θ app can be calculated by comparing the depth Δ of the contact line relative to the level of the bath to the depression of a static meniscus (Landau & Lifshitz 1984) Δ = 2(1 − sin θ app ) (once more in units of l γ ). We use the GL equation to evaluate the error in predictions of θ app between simulations where the correct value of c is implemented, to those where the lubrication value (c = 3) is used (figure 5). For finite viscosity ratios M, above a critical capillary number Ca cr the stationary meniscus bifurcates to a time-dependent solution, in which the outer fluid is entrained Kamal et al. 2018). In figure 5(a), we vary the viscosity ratio M as the contact angle is held constant at θ eq = π/2. It shows the considerable error in θ app , which increases rapidly as Ca cr is approached. In figure 5(b), we hold M = 0 constant, while varying θ eq between π/2 and π. As expected, the error in θ app increases substantially as θ eq increases. In fact, as seen in (2.10), the true value of c decreases exponentially as θ eq = π is approached.
In conclusion, we have analysed the Cox-Voinov theory with slip, for arbitrary contact angles and arbitrary viscosity ratios. Closed form expressions are provided using the generalized lubrication equation, which is shown to be a powerful tool to calculate interface flows involving moving contact lines in situations in which the interface slope is not necessarily small. The present paper provides a consistent version of the GL equation for arbitrary contact angles.