2.1. Drops with moving contact lines

Both the precursor-film and Navier-slip models are commonly used to regularise the viscous stress singularity at moving contact lines (Huh & Scriven Reference Huh and Scriven1971; Sibley et al. Reference Sibley, Nold and Kalliadasis2015). The precursor-film model is particularly suitable for problems involving free-surface topological changes such as breakup or rupture (Mhatre & Kumar Reference Mhatre and Kumar2024), whereas the Navier-slip model is more straightforward and conceptually simpler. For the present two-dimensional drop deformation and motion problem with moving contact lines, we therefore adopt the Navier-slip model (see, Hocking Reference Hocking1983; Eggers Reference Eggers2005; Snoeijer & Andreotti Reference Snoeijer and Andreotti2013; Qin, Lu & Gao Reference Qin, Lu and Gao2024). The slip length is denoted by $\lambda$ , which is typically a few nanometres. Then, following the procedure outlined by Oron et al. (Reference Oron, Davis and Bankoff1997), we derive the pressure and velocity distributions within the drop as follows:

(2.3a) \begin{align} p&=\rho g\!\left (h-z\right )-\sigma \frac {\partial ^2 h}{\partial x^2} +p_0, \\[-12pt]\nonumber \end{align}

(2.3b) \begin{align} u&=\frac {1}{2\mu }\frac {\partial p}{\partial x}\big (z^2-2hz-2\lambda h\big )+\frac {\tau }{\mu }(z+\lambda ), \end{align}

where $g$ is the gravitational acceleration.

Accordingly, the lubrication equation governing the deformation and motion of thin drops reads

(2.4) \begin{align} \frac {\partial h}{\partial t} +\frac {1}{3\mu }\frac {\partial }{\partial x} \left [h^2(h+3 \lambda )\left (\sigma \frac {\partial ^3 h}{\partial x^3}+\rho g \frac {\partial h}{\partial x}\right )\right ] +\frac {\tau }{2\mu }\frac {\partial }{\partial x}\left [h(h+2\lambda )\right ]=0. \end{align}

Through the scaled variables

(2.5) \begin{align} x=RX,\quad z=R\theta _eZ,\quad h=R\theta _eH,\quad \lambda =R\theta _e \varLambda ,\quad t=\frac {3\mu R}{\sigma \theta _e^3} T, \end{align}

we non-dimensionalise the governing equation (2.4) to

(2.6) \begin{align} \frac {\partial H}{\partial T} +\frac {\partial }{\partial X} \left [H^2(H+3 \varLambda )\left (\frac {\partial ^3 H}{\partial X^3}+Bo\frac {\partial H}{\partial X}\right )\right ] +S\frac {\partial }{\partial X}\left [H(H+2\varLambda )\right ]=0, \end{align}

where

(2.7) \begin{align} S=\frac {3\tau R}{2\sigma \theta _e^2} \end{align}

is the dimensionless parameter characterising the ratio of shear stress to the surface tension, and the Bond number $Bo=\rho gR^2/\sigma$ characterises the relative importance of gravity to surface tension. We focus on small drops such that the gravitational effect is negligible, i.e. $Bo=0$ .

After non-dimensionalisation, the initial parabolic drop profile reads

(2.8) \begin{align} H_{\textit{init}}=X-\frac {X^2}{2}, \end{align}

indicating that the contact lines are initially located at $X=0$ and 2.

We assume that the advancing contact angle equals its receding counterpart, i.e. the contact-angle hysteresis is neglected. The fourth-order equation (2.6) is supplemented by the following conditions:

(2.9) \begin{align} H=0,\quad \frac {\partial H}{\partial X}=1 \quad {\mathrm{at}} \quad X=L_{\textit{rec}}(T), \\[-12pt]\nonumber \end{align}

(2.10) \begin{align} H=0,\quad \frac {\partial H}{\partial X}=-1 \quad {\mathrm{at}} \quad X=L_{\textit{ad}v}(T), \end{align}

where $L_{\textit{rec}}(T)$ and $L_{\textit{ad}v}(T)$ denote the positions of the receding and advancing contact lines, respectively. As the contact lines are moving, additional boundary conditions are required

(2.11) \begin{align} c_{\textit{rec}}=3\varLambda \lim _{x\to L_{\textit{rec}}^+}H\frac {\partial ^3 H}{\partial X^3}+2\varLambda S ,\,\,c_{\textit{ad}v}=3\varLambda \lim _{x\to L_{\textit{ad}v}^-}H\frac {\partial ^3 H}{\partial X^3}+2\varLambda S , \end{align}

where $c_{\textit{adv/rec}} \equiv {\mathrm{d}} L_{\textit{adv/rec}}/{\mathrm{d}} t$ represents the advancing/receding contact-line velocity. These conditions guarantee that the contact lines are material lines and can be obtained by calculating the limit of the depth-averaged velocity towards the contact lines.

The governing equation (2.6) together with the initial and boundary conditions (2.8)–(2.11) is solved numerically, with details provided in Appendix A.

2.2. Drops with pinned contact lines

Due to chemical heterogeneity, surface roughness or deliberate experimental constraints, a drop may become pinned (see, Savva et al. Reference Savva, Groves and Kalliadasis2019; Lohse Reference Lohse2022; Mhatre & Kumar Reference Mhatre and Kumar2024). In this case, the viscous stress (associated with $H H_{\textit{XXX}}$ ) is not singular in the absence of slip, as inferred from the local asymptotic expansion of the film thickness in Appendix B. Therefore, the slip models are not implemented. Accordingly, the non-dimensional governing equation (2.6) becomes

(2.12) \begin{align} \frac {\partial H}{\partial T} +\frac {\partial }{\partial X} \left (H^3\frac {\partial ^3 H}{\partial X^3}\right ) +S\frac {\partial }{\partial X}\big(H^2\big)=0. \end{align}

