2. Model and numerical method

We investigate the local flows in the vicinity of the moving contact line, where the interface of fluid A and fluid B intersects the solid substrate with a microscopic contact angle $\theta _e$ , as sketched in figure 1. We adopt the Navier-slip condition to relieve the contact-line singularity and concentrate on the small region with significant slip effect. The motion of the fluids can be treated as a locally two-dimensional Stokes flow where inertial and gravitational forces are neglected. The material parameters of relevance include the viscosities of the fluids, ${\rm \mu} _A$ and ${\rm \mu} _B$ , and the interface tension $\gamma$ . It is convenient to use a Cartesian reference frame in which the interface evolves into a stationary state, while the substrate is moving at a constant velocity $U$ . The velocity $\boldsymbol{u}=(u_x,u_y)$ and pressure $p$ in each phase are governed by the Stokes equations

(2.1) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \quad -\boldsymbol{\nabla }p + {\rm \mu} {\nabla} ^2 \boldsymbol{u}=\boldsymbol{0}. \end{equation}

Figure 1.

Sketch of the geometry near a moving contact line. The computational domain is shown as the grey region.

The velocity satisfies the Navier-slip boundary condition at the substrate

(2.2) \begin{equation} u_x-U=\lambda \frac {\partial u_x}{\partial y}, \quad u_y=0, \end{equation}

where $\lambda$ is the slip length. At the fluid–fluid interface, the boundary conditions are

(2.3) \begin{align}&\qquad\,\, \boldsymbol{u}_A=\boldsymbol{u}_B, \end{align}

(2.4) \begin{align}& \gamma \kappa \boldsymbol{n}=(\boldsymbol{\sigma }_A-\boldsymbol{\sigma }_B)\boldsymbol{\cdot }\boldsymbol{n}. \end{align}

Here, $\boldsymbol{\sigma } = -p \boldsymbol{I} + {\rm \mu} [\boldsymbol{\nabla }\boldsymbol{u} +(\boldsymbol{\nabla }\boldsymbol{u})^T ]$ is the stress tensor, $\kappa =\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{n}$ is the curvature of the interface and $\boldsymbol{n}$ is the unit normal pointing into fluid A. Equations (2.3) and (2.4) represent the continuity of the velocity and tangential stress, while the normal stress has a jump due to the presence of the interface tension. For a general moving or deforming interface, we also require that there is no mass flow across the interface

(2.5) \begin{equation} \frac {\mathrm{d} \boldsymbol{x}_i}{\mathrm{d}t}=(\boldsymbol{u}_A\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}=(\boldsymbol{u}_B\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}, \end{equation}

where $\boldsymbol{x}_i$ represents the location of the interface and $t$ is time. Finally, a contact angle condition is needed at the contact line. For simplicity, we assume that the microscopic contact angle is independent of the contact-line speed and equals the equilibrium contact angle $\theta _e$ .

In addition to $\theta _e$ , there exist two dimensionless numbers governing the problem, i.e. the capillary number $Ca={\rm \mu} _A U/\gamma$ and the viscosity ratio $R={\rm \mu} _B/{\rm \mu} _A$ . We focus on the situations of $Ca\ll 1$ such that the interface is weakly deformed.

Before presenting an approximate theory of moving contact lines, we mention here the computational model to obtain accurate solutions of the full problem, which is important in validating the theoretical results. The flow domain of interest is truncated at a position sufficiently far from the contact line. We employ a semicircular region centred at the contact line as the computational domain, as shown in figure 1. Since the ‘far-field’ stress should be much smaller than that near the contact line, we set the stress-free condition at the new boundary, i.e. $\boldsymbol{\sigma }\boldsymbol{\cdot }\boldsymbol{n}=0$ . We also impose a vanishing curvature of the interface at the intersection $P$ . These conditions are expected to have no influence on the local wetting flow as long as the domain size is sufficiently large. To obtain a steady solution, the position of $P$ is fixed. The size of the domain can be determined by adjusting the distance between $P$ and the substrate, denoted by $H$ .

