1. Introduction
Understanding the way in which a contact line (the line along which an interface separating two immiscible fluids intersects a solid boundary wall) moves under an imposed flow is a problem of fundamental importance in wetting dynamics. It underpins a wide range of phenomena in nature, for example in the functional adaptation of many biological systems (Parker & Lawrence Reference Parker and Lawrence2001; Hu, Chan & Bush Reference Hu, Chan and Bush2003; Zheng, Gao & Jiang Reference Zheng, Gao and Jiang2007), as well as in technological applications: from oil recovery (Morrow Reference Morrow1990) to inkjet printing and microfluidics (Tabeling Reference Tabeling2010). From these examples, it has become clear that experimental progress needs to be accompanied by the development of physical models and computational methods that can accurately capture wetting dynamics, even in complex geometries.
Any purely macroscopic treatment – in which the two fluids are considered to be separated by a strictly sharp interface, and subject to a strictly no-slip condition on the fluid velocity at the solid wall – has long been recognised to predict a non-integrable stress singularity for any non-zero velocity of the contact line (Huh & Scriven Reference Huh and Scriven1971; de Gennes Reference de Gennes1985). This so-called moving contact line singularity is clearly at odds with the commonplace experience that the contact line does, in fact, move.
To regularise this singularity, the macroscopic picture just described must be supplemented by additional physics that intervenes on shorter, microscopic length scales. One approach is to relax the no-slip boundary condition on the fluid velocity, by inserting a slip boundary condition in a region of microscopic size
${\it\xi}$
in the vicinity of the contact line. This can be done by hand (Zhou & Sheng Reference Zhou and Sheng1990; Spelt Reference Spelt2005), for example by imposing a slip velocity
$v_{\text{s}}=U\exp (-|x|/{\it\xi})$
, where
$x$
is the distance along the wall away from the contact line, or by choosing
$v_{\text{s}}={\it\xi}{\it\sigma}/{\it\eta}$
, where
${\it\sigma}$
and
${\it\eta}$
are the wall stress and fluid viscosity, respectively. In any computational study, slip can also arise indirectly as a numerical artefact (Renardy, Renardy & Li Reference Renardy, Renardy and Li2001).
Another approach is to keep the no-slip boundary condition and instead to recognise that any interface between two fluids can never be perfectly sharp, but in practice must have some non-zero microscopic diffuse width
$l$
. Accordingly, the use of diffusive interface models (Anderson, McFadden & Wheeler Reference Anderson, McFadden and Wheeler1998; Jacqmin Reference Jacqmin2000) has gained increasing popularity in recent years. A diffuse interface model of liquid–gas coexistence then accommodates contact line motion by a condensation–evaporation mechanism which, in some small region within close proximity of the contact line, slowly transfers matter from one side of the interface to the other (Seppecher Reference Seppecher1996; Briant, Wagner & Yeomans Reference Briant, Wagner and Yeomans2004; Diotallevi et al.
Reference Diotallevi, Biferale, Chibbaro, Pontrelli, Toschi and Succi2009; Pooley, Kusumaatmaja & Yeomans Reference Pooley, Kusumaatmaja and Yeomans2009). Similarly, a diffuse interface model of coexisting immiscible binary fluids accommodates contact line motion via a slow intermolecular diffusion of the two fluids across the interface between them, again acting in a small region in the vicinity of the contact line (Jacqmin Reference Jacqmin2000; Briant & Yeomans Reference Briant and Yeomans2004; Pooley et al.
Reference Pooley, Kusumaatmaja and Yeomans2009).
Diffuse interface models (Anderson et al.
Reference Anderson, McFadden and Wheeler1998; Jacqmin Reference Jacqmin2000) are also convenient in computational practice. An order parameter (‘phase field’) is introduced to distinguish one fluid from another, with the length scale that characterises the variation of this field across the interface between the fluids prescribing the interfacial width,
$l$
. An advection–diffusion equation (the Cahn–Hilliard equation supplemented by an advective term) then determines the evolution of this order parameter, and hence of the interface position. This is solved in tandem with the Navier–Stokes equation for the fluid velocity field. A key advantage of this approach is that all computational grid points are treated on an equal footing, without any need to explicitly track the position of the interface. In this way, even complex flow geometries can be simulated conveniently.
It is important to emphasize, as first elucidated by Cox (Reference Cox1986), that although the physical origins of contact line motion may differ according to the detailed microscopic physics invoked (wall slip, intermolecular diffusion across the interface, etc.), the macroscopic hydrodynamics far from the contact line nonetheless converges to a universal solution that is informed by the microscopics only in the sense of being rescaled by a single parameter, the ‘slip length’
${\it\xi}$
, that emerges from this underlying microscopic physics. An important problem within any model of the microscopics is therefore to determine this emergent slip length
${\it\xi}$
: not only because it is a key variable that determines the macroscopic wetting and fluid dynamics, but also because it sets a scale to which other length scales (surface heterogeneities, droplet size, etc.) must be compared.
With that motivation in mind, the primary aim of this work is to determine the scaling of the slip length within a diffuse interface model of immiscible binary fluids. We assume the fluids to have matched viscosity
${\it\eta}$
, and denote by
$M$
the mobility parameter characterising the rate of intermolecular diffusion in the Cahn–Hilliard equation. In the existing literature, (at least) three apparently contradictory scalings have been proposed. Several authors (e.g. Jacqmin Reference Jacqmin2000; Yue, Zhou & Feng Reference Yue, Zhou and Feng2010) suggest that
${\it\xi}\sim l_{D}$
where
$l_{D}=(M{\it\eta})^{1/2}$
is the characteristic length scale below which intermolecular diffusion dominates advection and above which the opposite is true. For convenience we call this the diffusion length in what follows, and we emphasise that it is determined only by the macroscopic quantities
$M$
and
${\it\eta}$
. In contrast, Briant & Yeomans (Reference Briant and Yeomans2004) suggest a different scaling,
${\it\xi}\sim \sqrt{l_{D}l}$
, which is the geometric mean of the macroscopic diffusion length and the microscopic interface width. And while these two proposed scalings do not have any dependence on the imposed flow velocity, Chen, Jasnow & Viñals (Reference Chen, Jasnow and Viñals2000) in contrast suggested that the slip length does depend on the velocity as
${\it\xi}\propto V_{0}^{-1/2}$
.
In view of these apparently contradictory predictions, an important contribution of this work is to provide a coherent picture that fully reconciles two of these existing predictions, by showing each to apply in a different regime of parameter space, and that carefully appraises the third. This is achieved by extensive numerical studies performed across unprecedentedly wide ranges of the two relevant control parameters. Specifically, we explore four decades of the dimensionless ratio
$M{\it\eta}/l^{2}=l_{D}^{2}/l^{2}$
of the (squared) macroscopic diffusion length to the (squared) microscopic interface width; and three decades of the capillary number
$Ca$
(defined below), which adimensionalises the imposed flow velocity
$V_{0}$
in terms of an intrinsic interfacial velocity scale. Moreover, we use three different, independently coded, numerical methods and show their results to be in excellent quantitative agreement across all decades.
For small values of the capillary number
$Ca$
, corresponding to low imposed flow velocities, we demonstrate that the slip length converges to a well-defined value that is independent of the velocity
$V_{0}$
. Any dependence on the velocity is only observed for
$Ca>0.01$
, and we comment on the extent to which our results are consistent with the prediction
${\it\xi}\propto V_{0}^{-1/2}$
of Chen et al. (Reference Chen, Jasnow and Viñals2000) in this fast flowing regime.
The remainder of the paper focuses on the slow-flow regime,
$Ca<0.01$
, in which the slip length
${\it\xi}$
is independent of the flow rate
$V_{0}$
. Here we show that, in the limit of large
$l_{D}/l$
, the slip length scales as
${\it\xi}\propto l_{D}$
, informed only by the macroscopic length scale
$l_{D}$
. This agrees with the original prediction of Jacqmin (Reference Jacqmin2000), Yue et al. (Reference Yue, Zhou and Feng2010) and corresponds to the sharp interface limit of the diffuse interface model. In this limit the microscopic length
$l$
essentially drops out of the problem, apart from providing a ‘singular perturbation’ to regularise the contact line singularity. In contrast, for small
$l_{D}/l$
we find that the slip length scales as
${\it\xi}\propto (l_{D}l)^{1/2}$
, as proposed by Briant & Yeomans (Reference Briant and Yeomans2004). This corresponds to the diffuse interface limit in which the emergent slip length does depend on the underlying microscopic length
$l$
. In this way we reconcile these two previously apparently contradictory scalings, by showing each to apply in a different limiting regime of
$l_{D}/l$
. The cross-over between the two is furthermore consistent with the onset of the sharp interface regime discussed by Yue et al. (Reference Yue, Zhou and Feng2010).
The fact that the slip length
${\it\xi}$
is in general different from the interface width
$l$
, and indeed greatly exceeds it for large
$l_{D}/l$
, has a clear practical consequence for any simulation: the resolution of the full wetting dynamics appears only on a length scale corresponding to the larger of the interface width and the slip length. This, in turn, limits the region of parameter space that can be considered reliable in any numerical study.
We provide additional evidence for the distinction between the diffuse interface limit
$l_{D}/l\ll 1$
, in which the microscopic physics informs
${\it\xi}$
, and the sharp interface limit
$l_{D}/l\gg 1$
, in which it does not, by performing further simulations in which the we modify the interfacial (microscopic) contribution to the underlying free-energy functional of the binary fluid. We show that the slip length depends on this modification only in the diffuse interface regime. In the sharp interface regime of large
$l_{D}/l$
, the slip length continues to depend only on the macroscopic dynamical variables,
${\it\eta}$
and
$M$
, free of this modification to the microscopic interfacial term.
Finally, we test the validity of Cox’s formula (Cox Reference Cox1986) for the dependence of the dynamic contact angle on the Capillary number in the different slip length regimes. Our numerical results are in good agreement with Cox’s analytical result if we allow the slip length to be rescaled by a dimensionless parameter. Moreover this parameter appears, suggestively, to converge to unity in the sharp interface limit, but is smaller in the diffuse interface limit.
This paper is organized as follows. In the next section we describe the models, methods and set-ups we employ to compute the slip lengths and dynamic contact angles. We present our results in § 4. Finally, we summarise our key findings and discuss avenues for further work in § 5.
2. Model, geometry and boundary conditions
In this section we specify the model that we shall use throughout the paper, starting with the thermodynamics in §§ 2.1 and 2.2 then the dynamical equations of motion in § 2.3. We then specify the flow geometry and boundary conditions in § 2.4.
2.1. Thermodynamics
We consider a binary mixture of two mutually phobic fluids, A and B, and denote the volume fraction of fluid A by the continuum phase field
${\it\phi}(\boldsymbol{r},t)$
. The volume fraction of fluid B is then simply
$1-{\it\phi}$
by mass conservation and need not be considered separately. We consider a Landau free energy (Bray Reference Bray1994)
$$\begin{eqnarray}\displaystyle {\rm\Psi}=\int _{V}{\it\psi}_{b}\,\text{d}V=\int _{V}\left[\frac{G}{4}({\it\phi}^{2}-1)^{2}+\frac{Gl^{2}}{2}(\boldsymbol{{\rm\nabla}}{\it\phi})^{2}\right]\text{d}V, & & \displaystyle\end{eqnarray}$$
which allows a coexistence of two bulk phases: an A-rich phase with
${\it\phi}_{A}=1$
and a B-rich phase with
${\it\phi}_{B}=-1$
. The bulk constant
$G$
has dimensions of energy per unit volume (or equivalently of force per unit area, and so modulus). If the two fluids have different affinities to the solid surface, a surface (wetting) contribution to the free energy can be added (Cahn Reference Cahn1977; Pooley, Kusumaatmaja & Yeomans Reference Pooley, Kusumaatmaja and Yeomans2008; Pooley et al.
Reference Pooley, Kusumaatmaja and Yeomans2009) to the right-hand side of (2.1). Throughout this work we assume neutral wetting conditions, for which no such contribution is needed.
The chemical potential
${\it\mu}$
follows as a functional derivative of the free energy density
${\it\psi}_{b}$
with respect to
${\it\phi}$
, giving
$$\begin{eqnarray}{\it\mu}=-G{\it\phi}+G{\it\phi}^{3}-Gl^{2}\,{\rm\nabla}^{2}{\it\phi}.\end{eqnarray}$$
In equilibrium, the chemical potential
${\it\mu}=0$
. Assuming a flat interface of infinite extent in the
$y{-}z$
plane, with the interfacial normal in the
$x$
direction and the interface located at
$x=0$
, we then obtain an interfacial solution of the form
$$\begin{eqnarray}{\it\phi}=\tanh \left(\frac{x}{\sqrt{2}l}\right),\end{eqnarray}$$
with a homogeneous B-rich phase for
$x\ll -l$
and homogeneous A-rich phase for
$x\gg l$
. The interfacial constant
$l$
specifies the characteristic length scale over which
${\it\phi}$
varies in between these two phases, and so corresponds to the interfacial width. The surface tension associated with this interface is given by
$$\begin{eqnarray}{\it\sigma}=\frac{2\sqrt{2}Gl}{3}.\end{eqnarray}$$
2.2. Generalisation of the Landau theory: introduction of a curvature term
So far, we have specified the Landau free energy in the form most commonly used in the literature. It will also be instructive in what follows to consider the extent to which our results for the slip length do or do not depend on the microscopic details of the model used. Accordingly, we now generalise the free energy slightly to give
$$\begin{eqnarray}{\rm\Psi}^{\prime }=\int _{V}{\it\psi}_{b}^{\prime }\,\text{d}V=\int _{V}\left[\frac{G}{4}({\it\phi}^{2}-1)^{2}+\frac{Gl^{2}}{2}(\boldsymbol{{\rm\nabla}}{\it\phi})^{2}+{\it\alpha}\frac{Gl^{4}}{4}({\rm\nabla}^{2}{\it\phi})^{2}\right]\text{d}V.\end{eqnarray}$$
Compared to the original free energy defined above, this has an additional interfacial curvature term of amplitude set by
${\it\alpha}$
. The mapping between this form of Landau free energy and the (Helfrich) continuum elastic energy can be found in e.g. Gompper & Zschocke (Reference Gompper and Zschocke1991). The associated chemical potential is
$$\begin{eqnarray}{\it\mu}=G\left(-{\it\phi}+{\it\phi}^{3}-l^{2}\,{\rm\nabla}^{2}{\it\phi}-{\it\alpha}l^{4}\,{\rm\nabla}^{4}{\it\phi}\right).\end{eqnarray}$$
It is important to note that the bulk terms are unaltered, with the modification affecting only interfacial gradient terms involving powers of
$l\,\boldsymbol{{\rm\nabla}}$
. A careful comparison of the original model for which
${\it\alpha}=0$
, with this generalised model for which
${\it\alpha}>0$
, will allow us to demonstrate that the slip length is independent of the microscopic details, as specified by
${\it\alpha}$
, in the sharp interface regime
$l_{D}/l\gg 1$
. In contrast in the diffuse interface regime
$l_{D}/l\ll 1$
we find that the slip length does depend on the microscopics, via
${\it\alpha}$
.
2.3. Equations of motion
The dynamics of the order parameter field is specified by the Cahn–Hilliard equation (see e.g. Bray Reference Bray1994) generalised to include an advective term:
$$\begin{eqnarray}\left(\partial _{t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right){\it\phi}=M{\rm\nabla}^{2}{\it\mu}.\end{eqnarray}$$
Here
$M$
is the molecular mobility, which we assume constant. The fluid velocity and pressure fields,
$\boldsymbol{v}(\boldsymbol{r},t)$
and
$p(\boldsymbol{r},t)$
, obey the continuity and Navier–Stokes equations
$$\begin{eqnarray}\displaystyle & \partial _{t}{\it\rho}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\rho}=-{\it\rho}\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & {\it\rho}\left(\partial _{t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\right)\boldsymbol{v}={\it\eta}{\rm\nabla}^{2}\boldsymbol{v}-\boldsymbol{{\rm\nabla}}p-{\it\phi}\boldsymbol{{\rm\nabla}}{\it\mu}. & \displaystyle\end{eqnarray}$$
We denote by
${\it\rho}$
and
${\it\eta}$
the fluid density and viscosity respectively, assuming throughout that the two fluids are perfectly matched in both density and viscosity. In the Navier–Stokes equation, gradients in the chemical potential contribute an additional forcing term to the fluid motion,
$-{\it\phi}\boldsymbol{{\rm\nabla}}{\it\mu}$
(Jacqmin Reference Jacqmin1999).
Diffuse interface models are widely used in the computational fluid dynamics community and they have been implemented using various approaches. To ensure our results are free from numerical artefacts, here we have used three different numerical methods: (MI) the lattice Boltzmann method, (MII) the spectral method, and (MIII) the immersed boundary method to solve the equations of motion. They are detailed in § 3 below. For MII and MIII, we have assumed incompressible flow and taken the inertialess limit of zero-Reynolds-number Stokes flow. This corresponds to setting to zero the terms on the left-hand sides of (2.8) (incompressibility) and (2.9) (zero inertia). The lattice Boltzmann method (MI) intrinsically requires a small but non-zero inertia and compressibility, though our numerical results confirm the effects of this difference to be negligible for the problem considered here.
2.4. Geometry, initialisation and boundary conditions
Schematic illustration of typical steady states in our channel geometry. Marked are the key length scales: slip length
${\it\xi}$
, interface displacement
$d$
, channel dimensions
$L_{x},L_{y}$
and diffuse interface width
$l$
. In (a) a shear flow is imposed by moving the plates with equal and opposite velocities
$\pm V_{0}$
. In (b) a Poiseuille flow is imposed by applying a pressure drop along the length of the channel, sketched here in the frame in which the fluid bridge is stationary. The definitions for the dynamic (
${\it\theta}_{d}$
, macroscopic) and equilibrium (
${\it\theta}_{e}$
, microscopic) contact angles are shown.
We consider flow between infinite flat parallel plates a distance
$L_{y}$
apart, with plate normals in the
$y$
direction. In the flow direction
$x$
the cell is taken to have length
$L_{x}=2L_{y}$
, with periodic boundary conditions on both
${\it\phi}$
and
$\boldsymbol{v}$
. All quantities are assumed invariant in the
$z$
direction. The phase field is initialised in an equilibrium state with a bridge of A-rich phase, in which
${\it\phi}=+1$
, separated by two vertical diffuse interfaces of width
$\ell$
at
$x=L_{x}/4$
and
$3L_{x}/4$
from B-rich phases on either side, where
${\it\phi}=-1$
. The equilibrium contact angle is taken to be
${\it\theta}_{e}=90^{\circ }$
, corresponding to neutral wetting conditions. Throughout we define the position of the interface itself by the locus
${\it\phi}=0$
.
The fluid is taken to be initially at rest with
$\boldsymbol{v}=0$
everywhere. Starting from this initial condition, we then implement one of two common flow protocols: boundary-driven planar Couette flow and pressure-driven planar Poiseuille flow. In the first of these, a constant shear rate
$\dot{{\it\gamma}}$
is applied by moving the top and bottom plates at velocities
$V_{0}=\pm \dot{{\it\gamma}}L_{y}/2$
. This deforms the phase field and for small
$Ca\equiv {\it\eta}V_{0}/{\it\sigma}\lesssim 0.1$
a steady state is reached in which the interface has displaced a distance
$d\propto \pm Ca$
at the top and bottom walls, as shown in figure 1(a). The deformed bridge is then stationary, with the contact lines moving at a velocity
$\mp V_{0}$
relative to the top and bottom walls.
In the Poiseuille protocol, the flow is driven by an imposed pressure drop
${\rm\Delta}P$
along the length of the channel. A steady state then develops in which the bridge migrates along the channel at a constant speed, with the contact lines moving along each wall with a measured velocity
$V_{0}$
. In the reference frame of the contact line, therefore, both walls move with a velocity
$-V_{0}$
. By fitting the interface in the central region of the channel (between
$y=1/4\rightarrow 3/4$
) to a circle, see figure 1(b), we define the dynamic contact angle
${\it\theta}_{d}$
as the angle of intersection this circle makes with the wall. The dynamic contact angle is in general different from the equilibrium contact angle.
At the plates we assume boundary conditions of no slip and no permeation for the fluid velocity. For the phase field, we implement
$\partial _{y}{\it\mu}|_{y=0,L_{y}}=0$
in method MIII. In method MII it is instead more convenient to set
$\partial _{y}{\it\phi}|_{y=0,L_{y}}=\partial _{y}^{3}{\it\phi}|_{y=0,L_{y}}=\partial _{y}^{5}{\it\phi}|_{y=0,L_{y}}=0$
(where the last equality need only be imposed in the case
${\it\alpha}\neq 0$
). Note that while this condition automatically ensures zero gradient of the chemical potential, it is actually a stronger condition than that in demanding all the relevant odd derivatives of
${\it\phi}$
to vanish separately. We do, however, find no difference in our numerical results between these. For method MI, the bounce-back boundary conditions assure we have no slip and no permeation of either fluid across the boundary.
Illustration of our definition of the slip length
${\it\xi}$
, defined as the position of the first local maximum of
$\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}$
as measured along the interface (defined as the locus
${\it\phi}=0$
).
2.5. Definition of the slip length
Our definition of the slip length uses the fact that the slip mechanism in the binary fluid model is intermolecular diffusion. In steady state, the Cahn–Hilliard equation of motion for the order parameter reads
$$\begin{eqnarray}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}=M{\rm\nabla}^{2}{\it\mu}.\end{eqnarray}$$
In figure 2, we plot the way in which
$\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}$
typically varies along the fluid–fluid interface (taken as the locus
${\it\phi}=0$
, as noted above). We take the slip length
${\it\xi}$
to be given by the distance from the wall to the local maximum. In measuring this typical distance away from the wall over which
$\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\phi}=M{\rm\nabla}^{2}{\it\mu}$
remains appreciable,
${\it\xi}$
characterises the length scale over which intermolecular diffusion is important.
Other definitions of the slip length have been used in the literature to characterize the diffusive mechanism. For example, Briant & Yeomans (Reference Briant and Yeomans2004) defined the slip length as the point at which
$M{\rm\nabla}^{2}{\it\mu}$
reaches
$-10$
% of the value at the wall after passing the first maximum, as measured along the interface; Yue et al. (Reference Yue, Zhou and Feng2010) suggested the use of the location of the stagnation point in the flow. As we will show in the results section, we recover the same scaling laws for the diffuse and sharp interface limits as in Briant & Yeomans (Reference Briant and Yeomans2004) and Yue et al. (Reference Yue, Zhou and Feng2010), respectively, strongly suggesting these definitions are not independent, as expected if there is a genuine characteristic length scale in the problem.
2.6. Parameters and dimensionless groups
For the model and flow geometry that we have now defined, we have the following parameters: the fluid density
${\it\rho}$
, the fluid viscosity
${\it\eta}$
, the modulus
$G$
that sets the scale of the free energy density, the mobility
$M$
, the diffuse interfacial width
$l$
, the box dimensions
$L_{x}$
and
$L_{y}$
and the applied flow velocity
$V_{0}$
. (Also relevant numerically are the spatial mesh size and the time step; we have ensured convergence for all the results presented.) As noted above, the fluid density is set to zero (inertialess flow) for methods MII and MIII, and small for MI. This leaves
${\it\eta},G,M,l,L_{y},L_{x}$
and
$V_{0}$
as the relevant parameters. We are then free to choose units of mass, length and time, thereby reducing by three the number of parameters that must be explored numerically. Accordingly, we work in what follows with the following dimensionless groups:
$\sqrt{M{\it\eta}}/l$
,
$l/L_{y}$
,
$L_{x}/L_{y}$
and the capillary number
$Ca={\it\eta}V_{0}/{\it\sigma}$
. The last of these characterises the speed of the imposed flow compared to the intrinsic interfacial velocity scale
${\it\sigma}/{\it\eta}$
formed by comparing the surface tension to the fluid viscosity (recall that the surface tension
${\it\sigma}\propto Gl$
). Among these we set the aspect ratio of the box
$L_{x}=2L_{y}$
throughout. The ratio
$l/L_{y}$
is always set small, to ensure that the results are not contaminated by finite box size effects. This leaves
$\sqrt{M{\it\eta}}/l$
and
$Ca$
as the two dimensionless groups to be explored numerically.
3. Numerical methods
In this section we describe our three numerical methods. For convenience we will refer to these throughout the manuscript as MI, MII and MIII. We emphasize that each was coded independently of the others, one by each of the authors of the manuscript, thereby providing confidence that our results are free of algorithmic artefacts.
3.1. Method I: lattice Boltzmann method
We use a standard free energy lattice Boltzmann method to solve the binary fluid equations of motion. The basic idea behind the lattice Boltzmann algorithm is to associate distribution functions,
$\boldsymbol{f}(\boldsymbol{r},t)=\{f_{i}(\boldsymbol{r},t)\}$
and
$\boldsymbol{g}(\boldsymbol{r},t)=\{g_{i}(\boldsymbol{r},t)\}$
, discrete in time and space, to a set of velocity directions
$\boldsymbol{e}_{i}$
. Here we use a D2Q9 model, where
$\boldsymbol{e}_{i}=(0,0),(\pm 1,0),(0,\pm 1),(\pm 1,\pm 1)$
. The physical variables are related to the distribution functions by
$$\begin{eqnarray}{\it\rho}=\mathop{\sum }_{i}f_{i},\quad {\it\rho}v_{{\it\alpha}}=\mathop{\sum }_{i}f_{i}e_{i{\it\alpha}},\quad {\it\phi}=\mathop{\sum }_{i}g_{i},\quad {\it\phi}v_{{\it\alpha}}=\mathop{\sum }_{i}g_{i}e_{i{\it\alpha}}.\end{eqnarray}$$
The time evolution equations for the distribution functions can be broken down into two steps. In the collision step, we have used a single relaxation time approach for the order parameter and a multiple relaxation time approach for the fluid density,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\boldsymbol{f}^{\prime }(\boldsymbol{r},t)=\boldsymbol{f}(\boldsymbol{r},t)-\unicode[STIX]{x1D648}^{-\mathbf{1}}\unicode[STIX]{x1D64E}\unicode[STIX]{x1D648}[\boldsymbol{f}(\boldsymbol{r},t)-\boldsymbol{f}^{eq}(\boldsymbol{r},t)],\\ \displaystyle \boldsymbol{g}^{\prime }(\boldsymbol{r},t)=\boldsymbol{g}(\boldsymbol{r},t)-\frac{1}{{\it\tau}_{{\it\phi}}}\left[\boldsymbol{g}(\boldsymbol{r},t)-\boldsymbol{g}^{eq}(\boldsymbol{r},t)\right].\end{array}\right\} & & \displaystyle\end{eqnarray}$$
As shown by Pooley et al. (Reference Pooley, Kusumaatmaja and Yeomans2008), the multiple relaxation time approach is significantly more reliable to reduce spurious velocities at the contact line.
$\boldsymbol{f}^{eq}$
and
$\boldsymbol{g}^{eq}$
are local equilibrium distribution functions. The choice of the local equilibrium functions must recover the correct thermodynamics and hydrodynamics equations of motion in the continuum limit. The matrix
$\unicode[STIX]{x1D648}$
performs a change of basis to more physically relevant variables, including the density, momentum and viscous stress tensor. The matrix
$\unicode[STIX]{x1D64E}$
is a diagonal matrix and contains the information about how fast each of these physical variables relaxes to equilibrium at every time step. Detailed expressions for
$\boldsymbol{f}^{eq}$
,
$\boldsymbol{g}^{eq}$
,
$\unicode[STIX]{x1D648}$
, and
$\unicode[STIX]{x1D64E}$
can be found in Pooley et al. (Reference Pooley, Kusumaatmaja and Yeomans2008).
The relaxation parameter
${\it\tau}$
for the viscous stress tensor in
$\unicode[STIX]{x1D64E}$
can be related to the kinematic shear viscosity
${\it\nu}={\it\eta}/{\it\rho}$
via
$$\begin{eqnarray}{\it\nu}=({\rm\Delta}t({\it\tau}-1/2))/3.\end{eqnarray}$$
Similarly, the relaxation parameter
${\it\tau}_{{\it\phi}}$
for the order parameter is related to the mobility parameter in the Cahn–Hilliard equation through
$$\begin{eqnarray}M={\rm\Delta}t{\it\Gamma}({\it\tau}_{{\it\phi}}-{\textstyle \frac{1}{2}}).\end{eqnarray}$$
${\it\Gamma}$
is a constant and
${\rm\Delta}t$
is the simulation time step. In practice, the relaxation times can therefore be tuned to match the desired continuum dynamical variables,
$M$
and
${\it\eta}$
.
In the propagation step, the updated distribution functions are passed on to the neighbouring lattice points,
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\boldsymbol{f}(\boldsymbol{r}+\boldsymbol{e}_{i}{\rm\Delta}t,t+{\rm\Delta}t)=\boldsymbol{f}^{\prime }(\boldsymbol{r},t),\\ \boldsymbol{g}(\boldsymbol{r}+\boldsymbol{e}_{i}{\rm\Delta}t,t+{\rm\Delta}t)=\boldsymbol{g}^{\prime }(\boldsymbol{r},t).\end{array}\right\} & & \displaystyle\end{eqnarray}$$
At the two walls we use a bounce-back rule for both the
$\{f_{i}(\boldsymbol{r},t)\}$
and
$\{g_{i}(\boldsymbol{r},t)\}$
distribution functions. These ensure boundary conditions of no slip, no permeation and no diffusion for the fluid material across the boundary (Ladd & Verberg Reference Ladd and Verberg2001).
3.2. Method II: spectral method
Within this method, at each numerical time step we separately (i) solve the hydrodynamic sector of the dynamics to update the fluid velocity field
$\boldsymbol{v}$
at fixed phase field
${\it\phi}$
, then (ii) update the phase field at fixed fluid velocity field. In part (i) we use a streamfunction formulation to ensure that the incompressibility condition is automatically satisfied. After eliminating the pressure by taking the curl of the generalised Stokes equation ((2.9) with the left-hand side set to zero) we solve for the streamfunction using a Fourier spectral method.
For the phase field dynamics in (ii) we use a third order upwind scheme (Pozrikidis Reference Pozrikidis2011) to update the convective term on the left-hand side of (2.7). The gradient terms on the right-hand side are solved using a Fourier spectral method using both sine and cosine modes in the periodic
$x$
direction and only cosine modes in the
$y$
direction for consistency with the imposed boundary conditions
$\partial _{y}{\it\phi}=\partial _{y}^{3}{\it\phi}=\partial _{y}^{5}{\it\phi}=0$
at the walls
$y=0,L_{y}$
.
For the time stepping we adopt a method originally proposed by Eyre (Reference Eyre1998), and later studied in depth by Guillén-González & Tierra (Reference Guillén-González and Tierra2013). This splits the free energy term into an expansive part,
$-{\rm\nabla}^{2}{\it\phi}$
, and a contractive part,
${\rm\nabla}^{2}({\it\phi}^{3})$
,
$-{\rm\nabla}^{4}{\it\phi}$
. These are then distributed between the
$n$
th and
$(n+1)$
th time steps as
$$\begin{eqnarray}\frac{{\it\phi}^{n+1}-{\it\phi}^{n}}{{\rm\Delta}t}=MG\left[\left(-l^{2}{\rm\nabla}^{4}{\it\phi}+2{\rm\nabla}^{2}{\it\phi}\right)^{n+1}+\left({\rm\nabla}^{2}{\it\phi}^{3}-3{\rm\nabla}^{2}{\it\phi}\right)^{n}\right],\end{eqnarray}$$
where
${\rm\Delta}t$
is the time step. This method permits significantly larger time steps than, for example, a Crank–Nicolson algorithm (Press et al.
Reference Press, Teukolsky, Vetterling and Flannery1992) in which all terms are split equally between the
$n$
th and
$(n+1)$
th time steps. When present, for
${\it\alpha}>0$
, the sixth-order gradient term is treated fully implicitly.
3.3. Method III: immersed boundary method
In methods I and II, the walls of the flow cell were included by means of a simulation box that is closed in the flow gradient direction
$y$
. In method III we instead use a biperiodic box that conveniently enables us to solve both the hydrodynamic sector of the dynamics (again in the incompressible streamfunction formulation at zero Reynolds number) and the diffusive part of the concentration dynamics using fast Fourier transforms in both spatial dimensions.
To incorporate the walls of the flow cell at
$y=0,L_{y}$
, we then include a set of immersed boundary forces using smoothed (Peskin) delta functions as source terms in the Stokes equation along the desired location of each wall, as discussed in Peskin (Reference Peskin2002) and Lai & Peskin (Reference Lai and Peskin2000). The force required to ensure zero velocity at the location of each delta function is then evolved as
$$\begin{eqnarray}\partial _{t}\boldsymbol{F}=-{\it\kappa}\left(\boldsymbol{V}-\boldsymbol{V}_{0}\right),\end{eqnarray}$$
where
$\boldsymbol{V}_{0}$
is the prescribed wall velocity. Here
${\it\kappa}$
is a spring constant, and all the results shown below have been converged to the limit
${\it\kappa}\rightarrow \infty$
in successive runs.
For the phase field, we assume the boundary conditions
$\partial _{y}{\it\mu}=0$
at the plates
$y=0,L_{y}$
. This is implemented in an analogous way to the no-slip boundary condition, by means of adding an extra source term contribution
$$\begin{eqnarray}\partial _{t}{\it\phi}=-K\hat{\boldsymbol{n}}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}{\it\mu}\end{eqnarray}$$
to the equation of motion for the concentration at the location of the Peskin delta functions, where
$\hat{\boldsymbol{n}}$
is a unit vector normal to the wall.
4. Results
We now present our results. We start in § 4.1 by considering the standard Landau free energy with
${\it\alpha}=0$
, before in § 4.2 appraising the robustness of our findings to the inclusion of the additional interfacial gradient (curvature) term in the free energy, setting
${\it\alpha}>0$
. In § 4.3, we compare our numerical results for the dynamic contact angles to the analytical predictions by Cox (Reference Cox1986).
4.1. Standard Landau theory
As discussed previously, in the existing literature two apparently contradictory scaling laws have been proposed for the contact line’s slip length
${\it\xi}$
. Several authors, including Jacqmin (Reference Jacqmin2000) and Yue et al. (Reference Yue, Zhou and Feng2010), have proposed that the slip length is proportional to the diffusive length scale
$l_{D}$
, which describes the length below which intermolecular diffusion dominates advection and above which the opposite is true, giving
${\it\xi}\propto l_{D}=(M{\it\eta})^{1/2}$
. In contrast, the lattice Boltzmann simulations and scaling arguments of Briant & Yeomans (Reference Briant and Yeomans2004) suggest that the slip length depends not only on this macroscopic diffusive length scale
$l_{D}$
, but also on the microscopic interfacial width
$l$
, with
${\it\xi}\propto (l_{D}l)^{1/2}$
.
To resolve this apparent discrepancy we have performed extensive numerical simulations of the moving contact line problem across four decades of the relevant dimensionless control parameter
$M{\it\eta}/l^{2}=l_{D}^{2}/l^{2}$
. Moreover, to ensure that our results are free of algorithmic artefacts, we have used (across the entire range) the three different, independently designed and coded, numerical techniques just described. We focus the discussion in this section on the case of planar Couette flow, returning to consider planar Poiseuille flow later in the manuscript. We also focus on the slow flow limit
$Ca\rightarrow 0$
, returning later on to consider the effects of finite
$Ca$
.
The results are shown in figure 3. As can be seen, the three numerical techniques give virtually indistinguishable results across all four decades of
$M{\it\eta}/l^{2}$
. Over these four decades we can distinguish two distinct regimes: a (i) diffuse interface regime when
$M{\it\eta}/l^{2}\ll 1$
, in which
${\it\xi}\propto (l_{D}l)^{1/2}$
, and (ii) a sharp interface regime when
$M{\it\eta}/l^{2}\gg 1$
, in which
${\it\xi}\propto l_{D}$
. In between these two distinct regimes is a broad cross-over window that itself spans around the decade either side of
$M{\it\eta}/l^{2}\approx 0.2$
. In this way, importantly, our results encompass and unify both of the two previously apparently contradictory scalings put forward in the literature, by showing each to apply in a different regime of the relevant dimensionless control parameter
$M{\it\eta}/l^{2}$
.
In addition to the two basic length scales
$l_{D}$
and
$l$
present in the model equations, out of which the slip length scale
${\it\xi}$
emerges in the manner just described, there are two other length scales present in our simulations, set by the system size
$L_{y}$
and the discretisation scale
${\rm\Delta}x$
. We have ensured that the results presented in figure 3 are independent of possible finite size effects due to
$L_{y}$
and discretisation errors due to
${\rm\Delta}x$
. The former is a particular hazard in the limit of large
$M{\it\eta}/l^{2}$
, where the slip length
${\it\xi}$
becomes large, and the latter in the limit of small
$M{\it\eta}/l^{2}$
, where the slip length becomes small.
Data for the slip length
${\it\xi}$
from all three numerical methods for planar Couette flow clearly demonstrating both the diffuse- and sharp-interface limits for small and large
$M{\it\eta}/l^{2}$
respectively. All the results here are in the slow flow regime,
$Ca<0.01$
, where the slip length is independent of the shear velocity.
Having shown our numerical results to capture both the sharp interface limit of
${\it\xi}\propto l_{D}$
for large
$M{\it\eta}/l^{2}$
, and the diffusive interface limit of
${\it\xi}\propto (l_{D}l)^{1/2}$
for small
$M{\it\eta}/l^{2}$
, we are now in a position to reprise carefully the scaling arguments put forward in the earlier literature for each of these two scaling forms separately, and to discuss the validity of the assumptions made in arriving at these functional forms.
We will follow a similar line of arguments as in Yue et al. (Reference Yue, Zhou and Feng2010), where the Stokes and steady-state Cahn–Hilliard equations are integrated in the
$x$
-direction along the wall. For concreteness, let us focus on the contact line at the top wall, i.e. at
$y=L_{y}$
. The same scaling analysis applies for the bottom wall (
$y=0$
). For the
$x$
-component of the Stokes equation, this gives
$$\begin{eqnarray}\int _{-\infty }^{\infty }\left[{\it\eta}\frac{\partial ^{2}v_{x}}{\partial y^{2}}+{\it\eta}\frac{\partial ^{2}v_{x}}{\partial x^{2}}-\frac{\partial p}{\partial x}-{\it\phi}\frac{\partial {\it\mu}}{\partial x}\right]\text{d}x=0.\end{eqnarray}$$
The second term can be integrated to
${\it\eta}\partial v_{x}/\partial x|_{-\infty }^{\infty }=0$
. The contribution of the third term in the integral is also zero, because the pressure attains the same, constant value on either side of the interface, far from the contact line. The fourth (final) term scales as
${\it\mu}_{max}$
, which represents the magnitude of the chemical potential for a given flow set-up.
The first term in (4.1) needs to be analysed more carefully. In the sharp interface limit, the term
$\partial ^{2}v_{x}/\partial y^{2}$
scales as
$V_{0}/{\it\xi}^{2}$
, and this term is significant across a length scale given by the slip length
${\it\xi}$
. Thus, we have
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\it\eta}\frac{V_{0}}{{\it\xi}^{2}}{\it\xi}\sim {\it\mu}_{max},\\ {\it\mu}_{max}\sim {\it\eta}V_{0}/{\it\xi}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
In the diffuse interface limit, the fluid velocity still varies across a length scale
${\it\xi}$
in the
$y$
-direction, and as such, the term
$\partial ^{2}v_{x}/\partial y^{2}$
still scales as
$V_{0}/{\it\xi}^{2}$
. However, to capture the key physics, the spatial window of integration must be broadened to the interface width
$l$
, since here
${\it\xi}$
is less than the interfacial width
$l$
: recall the bottom leftmost data points in figure 3. Taking the suitable of window of integration into account, (4.1) leads to
$$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\it\eta}\frac{V_{0}}{{\it\xi}^{2}}l\sim {\it\mu}_{max},\\ {\it\mu}_{max}\sim {\it\eta}V_{0}l/{\it\xi}^{2}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$
Let us now carry out a similar analysis for the Cahn–Hilliard equation,
$$\begin{eqnarray}\int _{-\infty }^{\infty }\left[v_{x}\frac{\partial {\it\phi}}{\partial x}+v_{y}\frac{\partial {\it\phi}}{\partial y}-M\frac{\partial ^{2}{\it\mu}}{\partial x^{2}}-M\frac{\partial ^{2}{\it\mu}}{\partial y^{2}}\right]\text{d}x=0.\end{eqnarray}$$
The first term can be integrated to give
$V_{0}\,{\it\phi}|_{-\infty }^{\infty }=2V_{0}$
. The second term is zero simply because
$v_{y}=0$
at
$y=0$
. The third term integrates to
$-M\partial {\it\mu}/\partial x|_{-\infty }^{\infty }$
, which is again zero because
$\partial {\it\mu}/\partial x=0$
far from the interface.
The result for the fourth term depends on the slip length regimes.
$\partial ^{2}{\it\mu}/\partial y^{2}$
scales as
${\it\mu}_{max}/{\it\xi}^{2}$
. However, the window of integration is different in the sharp (
${\it\xi}$
) and diffuse (
$l$
) interface limits. For the sharp interface limit, the Cahn–Hilliard equation leads to the scaling
$$\begin{eqnarray}V_{0}\sim M{\it\mu}_{max}/{\it\xi}.\end{eqnarray}$$
On the other hand, for the diffuse interface limit, we obtain
$$\begin{eqnarray}V_{0}\sim M{\it\mu}_{max}l/{\it\xi}^{2}.\end{eqnarray}$$
We are now in the position to derive the scaling law predictions for the slip length. Combining (4.2) and (4.5) gives
$$\begin{eqnarray}{\it\xi}\sim (M{\it\eta})^{1/2}=l_{D}\end{eqnarray}$$
in the sharp interface regime. Similarly, combining (4.3) and (4.6) for the diffuse interface regime results in
$$\begin{eqnarray}{\it\xi}\sim (M{\it\eta}l^{2})^{1/4}=(l_{D}l)^{1/2}.\end{eqnarray}$$
As seen in figure 3, in between the sharp and diffuse interface regimes is a smooth cross-over where the behaviour varies smoothly from one scaling law to another. We define the cross-over point by fitting the power laws in the sharp and diffuse interface regimes separately, and finding where these two fits intersect each other. This occurs at
$l_{D}/l=0.38$
, in broad agreement with the sharp interface criterion
$l_{D}/l>0.25$
proposed by Yue et al. (Reference Yue, Zhou and Feng2010)
We note finally that the values of the slip length presented in figure 3 are independent of the wall velocity, and correspondingly of the Capillary number
$Ca={\it\eta}V_{0}/{\it\sigma}$
, provided
$Ca\lesssim 0.01$
(data not shown). We shall demonstrate this independence explicitly below for the case of planar Poiseuille flow (see the discussion associated with figure 5 below).
4.2. Influence of a curvature term
In this section we test the extent to which our results for the slip length do or do not depend on the microscopic details of the diffuse interface model. Intuitively we might expect the slip length to be independent of microscopics in the sharp interface regime, but dependent on microscopics in the diffuse interface regime. Our numerical results below provide evidence in favour of this intuition, to which we shall however also add a note of caution at the end of this section.
Numerical data demonstrating how slip lengths in the diffuse-interface limit are affected by the inclusion of a curvature free energy term. The results of the standard Landau theory (
${\it\alpha}=0$
; red circles) are shown for comparison. Inset: percentage change relative to the case without curvature for
${\it\alpha}=0.1$
(green squares),
${\it\alpha}=1$
(blue diamonds). These data are obtained using the spectral method (MII) in planar Couette flow.
There are many possible ways to modify the basic Landau theory for which we presented results in the previous subsection. We focus here on the simple but non-trivial extension made by introducing the additional interfacial curvature term with a strength set by
${\it\alpha}$
in (2.5). Figure 4 shows our numerical results for
${\it\alpha}=0.1$
and 1.0. Shown in the same plot for comparison are our original results, already presented in the previous subsection, for the original theory with
${\it\alpha}=0$
. As can be seen, the introduction of the curvature term affects the slip length quite strongly in the diffuse interface regime. This is to be expected: in this regime the physics of the problem is determined not only by macroscopic quantities, but by the microscopic gradient contributions to the free energy. In contrast, as the control parameter
$M{\it\eta}/l^{2}$
increases into the sharp interface regime the dependence of the slip length on this microscopic parameter
${\it\alpha}$
dramatically decreases and appears to become negligible in the sharp interface limit
$M{\it\eta}/l^{2}\rightarrow \infty$
.
Although these numerical results strongly suggest that the slip length is independent of the microscopic (gradient) contributions to the free energy in the sharp interface regime, we attach the following note of caution. In the case of equilibrium phase coexistence in the absence of flow, it is possible to show by integrating across the interface between the phases that macroscopic bulk quantities, such as the value of the chemical potential at phase coexistence, are independent of the structure of that interface, as determined by the spatial gradient terms in the free energy. However, the same is not true in general for the case of phase coexistence under the non-equilibrium conditions of an applied flow because the equations of motion are usually non-integrable in that case: they cannot be integrated across the interface to give a result that depends only on the bulk quantities on either side of the interface, but the result instead retains a dependence on the structure of the interface itself. This issue has been discussed in particular detail in the case of stress selection in shear banded flows (Olmsted & David Lu Reference Olmsted and Lu1999; Lu, Olmsted & Ball Reference Lu, Olmsted and Ball2000). Therefore, while our numerical results do strongly suggest that the slip length is independent of the gradient terms in the sharp interface limit, this question merits some further attention in future studies.
4.3. Relation to Cox’s result
The results presented so far have demonstrated that, in the sharp interface limit, the emergent slip length is determined by the macroscopic dynamical variables
${\it\eta}$
and
$M$
. While this is clearly an important finding in terms of our physical understanding of the moving contact line problem, in numerical practice this sharp interface limit
$M{\it\eta}/l^{2}\rightarrow \infty$
can be difficult to attain. Indeed, the time taken to attain a steady state after the inception of flow increases dramatically with increasing
$M{\it\eta}/l^{2}$
. Lengthy simulations were required to obtain the rightmost data points in figure 3, with the extreme cases requiring run times of the order of weeks with a parallelised code using 4 cores of Intel(R) Xeon(R) 2.40 GHz.
With this in mind, we now turn to address two questions of practical numerical importance: (i) whether simulations carried out in the diffuse interface regime still reproduce the expected macroscopic dynamics far from the contact line, and (ii) whether our definition of the slip length is still meaningful outside the sharp interface regime. We shall address these questions in the context of the Poiseuille flow protocol sketched in figure 1(b), in particular asking whether the dependence of the observed dynamic contact angle on flow velocity is the same in the diffuse interface regime as in the sharp interface regime.
Numerical data for the Poiseuille flow protocol, for
$M{\it\eta}/l^{2}=0.01$
(black circles ○), 0.1 (red squares ▫),
$1$
(green diamonds ♢),
$4$
(blue triangles ▵). (a) Demonstration that the slip length converges to a constant value as
$Ca\rightarrow 0$
. The inset shows that the measured value of
${\it\xi}_{Ca\rightarrow 0}$
using the Poiseuille flow protocol is in good agreement with the best fit power laws obtained for shear flow in figure 3. (b) Symbols mark the dynamic contact angle
${\it\theta}_{d}$
as measured from simulation. The solid lines show the analytical predictions of Cox (Reference Cox1986) where
${\it\delta}={\it\beta}{\it\xi}_{Ca}/L_{y}$
. The inset shows the best fit values of
${\it\beta}$
which appear to converge
${\it\beta}\rightarrow 1$
in the sharp interface limit,
$M{\it\eta}/l^{2}\rightarrow \infty$
. The simulation data are obtained using the spectral method (MII).
We start by checking that the slip length obtained for the Poiseuille flow protocol is equivalent to that for planar Couette flow. Figure 5(a) shows that the slip length for the Poiseuille case converges to a constant value as
$Ca\rightarrow 0$
, as for the planar Couette case. We have checked in this regime that the measured values of the slip length for the Couette and Poiseuille flows are within 0.5 % of each other over the whole range of
$M{\it\eta}/l^{2}$
. See the inset of figure 5(a), where the slip lengths for the Poiseuille flow are compared to the best fit power laws extracted for planar Couette flow in both the diffuse and sharp interface regimes. For larger
$Ca>0.01$
the slip length starts to depend on the fluid velocity, as in planar Couette flow.
We now turn to the dynamic contact angle, as defined in figure 1(b). It is known that this angle increases with increasing flow speed: for increasing pressure drop
${\rm\Delta}P$
in figure 1(b) the interface is deformed more significantly due to the larger flow speed mid-channel compared to the stationary wall, and correspondingly the dynamic contact angle deviates more from the equilibrium value. Associated with this is a larger gradient in the chemical potential, which indeed allows the contact line to move at a higher velocity.
An analytical prediction for this increase of the dynamic contact angle with flow speed, as characterised by the capillary number
$Ca$
, was given by Cox (Reference Cox1986):
$$\begin{eqnarray}g({\it\theta}_{d})=g({\it\theta}_{e})+Ca\ln ({\it\delta}^{-1}).\end{eqnarray}$$
Here
${\it\theta}_{d}$
is the dynamic contact angle and
${\it\theta}_{e}$
is the contact angle at equilibrium (here 90°). For fluids of matched viscosity the function
$g$
is defined as
$$\begin{eqnarray}g({\it\theta})=\int _{0}^{{\it\theta}}\frac{{\rm\pi}{\it\phi}({\rm\pi}-{\it\phi})+(2{\rm\pi}{\it\phi}-{\rm\pi}^{2})\sin {\it\phi}\cos {\it\phi}-{\rm\pi}\sin ^{2}{\it\phi}}{2{\rm\pi}^{2}\sin {\it\phi}}\,\text{d}{\it\phi}.\end{eqnarray}$$
The parameter
${\it\delta}$
in (4.9) is the ratio of the microscopic slip length to some characteristic macroscopic length scale of the system. We choose that to be the channel width
$L_{y}$
and set
${\it\delta}={\it\beta}{\it\xi}_{Ca}/L_{y}$
where
${\it\beta}$
is a constant fitting parameter.
To test Cox’s prediction, we plot our simulation results for the dependence
${\it\theta}_{d}$
on the capillary number
$Ca$
in figure 5(b). We do so for four different values of
$M{\it\eta}/l^{2}$
, which we ensured span both the diffuse and sharp interface limits, as well as the cross-over regime between them. Pleasingly, we find excellent agreement with Cox’s prediction across the the full range of
$Ca$
, for all values of
$M{\it\eta}/l^{2}$
. Furthermore, as shown in the inset of figure 5(b), the fitting parameter
${\it\beta}$
appears suggestively to converge to unity in the sharp interface limit,
$M{\it\eta}/l^{2}\rightarrow \infty$
.
These results suggest that simulations in the diffuse interface limit (and similarly the cross-over regime) still provide a reliable representation of the macroscopic dynamics. However, the effective slip length must be suitably corrected by an
$O(1)$
prefactor
${\it\beta}$
. Increasingly minor corrections are needed as one approaches the sharp interface limit, and we expect Cox’s result would be directly obtained for
$M{\it\eta}/l^{2}\rightarrow \infty$
, where
${\it\beta}\rightarrow 1$
.
5. Conclusions
In this paper, we have performed extensive numerical simulations to reconcile previously apparently contradictory scaling laws for the slip length associated with a moving contact line in diffuse interface models. We did so by exploring four decades of the relevant dimensionless control parameter
$M{\it\eta}/l^{2}=l_{D}^{2}/l^{2}$
, and three decades of the capillary number
$Ca$
. Moreover we used three independent numerical methods, demonstrating their predictions to be in excellent quantitative agreement across all decades.
In the limit of slow flows corresponding to typical capillary numbers
$Ca<0.01$
, we find that the slip length converges to a well-defined value that is independent of the flow speed. In this slow-flow regime we explored four decades in
$M{\it\eta}/l^{2}=l_{D}^{2}/l^{2}$
. Doing so enabled us to distinguish two distinct regimes for the scaling of the slip length, and a broad cross-over window in between. For
$M{\it\eta}/l^{2}\ll 1$
we found a diffuse interface regime in which the slip length scales as
${\it\xi}\sim (l_{D}l)^{1/2}$
, consistent with the previous result of Briant & Yeomans (Reference Briant and Yeomans2004), and is furthermore sensitive to the form of the microscopic interfacial terms in the free energy functional. In contrast for
$M{\it\eta}/l^{2}\gg 1$
we recover the sharp interface limit previously discussed by Jacqmin (Reference Jacqmin2000) and Yue et al. (Reference Yue, Zhou and Feng2010). Here, the slip length is proportional to the diffusion length scale,
${\it\xi}\propto l_{D}\sim (M{\it\eta})^{1/2}$
, depending only on the macroscopic dynamical variables of viscosity
${\it\eta}$
and molecular mobility
$M$
. Our numerical results further suggest the slip length to be insensitive to changes in the microscopic interfacial free energy terms in this sharp interface regime.
For
$Ca>0.01$
, the effective slip length decreases with increasing capillary number, in qualitative agreement with the work of Chen et al. (Reference Chen, Jasnow and Viñals2000). However, we are not able to reproduce their prediction of
${\it\xi}\propto V_{0}^{-1/2}$
in this fast flowing regime. More broadly, so far we are unable to find any reliable scaling law that is valid across a broad range of Capillary number, and for both the sharp and diffuse interface regimes.
Our numerical results also allow us to appraise suitable sets of simulation parameters. By comparing our simulation results for the dynamic contact angle to the analytical prediction by Cox (Reference Cox1986), we find excellent agreement if we allow the slip length to be rescaled by a dimensionless prefactor
${\it\beta}=O(1)$
, which pleasingly appears to converge to unity in the sharp interface limit. In numerical practice, however, the sharp interface limit can be very expensive to realise, partly because the time taken to attain a steady state increases dramatically with increasing
$M{\it\eta}/l^{2}$
, and partly because the slip length is
$O(10)$
larger than the interface width, with demanding implications on the overall size of the simulation box. In the diffuse interface limit, in contrast, the time scale and length scale requirements are less stringent, and simulations accordingly much less expensive to perform. Our fits to the Cox formula suggest that a reliable representation of the macroscopic physics is nonetheless attainable even in the diffuse regime, provided the lengths are suitable rescaled by this corrective prefactor
${\it\beta}=O(1)$
.
There are several avenues for future work. Firstly, it would be interesting to investigate how the two scaling regimes for contact line slip in diffuse interface models affect more complex flow problems, such as the macroscopic dynamics of a liquid droplet (Servantie & Muller Reference Servantie and Muller2008; Mognetti, Kusumaatmaja & Yeomans Reference Mognetti, Kusumaatmaja and Yeomans2010; Thampi, Adhikari & Govindarajan Reference Thampi, Adhikari and Govindarajan2013) or the criterion for the liquid entrainment (wetting failure) (Eggers Reference Eggers2004; Sbragaglia, Sugiyama & Biferale Reference Sbragaglia, Sugiyama and Biferale2008), and whether similar renormalization of the slip length is adequate to account for their differences, as for the Poiseuille flow considered here. In these cases, we need to implement non-neutral wetting boundary conditions. We expect the prefactors of the scaling laws to change with contact angles, but not the exponents themselves.
Secondly, given the cross-over we have reported in this study, it would be insightful to revisit the liquid–gas system and understand the sharp interface limit of this model. In the binary model the slip mechanism is diffusion, and as such, the slip length depends on the mobility parameter, in addition to the fluid viscosity and the interface width. In the liquid gas system, the mechanism is evaporation–condensation instead, and we expect the slip length to depend on the density ratio between the liquid and gas phases.
Thirdly, we note that we have taken the zero temperature limit in this work and have not accounted for thermal fluctuations. It will be interesting to study if such fluctuations simply broaden the effective interface width, or they give rise to new phenomena. For example, enhanced spreading due to thermal fluctuations has been reported by Davidovitch, Moro & Stone (Reference Davidovitch, Moro and Stone2005) and Gross & Varnik (Reference Gross and Varnik2014).
Finally, experimental systems consisting of colloid–polymer mixtures are known to phase separate into colloid-rich and colloid-poor domains, with an interface width that is of order
$1~{\rm\mu}\text{m}$
(Aarts et al.
Reference Aarts, Dullens, Lekkerkerker, Bonn and van Roij2004; Setu et al.
Reference Setu, Dullens, Hernndez-Machado, Pagonabarraga, Aarts and Ledesma-Aguilar2015), much larger than the typical value for molecular fluids. Such a system can potentially be exploited to realise the different contact line slip regimes discussed here experimentally.
Acknowledgements
The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 279365.