The initial drop profile is the same as (2.8). At the contact line, we expect the flux to vanish (Bernis et al. Reference Bernis, Hulshof and King2000; Eggers & Fontelos Reference Eggers and Fontelos2025), yielding the following boundary conditions:

(2.13) \begin{align} H(0,T)=H(2,T)=0,\quad \left (H^3\frac {\partial ^3 H}{\partial X^3}+\textit{SH}^2\right )\Bigg |_{X=0,2}=0. \end{align}

The governing equation (2.12) together with boundary conditions (2.13) is solved numerically, with details provided in Appendix B.

3.1. Low shear: $S\leq S_{\textit{cr}}$

Using the numerical algorithm described in Appendix B.1, we find that, under low shear, a drop evolves from its initial parabolic shape to a steady state. Focusing on the steady state of the drops with the profile $H(X)$ , the governing equation (2.12) can be integrated once to give

(3.2) \begin{align} \textit{HH}'''+S=0, \end{align}

where the conditions of vanishing flux have been used. The steady solution can be determined by the conditions of vanishing thickness at the contact lines and mass conservation

(3.3) \begin{align} H(0)=H(2)=0, \quad \int _0^2 H\,\mathrm{d}X=\frac 23. \end{align}

The steady-state problem (3.2) and (3.3) is solved using the numerical algorithm described in Appendix B.2. The solutions are computed for a sequence of shear numbers $S$ with a uniform increment $\Delta S=0.001$ , from which the steady drop profiles and the corresponding microscopic contact angles are obtained. As shown in figure 2, the steady-state drop profiles lean to the downstream side with increasing shear rate. We denote the upstream-side and downstream-side microscopic contact angles as $\theta _{\textit{up}}$ and $\theta _{{down}}$ , respectively; $\theta _{\textit{up}}$ decreases with increasing shear number $S$ and vanishes at $S=S_{\textit{cr}}$ and $\theta _{{down}}$ increases with $S$ , reaching a maximum value of 1.534 at $S=S_{\textit{cr}}$ .

Figure 2.

Steady-state pinned drops under low shear obtained from numerical solutions: (a) drop profiles; (b) variation of the upstream-side microscopic contact angle $\theta _{\textit{up}}$ and downstream-side microscopic contact angle $\theta _{{down}}$ with respect to the shear number $S$ .

Figure 3.

Force balance analysis of a steady drop with pinned contact lines under low shear.

To obtain the quantitative relationship among $\theta _{\textit{up}}$ , $\theta _{{down}}$ and $S$ , we perform a force analysis of a steady-state drop, as sketched in figure 3. Considering only the horizontal force balance, the drop is subject to shear stress $F_\tau$ , viscous force $F_\mu$ and capillary forces acting tangentially on the surface of the liquid film at the contact lines. Under the lubrication approximation, the total shear force is $F_\tau =2\tau R$ . Based on velocity distribution (2.3b ) (with slip length $\lambda =0$ ), we have

(3.4) \begin{align} \frac {\partial u}{\partial z}\Bigg |_{z=0}=\frac {\sigma }{\mu }hh'''+\frac {\tau }{\mu }=-\frac {\tau }{2\mu }. \end{align}

Here, we have used the dimensional steady-state governing equation (3.2), which reads $h h'''+3\tau /2\sigma =0$ . Thus the total viscous shear force induced by the substrate is

(3.5) \begin{align} F_\mu =-2\mu R \frac {\partial u}{\partial z}\Bigg |_{z=0}=\tau R. \end{align}

According to the horizontal force balance on the drop, we then arrive at

(3.6) \begin{align} \sigma \cos (\theta _{\textit{up}}\theta _e)-\sigma \cos (\theta _{{down}}\theta _e)=F_\tau +F_\mu =3\tau R. \end{align}

Equation (3.6) is consistent with the force analysis result by Dimitrakopoulos & Higdon (Reference Dimitrakopoulos and Higdon1997), who focused on the depinning condition of the drop under shear. Expanding (3.6) with respect to $\theta _e$ yields

(3.7) \begin{align} \theta _{{\textit{down}}}^2-\theta _{{\textit{up}}}^2-4S=0. \end{align}

The correctness of (3.7) is directly verified by the numerical results in figure 2, which shows that $\theta _{{down}}^2-\theta _{\textit{up}}^2$ grows linearly with $S$ . At the threshold, we have $\theta _{\textit{up}}=0$ and hence the maximum downstream angle

(3.8) \begin{align} \theta _{{down}}=2\sqrt {S_{\textit{cr}}}. \end{align}

As $\theta _{\textit{up}}$ vanishes at the threshold, the upstream-side asymptotic expansion for problem (3.2) and (3.3) has the following form:

(3.9) \begin{align} H \sim \sqrt {\frac {8S}{3}}X^{\frac 32}+P_1X^{\frac {5+\sqrt {13}}{4}}\quad \mbox{as}\quad X\to 0^+, \end{align}

where $P_1$ is a free parameter to be determined in the numerical process. Interestingly, (3.9) is formally similar to the near-contact-line behaviour of a drop on a completely wetting substrate (Benilov & Benilov Reference Benilov and Benilov2015). Moreover, asymptotic expansion (3.9) loses $\theta _{\textit{up}}$ relative to (B3), because the upstream contact angle vanishes at the critical case. To restore a closed boundary-value problem, the critical shear number $S_{\textit{cr}}$ is treated as an additional unknown parameter. The resulting system is then solved simultaneously using the numerical algorithm described in Appendix B.2, yielding a precise value of the critical shear number $S_{\textit{cr}} = 0.58857$ . The corresponding drop profile at $S_{\textit{cr}}$ is shown in figure 2(a).

From another perspective, the volumetric flux for a pinned drop is

(3.10) \begin{align} Q=\int _0^h u {\mathrm{d}} z=\frac {\sigma }{3\mu }h^3h'''+\frac {\tau }{2\mu }h^2. \end{align}