The governing equations in (2.1) together with the boundary conditions are solved using a boundary element method (BEM). This method has exhibited the ability to realise full-scale simulations of dynamic wetting with extremely small slip lengths, as demonstrated recently in bubble spreading (Ming et al. Reference Ming, Qin and Gao2023) and dip coating (Zhang & Gao Reference Zhang and Gao2024). The details of the numerical algorithm can be found in Ming et al. (Reference Ming, Qin and Gao2023). Starting from a flat interface with an inclination angle $\theta _e$ , the flow eventually evolves into a steady state. Figure 2 shows the steady solution of the interface slope angle $\theta$ as a function of the rescaled arc length $s/\lambda$ . Three representative sizes of the computational domain are selected, corresponding to $H/ \lambda =5\times 10^3$ , $1\times 10^4$ and $2\times 10^4$ . It is seen that the numerical results away from the outer boundary can hardly be distinguished as $H / \lambda$ increases. The numerical results presented hereafter were calculated at $H / \lambda =1\times 10^{4}$ , giving reliable comparison at $s/\lambda$ up to $O(10^3)$ .

Figure 2.

Variation of the slope angle $\theta$ as a function of $s/\lambda$ at $R=0$ , $\theta _e = 60^{\circ }$ , $Ca=0.01$ and three representative values of the domain size.

3. Generalised lubrication theory based on slippery wedge flows

The original form of the GL equation (Snoeijer Reference Snoeijer2006; Chan et al. Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013) relies on the approximation of the local flow by the no-slip wedge flow (Huh & Scriven Reference Huh and Scriven1971). We present alternatively a similarity solution of a slippery wedge flow, which can be used to approximate the local flow and to further develop an explicit form of the GL theory.

3.1. Slippery wedge flow

As a first approximation, the local geometry can be modelled by a two-phase wedge flow, in which the interface is assumed to be straight and tilted at an angle $\theta$ . Introducing the polar coordinates $(r,\phi )$ , the interface is given by $\phi = 0$ and the substrate is given by $\phi = -\theta$ and $\pi -\theta$ . In the presence of a Navier slip, Hocking (Reference Hocking1977) showed that closed-form solutions are not available. In contrast, we propose a similarity solution of the wedge flow by introducing a linearly varying slip length

(3.1) \begin{equation} \lambda =\beta r, \end{equation}

where $\beta$ is constant. Since this modified condition eliminates the sole length scale in the problem, a similarity solution in analogous to the no-slip case is retrieved. We note that the similarity solution is not the exact solution with a constant slip length and hence is not expected to describe the global flow near a moving contact line; nevertheless, we will show later that the local flow can be well approximated by this solution with $\beta$ appropriately adjusted using the radius of the local wedge (denoted by $r_s$ in figure 3), while the slip length remains fixed.

Figure 3.

Comparison of the velocity profile at $R=1$ , $Ca=0.005$ , $r_s/\lambda =1$ and $\theta _e=60^{\circ }$ in (a) and $120^\circ$ in (b). The slippery wedge flow represents (3.7) and (3.8) with $\beta =\lambda /r_s=1$ and the Huh & Scriven (Reference Huh and Scriven1971) solution corresponds to $\beta =0$ .

It is convenient to introduce the streamfunction $\varPsi$ such that

(3.2) \begin{equation} {u_r} = \frac {1}{r}\frac {\partial \varPsi }{\partial \phi },\quad {u_\phi } = - \frac {\partial \varPsi }{\partial r}. \end{equation}

Then the Stokes equations can be transformed into

(3.3) \begin{equation} \nabla^{4}\varPsi = 0, \end{equation}

which is complemented by the boundary conditions

(3.4) \begin{align}&\qquad\quad\,\, \frac {1}{r}\frac {\partial \varPsi _A}{\partial \phi } = U + \frac {\beta }{r}\frac {\partial ^2\varPsi _A }{\partial \phi ^2},\quad \varPsi _A = 0\quad \mbox{on}\quad \phi =-\theta , \end{align}

(3.5) \begin{align}&\qquad\, -\frac {1}{r}\frac {\partial \varPsi _B}{\partial \phi } = U + \frac {\beta }{r}\frac {\partial ^2\varPsi _B}{\partial \phi ^2},\quad \varPsi _B = 0\quad \mbox{on}\quad \phi =\pi -\theta , \end{align}