The steady state corresponds to $Q=0$ , implying that the shear-induced internal flow is exactly balanced by the capillary pressure gradient due to interface curvature. The critical shear number $S_{\textit{cr}}$ corresponds to the limit at which the drop can no longer generate a sufficient capillary pressure gradient to balance the shear-driven flux. In the present problem, this occurs when the upstream contact angle $\theta _{\textit{up}}$ approaches zero.

The definition of the shear number $S$ in (2.7) indicates an explicit relationship between the corresponding critical shear rate, drop size, surface tension, wettability and the viscosity of the displacing fluid. The present result shows that the critical shear rate is independent of the viscosity of the drop itself. Although the present model assumes a strictly pinned contact line corresponding to a $180^\circ$ contact-angle hysteresis, experiments on three-dimensional drops with finite hysteresis also show that the critical shear strength depends weakly on the drop viscosity (Fan et al. Reference Fan, Wilson and Kapur2011; Xu et al. Reference Xu, Xia, Huang and Li2024).

3.2. High shear: $S\gt S_{\textit{cr}}$

We take $S = 1 \gt S_{\textit{cr}}$ as an example to explore the evolution of the pinned drop under high shear. As shown in figure 4(a), the drop near the upstream-side contact line gradually flattens, forming a thinning film with vanishing thickness. The main body of the drop evolves towards a quasi-steady shape, which connects to the thin film at an apparent contact line.

Figure 4.

Unsteady evolution of a pinned drop under high shear $S=1$ : (a) drop profiles; (b) in the late stage of evolution ( $T \gg 1$ ), a close-up view of the vanishing film region, with the height magnified by $2ST$ .

Figure 4. Long description

Panel A: A line graph shows drop profiles under varying conditions. The horizontal axis is labeled X, and the vertical axis is labeled H. Different lines represent different time values (T = 1, T = 3, T = 5, T = 10, T = 50), with a legend indicating these values. The graph shows how the height (H) of the drop changes with respect to the horizontal position (X). The apparent contact line is marked on the graph. Panel B: Another line graph provides a close-up view of the vanishing film region in the late stage of evolution. The horizontal axis is labeled X, and the vertical axis is labeled 2STH. Different lines represent different time values (T = 50, T = 100, T = 150, T = 200, T = 250), with a legend indicating these values. The graph shows the relationship between 2STH and X, with a line indicating where 2STH equals X.

It is shown that the profile of the vanishing film becomes nearly straight in the late stage of the evolution, so that capillary force can be neglected in this region. Accordingly, the reduced governing equation and the corresponding asymptotic solution read

(3.11) \begin{align} \begin{aligned} &\frac {\partial H}{\partial T} +S\frac {\partial }{\partial X}\big(H^2\big)=0 \ \Longrightarrow \ H\sim \frac {X}{2ST} \quad \mbox{as}\quad T\to \infty . \end{aligned} \end{align}

The accuracy of asymptotic solution (3.11) is confirmed in figure 4(b). Therefore, the thickness of the thin film continuously decreases but can never reach a steady state.

Figure 5.

Unsteady evolution of a pinned drop under high shear $S=1$ : the variations of $\theta _{{\textit{up}}}$ and $\theta _{{\textit{down}}}$ with respect to the dimensionless time $T$ .

Figure 5 presents the temporal evolution of the microscopic contact angles on both sides of the drop. The upstream-side microscopic contact angle $\theta _{\textit{up}}$ gradually decreases and asymptotically approaches zero, while the downstream-side contact angle $\theta _{{down}}$ increases and eventually reaches 1.898. Figure 5(b) also confirms the validity of asymptotic solution (3.11), according to which we have $\theta _{\textit{up}}=1/(2ST)$ .

In the late stage of the evolution, almost all of the liquid is accumulated into the main body of the drop, which evolves in a quasi-steady manner. It can be shown that the main body of the drop is related to the critical steady solution developed above. The profile of the main body of the drop satisfies

(3.12) \begin{align} \textit{HH}'''+S=0,\quad H(X_{\textit{cl}})=H'(X_{\textit{cl}})=H(2)=0,\quad \int _{X_{\textit{cl}}}^{2}H\,{\mathrm{d}}X=\frac 23, \end{align}

where $X=X_{\textit{cl}}$ denotes the position of the apparent contact line. Here, we have neglected the film thickness at the apparent contact line and imposed the vanishing slope condition since the main drop should conform to the thin film. Introducing the rescaling

(3.13) \begin{align} X=\frac {2-X_{\textit{cl}}}{2}\tilde X+X_{\textit{cl}}, \quad H=\frac {2}{2-X_{\textit{cl}}}\tilde H, \end{align}

equation (3.12) is transformed into

(3.14) \begin{align} \tilde H\!\tilde H'''+\left (\frac {2-X_{\textit{cl}}}{2}\right )^5S=0,\quad \tilde H(0)=\tilde {H}'(0)=\tilde H(2)=0,\quad \int _0^2\tilde H\,{\mathrm{d}}\tilde X=\frac 23, \end{align}

which is formally equivalent to the steady problem (3.2) and (3.3) with $S$ taken as $S_{\textit{cr}}$ . A comparison of the constant terms thus yields

(3.15) \begin{align} X_{\textit{cl}}=2-2\left (\frac {S_{\textit{cr}}}{S}\right )^{1/5}. \end{align}

Furthermore, recalling that the critical downstream-side microscopic contact angle $\theta _{{down}}=2\sqrt {S_{\textit{cr}}}$ and using $H'(X)= (S/S_{\textit{cr}} )^{2/5}\tilde H'(\tilde X)$ , we obtain the late-stage downstream-side microscopic contact angle as

(3.16) \begin{align} \begin{aligned} \theta _{{down}}=-H'(2)=2S_{\textit{cr}}^{1/10}S^{2/5}. \end{aligned} \end{align}

The theoretical predictions (3.15) and (3.16) agree well with the numerical results, as shown in figure 6.

Figure 6.