(3.6) \begin{align}& \varPsi _A = \varPsi _B = 0, \quad \frac {\partial \varPsi _A }{\partial \phi } = \frac {\partial \varPsi _B}{\partial \phi },\quad \frac {\partial ^2\varPsi _A}{\partial \phi ^2}= R\frac {\partial ^2\varPsi _B}{\partial \phi ^2}\quad \mbox{on}\quad \phi =0. \\[9pt] \nonumber \end{align}

The solution of (3.3) has the form

(3.7) \begin{align} \varPsi _A&=Ur\left [ M_A\phi \sin (\phi + \theta ) + N_A( \phi + \theta )\sin \phi \right ]\!,\end{align}

(3.8) \begin{align} \varPsi _B&=Ur\left [M_B\phi \sin (\phi + \theta -\pi ) + N_B(\phi + \theta -\pi )\sin \phi \right ]\!. \end{align}

The coefficients are determined by the boundary conditions (3.4)–(3.6) as

(3.9) \begin{equation} \left . \begin{split} M_A&=\{R[ f_1(\pi -\theta ) - (2\beta +\pi )f_2(\pi -\theta )+2\beta (\pi -2\theta )\sin ^2\theta ]\\ &\quad - f_1(\pi -\theta )-2\beta f_2(\pi -\theta ) \}/\unicode{x1D6E5} ,\\[2pt] N_A&= -\cos \theta \{ R [2\beta f_2(\pi -\theta )-f_3(\pi -\theta )-(4\beta +\pi )\tan \theta -4\beta (\pi -\theta )] \\ &\quad -f_1(\pi -\theta )-2\beta f_2(\pi -\theta )\}/\unicode{x1D6E5} ,\\[2pt] M_B&= \{f_1(\theta )-(2\beta +\pi )f_2(\theta )-2\beta (\pi -2\theta )\sin ^2\theta -R[f_1(\theta ) + 2\beta f_2(\theta ) ] \}/\unicode{x1D6E5} ,\\[2pt] N_B&= \cos \theta \{2\beta f_2(\theta )-f_3(\theta ) + (4\beta +\pi )\tan \theta -4\beta \theta - R[f_1(\theta )+2\beta f_2(\theta )]\}/\unicode{x1D6E5} , \end{split} \right \} \end{equation}

where

(3.10) \begin{align} f_1(\theta )&=\theta ^2-\sin ^2\theta ,\quad f_2(\theta )=\theta -\sin \theta \cos \theta ,\quad f_3(\theta )=\theta (\pi -\theta )+\sin ^2\theta , \end{align}

(3.11) \begin{align} \unicode{x1D6E5} &= R\left [f_1(\theta )+ 2\beta f_2(\theta ) \right ]\big[ f_2(\pi -\theta )+2\beta \sin ^2\theta \big] \nonumber \\ &\quad + \left [ f_1(\pi -\theta )+2\beta f_2(\pi -\theta ) \right ]\big[f_2(\theta )+2\beta \sin ^2\theta \big]. \end{align}

It can be easily verified that this solution reduces to the no-slip wedge flow given by Huh & Scriven (Reference Huh and Scriven1971) when $\beta =0$ .

Figure 4.

Comparison of the velocity profile at $R=1$ , $Ca=0.005$ and $\theta _e=60^{\circ }$ . The slippery wedge flow represents (3.7) and (3.8), and the Huh & Scriven (Reference Huh and Scriven1971) solution corresponds to $\beta =0$ . The solid lines are produced by BEM with $r_s/\lambda =1$ (a, b), 10 (c, d) and 100 (e, f). Here, $\phi$ has been shifted such that the substrate is described by $\phi =0$ and $180^\circ$ .

To check how the flow can be approximated by this slippery wedge flow, a comparison with the BEM solution is performed. For an arbitrary point on the interface, a locally tangent wedge can be defined, as shown in figure 3. We then draw a semicircle centred at the wedge tip, which in general is different from the contact line. Denoting the radius of the semicircle as $r_s$ , the value of $\beta$ is calculated as $\lambda /r_s$ . This procedure enables us to plot the velocity of the slippery and no-slip wedge flows along the semicircle and to compare with the BEM results. In figure 3, we have selected a point with $r_s=\lambda$ , where the slip effect is important. For both $\theta _e=60^\circ$ and $120^\circ$ , the slippery wedge flow agrees well with the numerical results, which we believe is attributable to the matching of the local slip length, while significant deviations of the Huh–Scriven solution occur. A more quantitative comparison is illustrated in figure 4, which shows the variation of $u_r$ and $u_\phi$ as functions of $\phi$ for $\theta _e=60^\circ$ . We choose three representative locations corresponding to $r_s/\lambda =1$ , 10 and 100. It is seen that the slippery wedge flow captures satisfactorily the BEM results, especially for the radial component of the velocity $u_r$ . When $r_s/\lambda \gg 1$ , the slip effect can be neglected and the curves for the slippery wedge flow, the Huh-Scriven flow and the numerical results collapse.

Furthermore, we examine the asymptotic behaviour of the slippery wedge flow as the moving contact line is approached. On the one hand, when $\beta$ is fixed, it can be shown that both the shear stress $\tau _{r\phi }$ and the pressure $p$ have a $r^{-1}$ singularity as $r\to 0$ due to the vanishing slip length. On the other hand, for the flow with the local approximation, we have $\lambda =\beta r_s$ fixed and hence $\beta$ diverges since $r_s\approx r$ near the contact line. In this situation, it can be demonstrated that $\varPsi \propto Ur/\beta$ and $\boldsymbol{u}\propto U/\beta$ , both of which vanish at the contact line as expected. Accordingly, the viscous stress $\tau _{r\phi }\propto {\rm \mu} U/(\beta r)\approx {\rm \mu} U/\lambda$ and is regularised. The pressure gradient can be evaluated from the momentum equation as $\partial p/\partial r\propto {\rm \mu} U/(\beta r^2)\approx {\rm \mu} U/(\lambda r)$ , yielding $p-p_0\propto ({{\rm \mu} U}/{\lambda })\ln ({r}/{\lambda })$ , where $p_0$ is a reference pressure; therefore, the pressure still exhibits a logarithmic singularity, which has been known for moving contact lines with a Navier-slip condition (Buckingham, Shearer & Bertozzi Reference Buckingham, Shearer and Bertozzi2003; Benilov & Benilov Reference Benilov and Benilov2015; Sibley et al. Reference Sibley, Nold, Savva and Kalliadasis2015).

3.2. Generalised lubrication equation

We can now derive a novel form of the GL equation, following the procedure of Snoeijer (Reference Snoeijer2006) and Chan et al. (Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013), but using instead the slippery wedge solution to approximate the local flow. The interface profile is described by $y=h(s)$ . Introducing the wedge approximation as in figure 3, we have $r_s=h/\sin \theta$ and hence $\beta =\lambda /r_s=\lambda \sin \theta /h$ . The velocity field of this wedge flow can be used to calculate the local stress on the interface, yielding a GL equation for the interface geometry. To this end, the normal component of the dynamic boundary condition (2.4) can be written as

(3.12) \begin{equation} \gamma \frac {\mathrm{d}\theta }{\mathrm{d}s}=p_B-p_A, \end{equation}

where we have used $\kappa =\mathrm{d}\theta /\mathrm{d}s$ , and the viscous contribution on the stress can be proved to vanish, similar to the Huh–Scriven flow. Differentiating (3.12) with respect to $s$ and using the momentum equation in (2.1), we have

(3.13) \begin{equation} \gamma \frac {\mathrm{d}^2\theta }{\mathrm{d}s^2}=\frac {\partial p_B}{\partial s}-\frac {\partial p_A}{\partial s}=\boldsymbol{e}_s\boldsymbol{\cdot }\left (\boldsymbol{\nabla }p_B-\boldsymbol{\nabla }p_A\right )=\boldsymbol{e}_s\boldsymbol{\cdot }\left ({\rm \mu} _B\boldsymbol{\nabla} ^2 \boldsymbol{u}_B-{\rm \mu} _A\boldsymbol{\nabla} ^2 \boldsymbol{u}_A\right ), \end{equation}

where $\boldsymbol{e}_s$ denotes the unit vector tangent to the interface. The right-hand side of this equation can be calculated using the solution of the slippery wedge flow with $\beta$ given above. In the local polar coordinates