Late-stage behaviour of drop evolution under high shear: (a) theoretical prediction (3.15) and numerical results for $X_{\textit{cl}}$ ; (b) theoretical prediction (3.16) and numerical results for the late-stage downstream-side microscopic contact angle $\theta _{{down}}$ .

As $T\to +\infty$ , the local film profile near $X=X_{\textit{cl}}$ develops a bending structure connecting the thin film and the main body of the drop. We show here that this local film profile has a vanishing size as $T\to \infty$ and is governed by a similarity solution, which can be obtained by introducing the ansatz

(3.17) \begin{align} H(X,T)=\frac {X_{\textit{cl}}}{2ST}\mathcal{H}(\xi ), \quad \xi =S\left (\frac {2T}{X_{\textit{cl}}}\right )^{\frac 23}(X-X_{\textit{cl}}), \end{align}

where $\mathcal{H}(\xi )$ is the similarity profile. Away from the apparent contact line, this structure should match the main body of the drop and the thinning film, which are described locally by $H\approx \sqrt {8S/3}(X-X_{\textit{cl}})^{3/2}$ and $H\approx X_{\textit{cl}}/(2ST)$ according to (3.9) and (3.11), respectively. For the similarity variables, these conditions become

(3.18) \begin{align} \mathcal{H}(\xi ) \sim \begin{cases} 1, & \text{as} \quad \xi \to -\infty , \\[5pt] \sqrt {\dfrac {8}{3}}\xi ^{\frac 32}, & \text{as} \quad \xi \to +\infty . \end{cases} \end{align}

After transformation (3.17), the time derivative term $\partial H/\partial T$ in the governing equation (2.12) is $O(T^{-2})$ and can be neglected compared with the remaining terms that are $O(T^{-4/3})$ as $T\to \infty$ . Accordingly, after one integration, (2.12) reduces to

(3.19) \begin{align} \mathcal{H}^3\mathcal{H}'''+\mathcal{H}^2=1, \end{align}

where the integration constant is determined by the boundary condition as $\xi \to -\infty$ . The boundary-value problem (3.18) and (3.19) can be conveniently solved using a shooting method. Figure 7 shows the similarity profile $\mathcal{H}(\xi )$ and compares it with the numerical results of the full problem at $S=1$ . In the late stage of the evolution, the profile of the capillary transition region gradually approaches the theoretical prediction. It is noteworthy that the calculation becomes increasingly challenging since the curvature near the apparent contact line, given by $\partial ^2 H/\partial X^2=S(2T/X_{\textit{cl}})^{1/3}\mathcal{H}''(\xi )$ , diverges as $T\to \infty$ .

In summary, under the small-contact-angle assumption, when the shear number exceeds the critical value $S\gt S_{\textit{cr}}$ , the pinned state on the upstream side of the drop can be considered effectively lost, while the main body of the drop evolves in a quasi-steady manner. An idealised contact-angle hysteresis of $180^\circ$ is adopted in the present study, thus, the critical shear for ‘depinning’ is linked to the condition that upstream contact angle goes to zero. However, for realistic drops with finite contact-angle hysteresis, depinning can occur earlier, either at the upstream or downstream side, depending on which contact angle first exceeds the hysteresis range. For drops with large contact angles, exceeding the critical shear leads to pinch-off and eventual detachment on the downstream side, which do not occur in the present work (Schleizer & Bonnecaze Reference Schleizer and Bonnecaze1999; Spelt Reference Spelt2006; Ding & Spelt Reference Ding and Spelt2008; Ding et al. Reference Ding, Gilani and Spelt2010).

Figure 7.

Profile of the capillary transition region connecting the vanishing film and main body in the late stage of unsteady drop evolution under high shear: comparison between the theoretical result (black line) and numerical profiles with $S=1$ .

4.1. Low shear

The drops under low shear are characterised by small velocity $c$ , and an asymptotic theory can be developed following the procedure proposed by Hocking (Reference Hocking1983), who considered the spreading of a thin drop in the absence or presence of gravity. For the present problem, the profile of the steady-state drop consists of an outer region and two inner regions near the advancing/receding contact lines.

In the outer region, the slip effect is insignificant, and thus the governing equation (4.1) reduces to

(4.4) \begin{align} H^2H''' +\textit{SH}=c. \end{align}

Considering that both the shear number $S$ and the velocity $c$ are small, we expand $H$ to first order as

(4.5) \begin{align} \begin{aligned} H=H_0 + cH_1 + S H_2 + O(c^2, S^2, cS), \end{aligned} \end{align}

where $H_0$ , $H_1$ and $H_2$ represent the static drop profile and the first-order corrections due to flow and shear, respectively. Substituting (4.5) into (4.4), we obtain