(3.14) \begin{equation} \boldsymbol{e}_s\boldsymbol{\cdot }{\nabla} ^2 \boldsymbol{u}={\nabla} ^2u_r-\frac {u_r}{r^2}-\frac {2}{r^2}\frac {\partial u_\phi }{\partial \phi }=-2MU\frac {\sin \theta }{r^2}=-2MU\frac {\sin ^3\theta }{h^2}. \end{equation}

Accordingly, (3.13) becomes

(3.15) \begin{equation} \frac {\mathrm{d}^2\theta }{\mathrm{d}s^2}=\frac {2Ca\sin ^3\theta }{h^2}(M_A+RM_B). \end{equation}

After substituting (3.9) and making some elementary manipulations, we obtain the GL equation

(3.16) \begin{equation} \frac {\mathrm{d}^2\theta }{\mathrm{d}s^2} = - \frac {\textit{Ca}}{h}\frac {D(\theta ,R)h\sin ^2\theta + E(\theta ,R)\lambda \sin ^3\theta }{A(\theta ,R)h^2 + B(\theta ,R)\lambda h\sin \theta + C(\theta ,R)\lambda ^2\sin ^2\theta }, \end{equation}

where

(3.17) \begin{equation} \left . \begin{array}{ll} A(\theta ,R) = Rf_1(\theta )f_2(\pi -\theta )+ f_1(\pi -\theta )f_2(\theta ),\\[5pt] B(\theta ,R) = 2\sin ^2\theta [Rf_1(\theta ) + f_1(\pi -\theta )] + 2( R + 1)f_2(\theta )f_2(\pi -\theta ),\\[5pt] C(\theta ,R) = 4\sin ^2\theta [Rf_2(\theta ) + f_2(\pi - \theta ) ],\\[5pt] D(\theta ,R) = 2\sin \theta [R^2f_1(\theta ) + 2Rf_3(\theta ) + f_1(\pi - \theta )],\\[5pt] E(\theta ,R) = 4\sin \theta [R^2f_2(\theta ) + R\pi + f_2(\pi - \theta )]. \end{array}\right \} \end{equation}

The equation is closed by adding the geometry relation

(3.18) \begin{equation} \frac {\mathrm{d}h}{\mathrm{d}s} = \sin \theta . \end{equation}

The interface shape near a moving contact line is thus described by a two-equation system for the unknowns $h(s)$ and $\theta (s)$ . Both equations are expressed by elementary functions and free of any fit parameter. More importantly, the effect of slip has been incorporated. The equations are complemented by the boundary conditions

(3.19) \begin{equation} h(0) = 0, \quad \theta (0)=\theta _e, \quad \kappa (\infty ) = 0. \end{equation}

Several special forms of the GL equation can be obtained. First, setting $\lambda =0$ for the no-slip case, we have

(3.20) \begin{equation} \frac {\mathrm{d}^2\theta }{\mathrm{d}s^2} = -\frac {2Ca}{h^2}\frac {[R^2f_1(\theta ) + 2Rf_3(\theta ) + f_1(\pi - \theta )]\sin ^3\theta }{Rf_1(\theta )f_2(\pi -\theta )+ f_1(\pi -\theta )f_2(\theta )}, \end{equation}

which is the same equation proposed by Chan et al. (Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013) after correcting a sign error therein. Second, when $R=0$ , (3.16) can be simplified to

(3.21) \begin{equation} \frac {\mathrm{d}^2\theta }{\mathrm{d}s^2} = - \frac {CaF(\theta )}{h[h+F(\theta )\lambda ]},\quad \mbox{with}\quad F(\theta ) = \frac {2\sin ^3\theta }{\theta - \sin \theta \cos \theta }. \end{equation}

This equation has the same structure of the GL equation as proposed by Chan et al. (Reference Chan, Srivastava, Marchand, Andreotti, Biferale, Toschi and Snoeijer2013, Reference Chan, Kamal, Snoeijer, Sprittles and Eggers2020), but with the fit parameter $c$ replaced by $F(\theta )$ . Moreover, in the limit of small slope angles ( $\theta \rightarrow 0$ ), we have ${\mathrm{d}x}/{\mathrm{d}s} = \cos \theta \to 1$ , ${\mathrm{d}h}/{\mathrm{d}x}\to \theta$ and $F(\theta )\to 3$ ; accordingly, (3.21) reduces to

(3.22) \begin{equation} \frac {\mathrm{d}^3h}{\mathrm{d}x^3} = - \frac {3Ca}{h(h + 3\lambda )}, \end{equation}

which is the conventional lubrication equation with a Navier-slip condition (Hocking Reference Hocking1983).

4. Asymptotic solutions

We derive an explicit expression of $Q$ by an asymptotic analysis of the GL equation (3.16) under the condition $Ca\ll 1$ . Since the macroscopic region is not considered, the asymptotic structure consists of two regions depending on the importance of the slip effect, i.e. the inner region of size $O(\lambda )$ and the intermediate region with $h\gg \lambda$ .

For small $Ca$ , the inner solution can be expanded as

(4.1) \begin{equation} h(s)=s\sin \theta _e + Cah_1(s) + O({\textit{Ca}}^2),\quad \theta (s)=\theta _e + Ca\theta _1(s) + O({\textit{Ca}}^2). \end{equation}

The leading order simply represents a straight wedge geometry with contact angle $\theta _e$ . Linearising (3.16) around this static solution, we have

(4.2) \begin{equation} \frac {\mathrm{d}^2\theta _1}{\mathrm{d}s^2} = -\frac {Ds+E\lambda }{s(As^2+B\lambda s+C\lambda ^2)}. \end{equation}

The solution is derived as

(4.3) \begin{align} \theta _1&= \frac {s}{2C\lambda }\left (E\ln \frac {As^2+B\lambda s+C\lambda ^2}{As^2} + \frac {2\textit{CD} - \textit{BE}}{\varOmega }\ln \frac {2\textit{As} + B\lambda + \varOmega \lambda }{2As + B\lambda -\varOmega \lambda } \right ) \nonumber\\ &\quad + \frac {D}{2A}\ln \frac {As^2 + B\lambda s + C\lambda ^2}{C\lambda ^2} + \frac {BD-2AE}{2A\varOmega }\left (\ln \frac {2As + B\lambda + \varOmega \lambda }{2As+B\lambda - \varOmega \lambda } - \ln \frac {B + \varOmega }{B - \varOmega } \right ), \end{align}

where we have used the boundary condition $\theta _1(0)=0$ and suppressed the terms yielding finite curvature at large $s$ . The coefficients $A\sim E$ are given by (3.17) with $\theta =\theta _e$ , and $\varOmega =\sqrt {B^2- 4AC}$ . It can be verified that $B^2- 4AC\geqslant 0$ and hence $\varOmega$ is real. The asymptotic behaviour at $s\gg \lambda$ is obtained as

(4.4) \begin{equation} \theta -\theta _e\sim \frac {D}{A}\left (\ln \frac {\mathrm{e}s}{\lambda }+\frac {1}{2}\ln \frac {A}{C} + \frac {2\textit{AE} - \textit{BD}}{2D\varOmega }\ln \frac {B + \varOmega }{B - \varOmega } \right )\boldsymbol{Ca}. \end{equation}

Clearly, the inner solution is not expected to hold as $s\to \infty$ since $\theta _1$ diverges.

In the intermediate region ( $h\gg \lambda$ ), the slip effect can be neglected. The governing equation has been given by (3.20) and is rewritten as

(4.5) \begin{equation} \frac {A(\theta ,R)}{D(\theta ,R)}\frac {\mathrm{d}^2\theta }{\mathrm{d}s^2} = - \frac {\textit{Ca}}{h^2}\sin ^2\theta . \end{equation}

Following Snoeijer (Reference Snoeijer2006), we introduce $h\approx s \sin \theta$ for small $Ca$ . Substituting into the above equation and integrating from $s$ to $\infty$ , we have

(4.6) \begin{equation} \frac {A(\theta ,R)}{D(\theta ,R)}\frac {\mathrm{d}\theta }{\mathrm{d}s}-\int _s^{\infty }\left (\frac {\mathrm{d}\theta }{\mathrm{d}s}\right )^2\frac {\mathrm{d}}{\mathrm{d}\theta }\left [\frac {A(\theta ,R)}{D(\theta ,R)}\right ]\mathrm{d}s= \frac {\textit{Ca}}{s}. \end{equation}