(4.6) \begin{align} \begin{aligned} &H_0''' = 0, \quad H_0(0)=H_0(L)=0, \quad \int _0^L H_0 \, {\mathrm{d}}X = V,\\[3pt] &H_1''' =\frac {1}{H_0^2},\ H_1(0)=H_1(L)=0,\,\int _0^LH_1{\mathrm{d}}X=0,\\[3pt] &H_2''' =-\frac {1}{H_0},\ H_2(0)=H_2(L)=0,\,\int _0^LH_2{\mathrm{d}}X=0. \end{aligned} \end{align}

Solutions to (4.6) are

(4.7) \begin{align} \begin{aligned} &H_0=\frac {6V}{L^3}X(L-X),\\[3pt] &H_1=\frac {L^3(L-X)X}{36V^2}{\mathrm{ln}}\frac {L-X}{X},\,\\[3pt] &H_2=\frac {L^2}{12V}\left [(L-X)^2{\mathrm{ln}}\frac {L-X}{L}-X^2{\mathrm{ln}}\frac {X}{L}\right ]\!. \end{aligned} \end{align}

It can be shown that the slope of the outer solution exhibits logarithmic singularities near the moving contact lines. Following Eggers (Reference Eggers2005), for the purpose of matching, we provide the cube of the slope near the moving contact lines as

(4.8) \begin{align} \left (H_{\textit{out}}'\right )^3 &\to \begin{cases} \left (\dfrac {6V}{L^2}\right )^3 - 9 S \dfrac {V}{L} + 3c\ln \dfrac {L}{\mathrm{e}X}, & \text{for}\quad 0\lt X\ll 1, \\[11pt] -\left (\dfrac {6V}{L^2}\right )^3 - 9 S \dfrac {V}{L} + 3c\ln \dfrac {L}{\mathrm{e}(L-X)}, & \text{for}\quad 0\lt L-x\ll 1. \end{cases} \end{align}

The inner region is characterised by a length of $\varLambda$ , over which a significant variation of the surface slope occurs due to the large viscous force therein. The analysis of (2.6) in the inner region reveals that the shear stress can be neglected near the contact line. Therefore, as shown in many previous studies (see, Hocking Reference Hocking1983; Eggers Reference Eggers2005; Benilov & Benilov Reference Benilov and Benilov2015; Sibley et al. Reference Sibley, Nold and Kalliadasis2015; Savva, Rezgdanikov & Colinet Reference Savva, Rezgdanikov and Colinet2017; Wang, Qin & Gao Reference Wang, Qin and Gao2025), the inner solutions exhibits the following asymptotic behaviour in the far field:

(4.9) \begin{align} \left (H_{\textit{in}}'\right )^3 &\to \begin{cases} 1 - 3c\, \ln \!\dfrac {\mathrm{e} X}{3\varLambda }, & \text{for}\quad \varLambda \ll X\ll 1,\\[11pt] -1 - 3c \ln \!\dfrac {\mathrm{e} (L - X)}{3\varLambda }, & \text{for}\quad \varLambda \ll L-X\ll 1. \end{cases} \end{align}

A comparison of (4.8) and (4.9) shows that they can indeed match, yielding

(4.10) \begin{align} 3c\ln \frac {L}{3\varLambda } = 9 S \frac {V}{L} - \left (\frac {6V}{L^2}\right )^3 + 1, \\[-12pt]\nonumber \end{align}

(4.11) \begin{align} 3c\ln \frac {L}{3\varLambda } = 9 S \frac {V}{L} + \left (\frac {6V}{L^2}\right )^3 - 1 . \\[10pt] \nonumber \end{align}

Given $V=2/3$ , we finally obtain

(4.12) \begin{align} L=2, \quad c=\frac {S}{\ln \dfrac {2}{3\varLambda }}. \end{align}

Figure 10.

Steady-state drops under low shear: (a) numerical profiles (solid lines) and theoretical profiles (4.5) (dashed lines) for $S=0.1,\,0.3,\,0.5$ and $\varLambda =10^{-5}$ with the receding contact line shifted to $X=0$ , where the $L$ in (4.5) comes from the numerical solutions; (b) comparison between the numerical velocities $c$ (solid lines) and asymptotic predictions (4.12) (dashed lines) for slip lengths $\varLambda =10^{-4},\,10^{-5},10^{-6}$ .

As exhibited in figure 10, the analytical results for the steady-state drop velocity $c$ and drop profile $H(X)$ show good agreements with the numerical results in the low-shear regime. When $S$ is small, the steady drop length $L$ remains close to 2, and $c$ is approximately proportional to $S$ , with a proportionality coefficient given by $1/\ln (2/3\varLambda )$ , as shown in figure 10(b). Ding et al. (Reference Ding, Gilani and Spelt2010) studied steadily moving drops in shear flow and reported that the drop speed increases linearly with shear rate in the low-Reynolds-number regime, consistent with the present predictions.

Asymptotic results (4.12) indicate that the drop essentially retains its original length in the low-shear regime. Although this prediction shows slight deviations from the numerical results in figure 8(b), the agreement remains reasonably good. Since the direction of shear (i.e. the sign of $S$ ) does not affect the drop length because of symmetry, a more refined theoretical description of $L$ would require including higher-order even powers of $S$ . Deriving these higher-order corrections involves considerably more algebraic complexity, and is beyond the scope of the present work.

Moreover, in the low-shear regime where the drop deformation is small and the drop profile remains close to the parabolic shape, the apparent contact angle $\theta _{\textit{ap}}$ can be approximated by the interface slope of the leading-order outer solution evaluated at the contact line, i.e. $6V/L^2$ . Under this approximation, the matching results (4.10) and (4.11) can be extended to the unsteady drop evolution and lead to the following Cox–Voinov-like law:

(4.13) \begin{align} \theta _{\textit{ap,rec}}^3 -1= -3c\,{\mathrm{ln}}\frac {L}{3\lambda } + 9S\frac {V}{L},\,\theta _{\textit{ap},adv}^3 -1 = 3c\,{\mathrm{ln}}\frac {L}{3\lambda } - 9S\frac {V}{L}. \end{align}

Equation (4.13) retains the characteristic cubic dependence of the Cox–Voinov law, with an effective macroscopic scale $L/3$ and microscopic scale $\lambda$ , while the additional correction term $9SV/L$ reflects the influence of the external shear on the outer-region drop profile.

4.2. Intermediate shear

Figure 11 shows that the drops under intermediate shear can no longer maintain the parabolic-like shape observed in the low-shear region. Instead, a flat film rapidly forms in the central region of the drop, and the drop length $L$ experiences a sharp increase. Moreover, the non-dimensional time required for the drop to reach steady state increases by an order of magnitude.

Figure 11.

Steady-state drop profile under intermediate shear with $\varLambda =10^{-5}$ and representative values of $S$ . The receding contact line is shifted to $X=0$ .

As shown in figure 8(a), the steady-state velocity $c$ of the drop exhibits a turning point at $S=0.535$ (for $\varLambda =10^{-5}$ ), where $c$ begins to temporarily decrease despite increasing $S$ . Comparing this with figure 9, we observe that, for $S = 0.535 \lt 0.545$ , the shear driving effect, capillary resistance effect and viscous resistance effect still grow linearly with $S$ , i.e. $cL \propto S.$ Since the drop length $L$ increases steeply with $S$ in this regime, the increase in drop motion cost is slower than that of $L$ , leading to a temporary decline in the drop speed  $c$ . Direct experimental measurements of drops under well-controlled constant shear are scarce. Nevertheless, qualitatively similar non-monotonic behaviour of the drop velocity has been reported in other configurations, i.e. gravity-driven drops on an inclined plane (Benilov & Benilov Reference Benilov and Benilov2015) and drop motion on strained soft substrates (Chao, Jeon & Karpitschka Reference Chao, Jeon and Karpitschka2025).

4.3. High shear: steady-state profile

In the high-shear regime, figure 12 shows the steady-state drop profiles at several representative shear numbers, illustrating the emergence of a pronounced flat film region in the drop rear. In the present study, the thin film effectively decouples the steady drop into three distinct zones: entrained film, dimple and capillary ridge. Similar elongation and structural features, including a thin trailing film, capillary ridge and the connecting dimple, have been observed qualitatively in experiments with moving drops (Fan et al. Reference Fan, Wilson and Kapur2011; Wang et al. Reference Wang, Chang, Zhao and Yang2022; Xu et al. Reference Xu, Xia, Huang and Li2024).

Figure 12.

Steady-state drop profiles under high shear with $\varLambda =10^{-5}$ . The receding contact line is shifted to $X=0$ .

4.3.1. Entrained film

The drop rear is characterised by the entrainment of thin film with a receding contact line. Its profile can also be determined via asymptotic analysis. We divide the drop rear region into an outer region and an inner region. In the outer region, the slip effect is insignificant and thus neglected, the reduced governing equation is again (4.4). The corresponding boundary conditions read

(4.14) \begin{align} H(0) \to 0, \quad H(+\infty ) \to \frac {c}{S}, \end{align}

where the far-field film thickness can be easily obtained by neglecting the capillary force term in (4.4). Using the coordinate rescalings

(4.15) \begin{align} H = \frac {c}{S} \bar H, \quad X = \frac {c^{2/3}}{S} \bar X, \end{align}

we transform problem (4.4) and (4.14) into

(4.16) \begin{align} \begin{aligned} &\bar H^2 \bar H''' + \bar H = 1, \quad \bar H(0) \to 0, \quad \bar H(+\infty ) \to 1. \end{aligned} \end{align}