We anticipate $\mathrm{d}\theta /\mathrm{d}s\sim O(\textit{Ca}/s)$ , thereby the integral is of $O(Ca^2/s)$ and can be neglected. Integrating (4.6) once again, we obtain the Cox–Voinov relation

(4.7) \begin{equation} g(\theta ,R) = g(\theta _e,R)+\textit{Ca}\left (\ln \frac {s}{\lambda } + Q\right ), \end{equation}

where

(4.8) \begin{equation} g(\theta ,R) = \int _0^\theta \frac {A(\varphi ,R)}{D(\varphi ,R)} \mathrm{d}\varphi , \end{equation}

is the function first proposed by Cox (Reference Cox1986), and $Q$ is the remaining integral constant.

The constant $Q$ is determined by matching the inner and intermediate solutions. For small deviations of $\theta$ from $\theta _e$

(4.9) \begin{equation} \theta -\theta _e=\frac {g(\theta,R)-g(\theta _e,R)}{\partial_{\theta} g(\theta_e,R)}=\frac {D}{A}\left (\ln \frac {s}{\lambda }+Q\right )\textit{Ca}. \end{equation}

A comparison of (4.4) and (4.9) yields the central results of the present work

(4.10) \begin{equation} Q(\theta _e,R) = \frac {1}{2}\ln \frac {A}{C} + \frac {2\textit{AE} - \textit{BD}}{2D\varOmega }\ln \frac {B+\varOmega }{B-\varOmega } + 1. \end{equation}

In particular, in the presence of a single phase only ( $R=0$ )

(4.11) \begin{equation} Q(\theta _e,0) = \ln \frac {\theta _e-\sin \theta _e\cos \theta _e}{2\sin ^2\theta _e} + 1, \end{equation}

which can be alternatively obtained from a similar analysis of the single-phase equation (3.21).

For comparison, the procedure of matched asymptotic expansions can be also applied to the GL equation (1.2). As a result, we obtain

(4.12) \begin{equation} Q(\theta _e,R)=\ln \frac {\sin \theta _e}{c(\theta _e,R)}+1. \end{equation}

The value of $c$ can be obtained explicitly only under special conditions (Chan et al. Reference Chan, Kamal, Snoeijer, Sprittles and Eggers2020). When $R=0$ , we have

(4.13) \begin{equation} Q = \frac {\theta _e-\sin \theta _e\cos \theta _e}{2\sin ^2\theta _e}h_1(\theta _e)+1, \end{equation}

which has also been reported by Hocking & Rivers (Reference Hocking and Rivers1982). The function $h_1(\theta _e)$ still needs to be calculated by solving numerically an integral equation. Furthermore, an explicit form for $R\neq 0$ is available for $\theta _e=90^\circ$ and reads

(4.14) \begin{equation} Q=\frac {\pi (\pi ^2-4)\left [(R-1)^2h_a+4Rh_b\right ]}{4\pi ^2(R+1)^2-16(R-1)^2}+1, \end{equation}

where $h_a=4(\gamma _E-\ln 2)/\pi$ and $h_b=-1.539$ with $\gamma _E$ denoting the Euler constant. (Note that the coefficient of $h_b$ should be $2$ rather than $2M$ in (2.15b) of Chan et al. (Reference Chan, Kamal, Snoeijer, Sprittles and Eggers2020), which was attribute to a typo in (4.11b) of Hocking (Reference Hocking1977)).

Figure 5.

Variation of $Q$ as a function of $\theta _e$ for $R=0$ . The solid line denotes the present theory (4.11); the circles represent (4.13) with the values listed in table 1 of Hocking & Rivers (Reference Hocking and Rivers1982). The directory including the data and the Jupyter notebook that generated this figure can be accessed at https://multimediafluids.uam.elogim.com/S0022112025105879/JFM-Notebooks/files/figure_5/figure5.ipynb.

Figure 6.

Variation of $Q$ as a function of $R$ for representative values of $\theta _e$ . The solid lines denote the present theory (4.10); the dashed line represents (4.14) valid for $\theta _e=90^\circ$ . The directory including the data and the Jupyter notebook that generated this figure can be accessed at https://multimediafluids.uam.elogim.com/S0022112025105879/JFM-Notebooks/files/figure_6/figure6.ipynb.