When $\bar X \to 0^+$ , $\bar {H}$ has the following asymptotics:

(4.17) \begin{align} \bar H \sim \bar X \left [-3 \ln (k_0 \bar X)\right ]^{1/3}. \end{align}

We obtain $k_0 = 2.19$ by solving (4.16) with a shooting method and fitting the solution in the interval $\bar X \in [10^{-6}, 10^{-4}]$ to (4.17). Accordingly, we have

(4.18) \begin{align} \left (H_{\textit{out}}'\right )^3 \sim -3c \ln \frac {\mathrm{e} S k_0 X}{c^{2/3}}. \end{align}

Matching (4.18) with the inner asymptotic behaviour in (4.9) yields the relation

(4.19) \begin{align} c^{\frac 23} e^{-\frac {1}{3c}} = 3k_0 \varLambda S \quad \Rightarrow \quad c = \frac {1}{2 \mathscr{W} \left [ \dfrac {1}{2}(3k_0\varLambda S)^{-3/2} \right ]}, \end{align}

where $\mathscr{W}$ denotes the principal branch of the Lambert $\mathscr{W}$ function defined as the solution to $\mathscr{W}(x) e^{\mathscr{W}(x)} = x$ with $x\gt 0$ , and admits the following asymptotic expansion:

(4.20) \begin{align} \mathscr{W}(x) \sim \ln x - \ln \ln x + \frac {\ln \ln x}{\ln x} + {\cdots} , \quad {\mathrm{as}} \quad x \to +\infty . \end{align}

As shown in figure 13, both the drop rear profiles obtained via the shooting method and the steady velocity $c$ predicted by (4.19) agree well with the numerical results.

Figure 13.

(a) Comparison between numerical profiles (solid lines) and shooting results (dashed lines) of the entrained film. (b) Comparison of steady-state drop velocity $c$ as a function of $S$ between numerical results (solid line) and theoretical prediction (4.19) (dashed line). The slip length $\varLambda =10^{-5}$ .

4.3.2. Capillary ridge

For small drop velocity $c$ , the capillary ridge is primarily determined by the balance between the capillary and shear forces. In what follows, we only compute the leading-order approximation of the drop height $H$ in the capillary ridge region, omitting the effect of contact-line movement. For the macroscopic profile of the capillary ridge, the governing equation (4.1) reduces to

(4.21) \begin{align} \textit{HH}'''+S=0, \end{align}

which is the same as (3.2) for pinned drops.

Figure 14.

(a) Comparison between the leading-order prediction (4.30) with $\alpha =20$ and numerical results for the maximum height of the capillary ridge. (b) Drop profile in the capillary ridge region with $\varLambda = 10^{-5}$ : comparison between the leading-order outer solutions and the full numerical results. In both panels, dashed curves represent asymptotic results.

For convenience, we shift the advancing contact line to $X=0$ , so that the right boundary condition reads

(4.22) \begin{align} \begin{aligned} H(0)=0. \end{aligned} \end{align}

Regarding the left-side boundary conditions, Benilov & Benilov (Reference Benilov and Benilov2015) and Wang et al. (Reference Wang, Qin and Gao2025) – in their studies on drop sliding and spin coating of liquid films – assumed that the capillary ridge approaches the $X$ -axis with zero contact angle on the left. We adopt the same assumption here and neglect the thickness of the entrained film, as illustrated by the dashed line in figure 14(b). Thus, the left boundary conditions of (4.21) are prescribed as

(4.23) \begin{align} \begin{aligned} H(-L_{\textit{ridge}})=0,\quad H'(-L_{\textit{ridge}})=0, \end{aligned} \end{align}

where $L_{\textit{ridge}}$ denotes the length of the capillary ridge.

It can be shown that the leading-order solution of the capillary ridge resembles the profile of a pinned drop under critical shear $S_{\textit{cr}}$ , as discussed in § 3.2. Supposing the leading-order capillary ridge has a volume $V_{\textit{ridge}}$ , we introduce

(4.24) \begin{align} X=\frac {L_{\textit{ridge}}}{2}(\tilde X-2), \quad H=\frac {3V_{\textit{ridge}}}{L_{\textit{ridge}}} \tilde H, \end{align}

to transform problem (4.21)–(4.23) into

(4.25) \begin{align} \tilde H\!\tilde H'''+\frac {L_{\textit{ridge}}^5S}{72V_{\textit{ridge}}^2}=0,\,\tilde H(0)=\tilde H(2)=0,\,\tilde H'(0)=0,\,\int _0^2\tilde H\,{\mathrm{d}}\tilde X=\frac 23, \end{align}

which is equivalent to the critical problem of pinned drops described by (3.2) and (3.3) with $S=S_{\textit{cr}}$ . A comparison of the constant terms yields

(4.26) \begin{align} \frac {L_{\textit{ridge}}^5S}{72V_{\textit{ridge}}^2}=S_{\textit{cr}}. \end{align}

To fully determine the solution, we still need a second relation between $L_{\textit{ridge}}$ and $V_{\textit{ridge}}$ , which turns out to be the apparent contact angle at the advancing contact line, $\theta _{\textit{ap}}\equiv -H'(0)$ . Strictly speaking, $\theta _{\textit{ap}}$ should be determined by matching the macroscopic profile of the ridge to the inner solution near the advancing contact line. However, this approach is infeasible since even the leading-order problem (4.21) should be solved numerically. In the drop sliding problem driven by gravity, Benilov & Benilov (Reference Benilov and Benilov2015) studied the shape of the capillary ridge using empirically $\theta _{\textit{ap}}= [1 + 3c \ln (1/\varLambda ) ]^{1/3}$ , which, however, holds for the low-speed regime without film entrainment only. Instead, we introduce the Cox–Voinov relation $\theta _{\textit{ap}}= [1 + 3c \ln (\alpha H_{\textit{m}}/\varLambda ) ]^{1/3}$ , where $H_{\textit{m}}$ is a macroscopic length and $\alpha$ is a fit parameter (Voinov Reference Voinov1976; Cox Reference Cox1986). As suggested by Ming, Qin & Gao (Reference Ming, Qin and Gao2023), it is appropriate to take $H_{\textit{m}}$ to be the ridge height, which is approximately proportional to $1/S$ , as shown below. Accordingly, we have

(4.27) \begin{align} \theta _{\textit{ap}}=\left (1 + 3 c \ln \frac {\alpha }{S\varLambda } \right )^{\frac {1}{3}}, \end{align}

which further gives rise to

(4.28) \begin{align} H'(0) =\frac {6V_{\textit{ridge}}}{L_{\textit{ridge}}^2}\tilde H'(2) =-\frac {12V_{\textit{ridge}}}{L_{\textit{ridge}}^2}\sqrt {S_{\textit{cr}}} =- \left (1 + 3c \ln \frac {\alpha }{S\varLambda } \right )^{\frac {1}{3}}. \end{align}

Here, we have used $\theta _{{down}}=-\tilde H'(2)$ given in (3.8).

Solving (4.26) and (4.28) simultaneously gives the leading-order predictions for $L_{\textit{ridge}}$ and $V_{\textit{ridge}}$ as

(4.29) \begin{align} L_{\textit{ridge}}=\frac {\left (1 + 3c \ln \dfrac {\alpha }{S\varLambda } \right )^{\frac {2}{3}}}{2S},\quad V_{\textit{ridge}}=\frac {\left (1 + 3c \ln \dfrac {\alpha }{S\varLambda } \right )^{\frac {5}{3}}}{48S^2\sqrt {S_{\textit{cr}}}}. \end{align}

The corresponding maximum height of the capillary ridge at leading order is

(4.30) \begin{align} H_{\textit{max}} = \frac {3V_{\textit{ridge}}}{L_{\textit{ridge}}} \tilde {H}_{\textit{max}}=\frac {1 + 3c \ln \dfrac {\alpha }{S\varLambda }}{8S\sqrt {S_{\textit{cr}}}}\tilde {H}_{\textit{max}}, \end{align}

where $\tilde {H}_{\textit{max}}=0.534$ can be inferred from figure 2. The maximum height $H_{\textit{max}}$ as a function of $S$ is plotted in figure 14(a) for representative values of $\varLambda$ . It can be seen that the numerical results agree perfectly with (4.30) with the fit parameter $\alpha =20$ . Figure 14(b) compares the leading-order ridge profiles with the numerical results and shows notable deviations in the upstream region, which are associated with the non-negligible thickness of the entrained film.

4.3.3. Dimple

Benilov et al. (Reference Benilov, Chapman, McLeod, Ockendon and Zubkov2010) and Benilov & Benilov (Reference Benilov and Benilov2015) noted that, in similar gravity-driven drop motion problem, the dimple region may theoretically consist of an infinite hierarchy of asymptotic zones. Nonetheless, numerical solutions usually reveal only one or two discernible ripples, while the remaining structures are idealised mathematical constructs. Consequently, a detailed investigation of these higher-order features is beyond the scope of this study.

For the dimple region connecting the entrained film and capillary ridge, the governing equation remains the same as (4.4). Similar to the approach outlined in § 4.3.2, we assume the local film thickness is much smaller than the height of the capillary ridge. Away from the dimple, the film thickness approaches a constant $c/S$ to the left and $\sqrt {{8S}/{3}}X^{({2}/{3})}$ to the right. Using the rescalings (4.15), we transform problem into

(4.31) \begin{align} \begin{aligned} &\bar H^2 \bar H''' + \bar H = 1, \quad \bar H(-\infty ) \to 1, \quad \bar H(+\infty ) \to \sqrt {\frac {8}{3}}X^{\frac 23}. \end{aligned} \end{align}

Using the shooting method, we obtain the scaled dimple profile and compare it with the numerical solution, as shown in figure 15. It is seen that the numerical results of the full problem at different $S$ values do collapse onto a single curve using the similarity variables (4.15). The similarity profile captures the oscillatory behaviours but shows large deviations, which are due to the fact that the thicknesses of the dimple and the capillary ridge are not well separated.

Figure 15.

Comparison of the scaled dimple profiles between the numerical results and theoretical solution of problem (4.31), where the slip length is $\varLambda = 10^{-5}$ .

4.4. High shear: unsteady evolution

In the high-shear region, the drop exhibits interesting transient film structures before the final steady state is approached. Taking the case of $S=1$ and $\varLambda = 10^{-5}$ as an example, figure 16(a) shows the drop profiles at several times, while figure 16(b) presents the time evolution of the advancing and receding contact-line velocities. Throughout the evolution, the receding contact-line velocity remains approximately constant and equals the steady-state value, and the drop rear profile also remains unchanged. As the drop evolves, liquid is transported from the capillary ridge into the drop rear, resulting the elongation of the entrained film.

Figure 16.

Unsteady evolution of a drop with $S=1$ and $\varLambda = 10^{-5}$ : (a) drop profiles with the receding contact lines shifted to $X=0$ ; (b) time evolution of the advancing and receding contact-line velocities.

Additionally, we observe the emergence of an even thinner film near the capillary ridge, as shown in figure 16(a). The appearance of this thinner film bears resemblance to the classical thin films observed in dip-coating problems, where a liquid layer is entrained by withdrawing a solid plate from a bath, as first described by Landau & Levich (Reference Landau and Levich1942) and Derjaguin (Reference Derjaguin1943). Similar thin-film structures have also been reported in studies of liquid entrainment associated with moving contact lines (Snoeijer et al. Reference Snoeijer, Delon, Fermigier and Andreotti2006; Gao et al. Reference Gao, Li, Feng, Ding and Lu2016; Xia, Qin & Gao Reference Xia, Qin and Gao2020; Wang et al. Reference Wang, Qin and Gao2025). Such thin films generated by liquid entrainment are known as Landau–Levich–Derjaguin (LLD) films, which scale with the $2/3$ power of the local meniscus velocity relative to the substrate. By analogy, we interpret this thinner film as being entrained from the trailing edge of the capillary ridge through a similar mechanism, although its scaling with the meniscus velocity differs from the classical LLD law (as discussed later). Accordingly, we refer to it as an LLD-like film.

During the evolution, the length of LLD-like film first grows and then decreases. The LLD-like film eventually vanishes and evolves into the steady dimple structure discussed above, inducing a fluctuation in the advancing contact-line velocity at roughly $T=400$ , as shown in figure 16(b). Interestingly, the LLD-like film is characterised by a non-uniform thickness, while the profiles remain approximately aligned along the same line throughout its evolution. These findings suggest that the film is produced through the local dynamics; once generated, the film itself hardly varies with time.

To estimate the thickness of the LLD-like film, we focus on the dynamic meniscus connecting the LLD-like film and the capillary ridge. Since the capillary ridge varies slowly in time (i.e. in a quasi-steady state), the dynamic meniscus can therefore be regarded as moving with the advancing contact-line speed $c_{\textit{ad}v}$ . As mentioned before, the liquid is transported from the capillary ridge into the drop rear, which implies that the flux through the dynamic meniscus cannot be neglected. In a reference frame comoving with the meniscus, the local film is approximately steady and governed by

(4.32) \begin{align} H^3H'''+S(H^2-H_L^2)=c_{\textit{ad}v}(H-H_L), \end{align}

where $H_L$ denotes the thickness of the LLD-like film with its weak non-uniformity neglected.

Since the thickness of the LLD-like film is small compared with that of the ridge, we can also assume that the capillary ridge contacts the substrate at the left end with a zero contact angle. Therefore, similar to § 4.3.3, the right-side far-field boundary condition of the dynamic meniscus reads $H\to \sqrt {8S/3}X^{3/2}$ . The left end of the dynamic meniscus connects to the flat LLD-like film, leading to the left-side far-field boundary condition $H\to H_L$ .

We introduce $H=H_L\bar H$ and $X=H_L\bar X/c_{\textit{ad}v}^{1/3}$ to transform (4.32) and the corresponding boundary conditions into

(4.33) \begin{align} \begin{aligned} \bar H^3\bar H''' +K(\bar H^2-1)=\bar H-1, \quad \bar H \sim \begin{cases} 1, & \text{as} \quad \bar X \to -\infty , \\[10pt] \sqrt {\dfrac {8K}{3}} \bar X^{\frac 32} , & \text{as} \quad \bar X \to +\infty , \end{cases} \end{aligned} \end{align}

where $K=\textit{SH}_L/c_{\textit{ad}v}.$

Problem (4.33) is solved using the shooting method. In particular, we determine $K = 0.2234$ , and hence

(4.34) \begin{align} \begin{aligned} H_L = 0.2234\,\frac {c_{\textit{ad}v}}{S}. \end{aligned} \end{align}

This relation implies that the thickness of the LLD-like film is proportional to the advancing contact-line velocity $c_{\textit{ad}v}$ and inversely proportional to the shear number $S$ , which differs from the classical LLD law. As shown in figure 17, a better agreement with numerical results is achieved when the coefficient is modified as $K=0.36$ . We attribute this deviation to the approximations made in the theoretical derivation. For instance, the height of the capillary ridge is not significantly larger than the thickness of the LLD-like film (only approximately $13\sim 17$ times greater), which violates the assumption of clear scale separation.

Figure 17.

Comparison of LLD-like film thickness with the asymptotic theory $H_L=0.36 c_{\textit{ad}v}/S$ with $\varLambda =10^{-5}$ .