The performance of the present theory is first validated in the presence a single phase ( $R=0$ ). In this case, the values of $Q$ have been tabulated in Hocking & Rivers (Reference Hocking and Rivers1982), who calculated $Q$ by solving numerically an integral equation. As shown in figure 5, the present theory given by (4.11) agrees well with Hocking & Rivers (Reference Hocking and Rivers1982) when $\theta _e\lesssim 120^\circ$ . For larger contact angles, pronounced deviations can be observed. The reason is that the flow near a moving contact line becomes poorly approximated by the slippery wedge flow when the opening angle of the wedge is large. We note that this situation is of less interest since in practice the contact angle is usually less than $120^\circ$ for a smooth substrate.

Figure 6 shows the behaviour of $Q$ given by (4.10) over a wide range of $R$ and $\theta _e$ . It is seen that $Q(\theta _e,R)=Q(\pi - \theta _e,R^{-1})$ , which reflects obviously the mirror symmetry of the set-up. A further comparison can be made at $\theta _e=90^\circ$ by plotting (4.14) as the dashed line. It is seen that the present theory exhibits an overall underestimation of $Q$ for $\theta _e=90^\circ$ . Nevertheless, the agreement is still satisfactory considering the approximations used in the derivation. It can be shown below that the error of $Q$ does not produce significant deviations from the numerical results of the Stokes equation.

Figure 7.

Comparison of the present theory with the BEM results at $Ca=0.01$ : (a–c) $R=0$ ; (d–f) $R=1$ . The contact angles are (a, d) $\theta _e=30^{\circ }$ , (b, e) $\theta _e=60^{\circ }$ and (c, f) $\theta _e=90^{\circ }$ . The (black) solid curves represent the full numerical results of the BEM. The (blue) dashed curves represent the numerical solutions of the present GL equation (3.16). The (red) long-dashed curves represent the asymptotic solution (4.7) with $Q$ given by (4.10). The (green) dot-dashed curves in (a–c) correspond to (4.13), first reported by Hocking & Rivers (Reference Hocking and Rivers1982) and available for $R=0$ only. The (green) dot-dashed curve in (f) corresponds to (4.14), which is based on Chan et al. (Reference Chan, Kamal, Snoeijer, Sprittles and Eggers2020) and valid for $\theta _e=90^{\circ }$ . The directory including the data and the Jupyter notebook that generated this figure can be accessed at https://multimediafluids.uam.elogim.com/S0022112025105879/JFM-Notebooks/files/figure_7/figure7.ipynb.

5. Comparison with numerical results

To validate the present GL equation and asymptotic solutions, we show in figure 7 a comparison between the solutions of the GL equation (3.16) and the Stokes equation by plotting the variation of $[g(\theta ,R)-g(\theta _e,R)]/Ca$ with $s/\lambda$ . The results are illustrated for $Ca=0.01$ together with two viscosity ratios $R=0$ and 1 and three contact angles $\theta _e=30^\circ$ , $60^\circ$ and $90^\circ$ ; the employment of different yet small values of $Ca$ produced indistinguishable curves. The numerical solutions of the Stokes equations are obtained by the BEM and represented by the solid curves. The numerical solutions of the GL equation are calculated using a shooting method. The two sets of numerical results exhibits good agreement, especially for $R=1$ .

The Cox–Voinov law is represented by the straight lines, where the values of $Q$ are given by (4.10) for the present theory, (4.13) for $R=0$ (Hocking & Rivers Reference Hocking and Rivers1982) and (4.14) for $R=1$ and $\theta _e=90^\circ$ (Chan et al. Reference Chan, Kamal, Snoeijer, Sprittles and Eggers2020). All of the theoretical curves agree reasonably well with the numerical solutions. A close examination of the results show slight deviations in figure 7(c, d). For $R=0$ and $\theta _e=90^\circ$ in figure 7(c), the curve of the present asymptotic theory lies below that of Hocking & Rivers (Reference Hocking and Rivers1982), which is due to the underestimation of $Q$ shown in figure 6. This deviation is reduced at $R=1$ , leading to an improved comparison in figure 7(f). The slight deviation in figure 7(d) is associated with the small contact angle and can be improved at smaller $Ca$ .