Gravity-driven motion of a highly viscous non-wetting drop on an inclined wall

Abstract

Extensive three-dimensional boundary-integral simulations are presented for the steady-state, low-Reynolds-number motion of a non-wetting deformable drop in another liquid on an inclined solid wall. The drop remains separated from the wall by a lubricating film. The boundary-integral formulation is based on the half-space Green function. The focus is on the challenging case of small tilt angles $\theta$ combined with high drop-to-medium viscosity ratios $\lambda$, when the drop travels with strong hydrodynamical interaction very close to the wall. Simulations to steady state have required ultrahigh drop surface resolutions (to $3\times 10^5$ boundary elements) achieved through multipole acceleration and combined with novel regularization to fully eliminate near-singular behaviour of the double-layer integrals due to small clearance. Non-dimensional drop speed $U$ is presented for $\theta \geq 7.5^\circ$, $\lambda \leq 300$ and in a broad range of Bond numbers $B$, covering from nearly spherical to strongly pancaked drops. The results are consistent with published experiments on liquid–liquid systems. At small $\theta$ and $\lambda \gg 1$, $U$ is a strong, decreasing function of $B$; the asymptotic regime $U\to 0$ at $B\to 0$ is not observed in the simulated range. For small $B$, the dimpled thin-film geometry is insensitive to $\lambda =1\unicode{x2013}300$. For pancaked drops, the lubrication film is much thicker for $\lambda =1$ than for $\lambda \gg 1$ drops. Approximate thin-film uniformity in the drop motion direction is confirmed for pancaked, but not for $B\ll 1$, drops. Kinematics of drop motion shows that neither perfect tank treading, nor perfect rolling can be approached for liquid–liquid systems in the purely hydrodynamical formulation. The methodology is applicable to other problems and can allow for direct inclusion of short-range colloidal forces in three-dimensional boundary-integral simulations.

1. Introduction

Gravity/buoyancy-driven low-Reynolds-number motion of a deformable drop on an inclined wall is a classical fluid mechanics problem, relevant to many technological processes and everyday life, and much effort has been reported in the literature to address this problem experimentally, analytically or numerically. A ubiquitous phenomenon with liquid drops in the air is their ability to stick to non-horizontal solid surfaces due to wetting, as long as the gravity force cannot overcome the contact angle hysteresis (Dussan & Chow 1983; Dussan 1985). However, highly viscous drops (of glycerol) in the air were observed to descend easily on a super-hydrophobic surface, where the hysteresis is made small and the contact angles (advancing and receding) are both close to $180^\circ$ (Richard & Quéré 1999). A similar non-wetting effect is achieved with the so-called liquid marbles, where a hydrophobic powder resides on the drop surface (Aussillous & Quéré 2001). Those authors also offered an approximation for the settling speed of a strongly pancaked drop (the Bond number $B\gg 1$) on the assumption that the drop fluid rolls, in perfect no-slip contact with the wall along the flat spot (in some other literature, this mode of motion for pancaked drops is labelled as tank treading). Close agreement of the theory and experiment confirmed this assumption and justified the neglect of the surrounding air for such physical systems. In the opposite limit $B\ll 1$ of a slightly deformed, highly viscous non-wetting drop rolling without slip in perfect contact with the wall, Mahadevan & Pomeau (1999) derived a simple asymptotic relation (to within a factor) for the settling drop speed $U$ along the wall. Their result $U\sim \sigma B^{-1/2}\sin (\theta )/\mu _d$ (with $\sigma$ being the drop–air surface tension, $\mu _d$ the drop dynamic viscosity and $\theta$ the plane inclination angle to horizontal) stems from the balance of the viscous dissipation rate inside the drop near the contact spot and the rate of change of the drop gravitational energy; again, the effect of the surrounding medium was fully neglected. The missing prefactor in the theory of Mahadevan & Pomeau (1999) was recently derived (Schnitzer, Davis & Yariv 2020) through a far more involved asymptotic analysis in the contact spot area. Note that the well-known contact line paradox (non-integrable stress singularity) did not appear in these theoretical studies due to the special contact angle of $180^\circ$. Surprisingly, the theories of Mahadevan & Pomeau (1999) and Schnitzer et al. (2020) predict the drop settling speed to be a decreasing function of the drop size through $B^{-1/2}$, which was qualitatively confirmed in some range by experiments (Richard & Quéré 1999; Aussillous & Quéré 2001; Aussillous 2002) for a drop (a mixture of water and glycerol) on a super-hydrophobic wall in the air (see also Quéré (2005) for a comprehensive review). There are, however, large fluctuations of experimental data in that range, indicating possible uncontrolled effects of numerous physical factors, e.g. disjoining pressure due to double-layer electrostatic repulsion (Del Castillo et al. 2011), surface roughness (Quéré 2005) etc. making the assumption of drop no-slip rolling in perfect contact with the wall less accurate. In particular, the drop speed, when scaled with $\sigma \sin \theta /\mu _d$, does not collapse on the same curve when the drop viscosity is reduced four times (Aussillous 2002).

Hodges, Jensen & Rallison (2004) developed a leading-order asymptotic analysis for an immiscible, deformable drop moving along a gently inclined wall, based on a different, purely hydrodynamical formulation, with the viscosity $\mu _e$ of the carrier fluid taken into account. In the absence of singular (but short-range) adhesive forces and surface roughness, such a drop is not able to reach perfect contact with the wall, and will always remain separated by a lubrication film due to drop motion, making the contact angle assumptions irrelevant in this case. They demarcated 11 asymptotic regimes (for two-dimensional and realistic three-dimensional drops), depending on the relations between the Bond number, viscosity ratio $\lambda =\mu _d/\mu _e$ (either small or large) and the inclination angle $\theta \to 0$. For a pancaked shape $B\gg 1$ and $(\sin \theta )^{-1/2}\ll \lambda \ll (\sin \theta )^{-2}$, their prediction is tank-treading motion with the settling drop speed $U\approx 4\sigma \sin \theta /(3\mu _d)$, which is identical to that from Richard & Quéré (1999) at the contact angle of $180^\circ$, i.e. the drop speed is unaffected by the carrier fluid viscosity. However, according to Hodges et al. (2004), the increase in $\lambda$ beyond ${\sim }(\sin \theta )^{-2}$ and to $\infty$ should lead to drop sliding as a rigid body, with small but finite film thickness and much larger drop speed (than for tank treading) controlled solely by the carrier fluid viscosity. As noted by Hodges et al. (2004), this behaviour is in stark contrast to the experiments of Richard & Quéré (1999), where the tank-treading regime with $U\approx 4\sigma \sin \theta /(3\mu _d)$ was observed even for $\lambda \approx 50\,000$, far above the theoretical bound of ${\sim }(\sin \theta )^{-2}$. Presumably, no attempt should be made to reconcile these differences: the model of Richard & Quéré (1999) simply postulates the kinematics of tank treading, with perfect drop–wall contact and no slip as a starting point of their analysis. Practically, such a behaviour for highly viscous drops can be achieved in the air environment due to the small–medium viscosity and especially small density. Namely, the drop is able to quickly reach close contact with the wall at the initial stage of drop deposition, after which strong adhesive forces, acting in concert with the surface roughness, can ‘pin’ the drop to the wall, resulting in the tank-treading kinematics.

The purely hydrodynamical formulation, adopted in Hodges et al. (2004) and in the present study, is obviously more realistic when the carrier fluid is liquid, not air (or gas). Indeed, sticking to a gently inclined wall was never reported for a bubble or a drop (made of water, silicon, paraffin oil or glycerol–water mixtures) in many viscous carrier liquids (Aussillous & Quéré 2002; Griggs, Zinchenko & Davis 2009; Rahman & Waghmare 2018). The apparent ‘contact’ angle in these experiments was always $180^\circ$, which points to the existence of a lubricating film and a negligible role of non-hydrodynamic forces that could pin the drop to the substrate.

Griggs, Zinchenko & Davis (2008) and Griggs et al. (2009) also performed the first (and the only so far, to the author's knowledge) three-dimensional boundary-integral (BI) simulations for gravity/buoyancy-driven motion of a drop along the wall to a steady state (or incipient drop breakup) in this formulation. Their analysis was based on the Green function of Blake (1971) for the half-space to reduce the problem to a BI equation on the deforming drop surface only, with direct node-to-node summations to solve this equation at each time step. Although a number of results were obtained for the steady-state drop speed and shape, those were only for generic cases (typically, moderate-to-large inclination angles and only small-to-moderate viscosity ratios $\lambda$; small Bond numbers also had to be excluded due to prohibitive computational difficulties). These shortcomings severely limited the parameter space, not allowing for comparisons with the asymptotic theories. The viscosity ratio $\lambda$ for different liquid–liquid combinations can reach $O(10^2\unicode{x2013}10^3)$, and it would be essential to handle this range by simulations. However, such values of $\lambda$ were totally unreachable with the algorithm of Griggs et al. (2008, 2009), even for not small inclination angles $\theta$. It turns out that, to overcome limitations of these earlier studies and work in a wider parameter range than in Griggs et al. (2008, 2009) (specifically, for most difficult combinations of small tilt angles $\theta$ with high viscosity ratios), it requires a practically all-new and more advanced approach, with extreme surface resolutions and novel desingularization tools. The goal of the present paper is to develop such an approach and apply it for a comprehensive, and highly accurate, study of the steady-state drop motion and related characteristics from first principles.

The problem is formulated in § 2. The solution method is described in § 3, including the BI equation with the half-space Green function and related stresslet, BI desingularizations, drop-mesh control (with different schemes necessary for small tilt angles, depending on the Bond number) and multipole acceleration. Although, for $\theta \ll 1$, the outer surface geometry (outside the near-contact spot) is close to axisymmetric, the lubrication film profile is essentially three-dimensional, and so it was computationally productive not to exploit this distinction, but rigorously consider the whole problem as three-dimensional. A novel and universal full desingularization for double-layer integrals in § 3 is noteworthy, since it resolves a long-standing issue in the BI simulations, is applicable to arbitrary drop shapes and can be potentially used in many other problems with strong, near-contact drop–wall, drop–drop, drop–particle or particle–particle hydrodynamical interactions. Even more crucial for a successful solution here is multipole acceleration, allowing for long-time simulations to steady state with ultrahigh surface discretizations. The complexity of the domain Green function and related stresslet, however, have required considerable algebraic effort for this element (compared with multipole acceleration schemes with the free-space Green function and stresslet). In § 4, we discuss the code validation, beneficial effects of the novel, full BI desingularization and the convergence analysis for the drop speed in the most difficult cases of small tilt angles and high viscosity ratios. In § 5, a systematic and highly accurate analysis is presented for the steady-state drop speed, geometry of the lubrication space and kinematics of the drop fluid motion in a wide range of parameters, including comparisons with the asymptotic theory of Hodges et al. (2004) and the semi-empirical drop-speed relation. The narrative in §§ 4 and 5 is practically independent of the method description in § 3. Conclusions are formulated in § 6, with a discussion of other fluid mechanics problems that can potentially benefit from the methodology developed in this work. The appendices mostly include mathematical details of the algorithm outlined in § 3 and present additional code validations. The numerical code developed for this problem is available from the author upon request.

Because of close drop–wall contact (resulting from small tilt angles) combined with strong drop–wall interactions at $\lambda \gg 1$ and necessitating extreme surface resolutions, a strong preference is given here to the BI method with its sharp interface treatment. It would be highly problematic to use instead more general three-dimensional (3-D) algorithms of computational fluid dynamics, with volume discretization and diffuse interface for the same purpose. With those methods, it would be particularly difficult to discern physical and numerical effects.

2. Problem formulation

Consider the gravity/buoyancy-induced motion of a 3-D deformable drop near an infinite, inclined plane solid wall under creeping flow conditions. Both the drop and external liquids are Newtonian, isothermal and free from surfactants. The drop is non-wetting, and so, it is always separated by a lubrication film from the wall due to drop motion. The wall tilt angle to horizontal is $\theta$, the non-deformed drop radius is $a$, the drop and the medium dynamic viscosities are $\mu _d$ and $\mu _e$, respectively, with the viscosity ratio $\lambda =\mu _d/ \mu _e$. Without a loss of generality, the drop is assumed to be heavier than the medium, with the density contrast $\Delta \rho >0$, and so both phases are in the upper half-space of figure 1.

Figure 1. Schematic for drop motion on an inclined wall.

For steady-state drop sedimentation down the wall, the non-dimensional parameters are $\theta$, $\lambda$ and the Bond number

(2.1)\begin{equation} B=\frac{\Delta\rho g a^2}{\sigma}, \end{equation}

with $\sigma =\mbox {const}$ being the interface surface tension and g the standard gravity acceleration. The main goal is a rigorous study of the drop steady-state speed $U$ along the wall, the geometry of the lubrication space and the mode of drop motion in a wide range of parameters $\theta$, $\lambda$ and $B$. Specifically, we are most interested in extreme viscosity ratios $\lambda \gg 1$ and small inclination angles – the most challenging combinations that could not be approached by the more basic BI algorithms of Griggs et al. (2008, 2009) and have required the development of many novel tools herein, both numerical and semi-analytical, to handle unavoidable extreme surface resolutions in dynamical simulations. The present solutions also cover a wide range of $B$ – from nearly spherical to strongly pancaked drops (but below critical for drop breakup). As in Griggs et al. (2008, 2009), the drop velocity is defined as the volume average of the fluid velocity inside the drop reduced to surface integration by Gauss’ theorem, and so only the interfacial velocity has to be determined. Below, where necessary to distinguish between the transient and steady-state values for the drop speed (or other quantities), the index ‘st’ will be added to steady-state values.

A convenient, although somewhat empirical, reference velocity scale is

(2.2)\begin{equation} U_{ref}=\frac{2}{9}\frac{(\lambda+1)}{(\lambda+2/3)} \frac{\Delta\rho g a^2}{\mu_e} \sin\theta, \end{equation}

which is the steady-state speed of an isolated spherical drop under reduced gravity $g\sin \theta$. To make the equations and results below non-dimensional, all the velocities will be scaled with $U_{ref}$, lengths with $a$ and times with $a/U_{ref}$. Note that a different velocity scale, $\sigma /\mu _e$, is used in the asymptotic theory of Hodges et al. (2004); the ratio of the two scales is $O(B\sin \theta )$. The scale (2.2) was chosen in the present work, since it makes the non-dimensional, steady-state drop speed $U$ much less sensitive to $\theta$ and $B$.

3. Method

3.1. Boundary-integral equation

A fixed, right Cartesian coordinate system $(x_1,x_2, x_3)$ is chosen with the origin on the wall, the $x_3$ axis normal to it and the $x_2$ axis along the projection of gravity vector $\boldsymbol {g}$ on the wall (figure 1). As in Griggs et al. (2008, 2009), the BI equation for the non-dimensional Stokes flow problem at any instantaneous, transient drop–wall configuration is based on the Green tensor $\boldsymbol {G}(\boldsymbol {x};\boldsymbol {y})$ and related fundamental stresslet $\boldsymbol {\tau }(\boldsymbol {x};\boldsymbol {y})$ for the half-space ($x_3, y_3>0$). A crucial advantage of using these functions (instead of their simple free-space counterparts $\boldsymbol {G}^{FS}$ and $\boldsymbol {\tau }^{FS})$ is that they work to exclude the wall BIs due to no slip, and so the problem can be formulated as a BI equation on the drop surface $S$ only for the fluid velocity $\boldsymbol {u}(\boldsymbol {x})$. Such an equation corresponds to the original, infinite half-space solution domain. An alternative of additionally discretizing the wall (with new unknown distributions) and embedding the drop–wall configuration into a finite computational domain (with slowly decaying finite-size effects) would make the solution much less efficient and accuracy far more difficult to achieve. Plus, exclusion of the wall from the equations greatly facilitates the logic of the multipole-accelerated BI algorithm paramount in the present superhigh resolution simulations.

The second-rank Green tensor $\boldsymbol {G}$ was derived by Blake (1971) and reviewed in Pozrikidis (1992), who also presented $\boldsymbol {\tau }$. With the normalization used herein, $\boldsymbol {G}^k= (G^k_1,G^k_2,G^k_3)(\boldsymbol {x};\boldsymbol {y})$ (for any fixed $k=1,2,3$) is, by definition, the unit-viscosity Stokes flow velocity at $\boldsymbol {x}$ due to the point force $-\boldsymbol {e}_k$ applied to the fluid at $\boldsymbol {y}$, subject to the no-slip boundary conditions $G^k_j(\boldsymbol {x};\boldsymbol {y})=0$ $(j=1,2,3)$ when $x_3=0$. The corresponding stress tensor components at $\boldsymbol {x}$ are $\tau ^k_{ij}(\boldsymbol {x};\boldsymbol {y})$.

Based on general theory (Rallison & Acrivos 1978; Pozrikidis 1992), the non-dimensional BI equation for $\boldsymbol {u}(\kern0.09em \boldsymbol {y})$ on the drop surface can be written as

(3.1)\begin{equation} \boldsymbol{u}(\kern0.09em \boldsymbol{y})= \boldsymbol{F}(\kern0.09em \boldsymbol{y})+\varkappa\left[2\int_S Q_i(\boldsymbol{x})\boldsymbol{\tau}_{ij} (\boldsymbol{x};\boldsymbol{y})n_j(\boldsymbol{x})\,\mbox{d}S_x+\boldsymbol{u}'(\kern0.09em \boldsymbol{y})\right],\quad \boldsymbol{y}\in S. \end{equation}

Here, $\varkappa =(\lambda -1)/(\lambda +1)$, $\boldsymbol {\tau }_{ij}=(\tau _{ij}^1,\tau _{ij}^2,\tau _{ij}^3)$, $\boldsymbol {n}(\boldsymbol {x})$ is the outward unit normal to $S$, and $\boldsymbol {u}'(\boldsymbol {x})$ is the projection of $\boldsymbol {u}(\boldsymbol {x})$ on the subspace of rigid-body motions; this projection is easy to calculate without the Gram–Schmidt orthogonalization, as noted in Zinchenko, Rother & Davis (1997). The form (3.1), with $\boldsymbol {Q}(\boldsymbol {x})= \boldsymbol {u}(\boldsymbol {x})-\boldsymbol {u}'(\boldsymbol {x})$ under the double-layer integral and the added-back term $\boldsymbol {u}'(\kern0.09em \boldsymbol {y})$, serves to generally reduce the magnitude of the integrand and increase the efficiency of the multipole-accelerated solution for (3.1). With the reference capillary number $Ca_{ref}=\mu _e U_{ref}/\sigma$, the non-dimensional inhomogeneous term, due to stress jump on the interface, can be written as

(3.2)\begin{equation} \boldsymbol{F}(\kern0.09em \boldsymbol{y})=\frac{2}{(\lambda+1)}\frac{B}{Ca_{ref}} \int_S\left[\frac{2\Delta k(\boldsymbol{x})}{B}+ \Delta x_3\cos\theta -\Delta x_2\sin\theta\right] n_j(\boldsymbol{x})\boldsymbol{G}_j(\boldsymbol{x};\boldsymbol{y})\,\mbox{d}S_x. \end{equation}

Here, $\boldsymbol {G}_j=(G_j^1,G_j^2,G_j^3)$ (note the difference from the above definition of $\boldsymbol {G}^k$), $\Delta k(\boldsymbol {x})= k(\boldsymbol {x})-\langle k\rangle$, $\Delta \boldsymbol {x}=\boldsymbol {x}-\boldsymbol {x}^c$, $k(\boldsymbol {x})=(k_1+k_2)/2$ is the local surface curvature (half-sum of principal curvatures) with the average $\langle k\rangle$ over the whole surface and $\boldsymbol {x}^c$ is the drop surface centroid; the centroid components and $\langle k\rangle$ are subtracted to generally reduce the magnitude of the integrand in (3.2) without changing the integral.

A major challenge in the simulations with $\theta \ll 1$ and $\lambda \gg 1$ was how to avoid non-convergence of iterations for the BI equation (3.1), stalling dynamical simulation well before the steady-state drop speed is reached, even with full desingularization (see below) and a powerful orthogonalized GMRES (Generalized Minimal Residual) method (Arnoldi iteration). For $\lambda \gg 1$, $(\lambda -1)/(\lambda +1)$ is close to the theoretical marginal spectral value $\varkappa =1$ of the right-hand side operator (3.1), and so small discretization errors growing in dynamical simulation sufficiently distort the spectrum and cause this non-convergence. Theoretically, partial Wielandt deflation (e.g. Kim & Karilla 1991; Pozrikidis 1992) purges unity from the spectrum by working with an auxiliary vector field $\boldsymbol {\omega }= \boldsymbol {u} -\varkappa \boldsymbol {u}'$. The equation for $\boldsymbol {\omega }(\kern0.09em \boldsymbol {y})$ is the same as (3.1), with $\boldsymbol {Q}(\boldsymbol {x})=\boldsymbol {\omega }(\boldsymbol {x})-\boldsymbol {\omega }'(\boldsymbol {x})$ (which is still $\boldsymbol {u}-\boldsymbol {u}'$), but without the added-back rigid-body projection term in the brackets. The physical velocity $\boldsymbol {u}$ is recovered as $\boldsymbol {\omega }+\varkappa \boldsymbol {\omega }'/(1-\varkappa )$. Against expectations, partial deflation never resolved the issue (presumably, because the spectrum is nearly continuous for $\theta \ll 1$, and purging just one marginal value does little); instead, using high/ultrahigh resolutions (through multipole acceleration) was imperative to make BI iterations convergent all the way to steady state, for both non-deflated and partially deflated versions. However, partial deflation provided a somewhat more accurate steady-state speed for pancaked drops and was the method of choice in that range; for $B\leq O(1)$, both the non-deflated (3.1) and partially deflated forms could be used with indistinguishable results. Note finally that full deflation (which also eliminates -1 from the spectrum and is often practiced in BI simulations) was irrelevant in our case and even detrimental, making the steady-state speed of a pancaked drop much slower convergent with respect to resolutions when $\theta \ll 1$ and $\lambda \gg 1$.

The full Green tensor and fundamental stresslet are split as

(3.3a,b)\begin{equation} \boldsymbol{G}= \boldsymbol{G}^{FS}+\boldsymbol{G}^{C},\quad \boldsymbol{\tau}= \boldsymbol{\tau}^{FS}+\boldsymbol{\tau}^{C}, \end{equation}

where $\boldsymbol {G}^{C}$ and $\boldsymbol {\tau }^{C}$ are the wall-correction parts, with explicit expressions given in Blake (1971) and Pozrikidis (1992) (note a different normalization used herein in the definitions of $\boldsymbol {G}$ and $\boldsymbol {\tau }$). Instead of the cumbersome forms for $\boldsymbol {G}^{C}$ and $\boldsymbol {\tau }^{C}$, it is more relevant to (3.1) and (3.2) to know how these tensors act on arbitrary vectors $\boldsymbol {W}$ and $\boldsymbol {Q}$. Namely,

(3.4)\begin{align} & {-}8{\rm \pi} W_j\boldsymbol{G}^{C}_j(\boldsymbol{x};\boldsymbol{y})= \left[\left(\frac{6x_3y_3}{R^2}-1\right)(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{R}) +2y_3 W_3\right]\frac{\boldsymbol{R}^{IM}}{R^3} \nonumber\\ &\quad -\left(1+\frac{2x_3y_3}{R^2}\right)\frac{\boldsymbol{W}^{IM}}{R} -2\left[\frac{x_3(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{R})}{R^2}+W_3\right]\frac{\boldsymbol{e}_3}{R} , \end{align}

and

(3.5)\begin{align} & \frac{4{\rm \pi}}{3}Q_iW_j\boldsymbol{\tau}^{C}_{ij}(\boldsymbol{x};\boldsymbol{y}) = \left[2y_3^2(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{Q})+\left(\frac{10x_3y_3}{R^2}-1\right) (\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{R})(\boldsymbol{Q}\boldsymbol{\cdot}\boldsymbol{R})\right]\frac{\boldsymbol{R}^{IM}}{R^5} \nonumber\\ &\quad -\frac{2x_3 y_3}{R^5}[(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{R})\boldsymbol{Q}^{IM}+ (\boldsymbol{Q}\boldsymbol{\cdot}\boldsymbol{R})\boldsymbol{W}^{IM}] -\frac{2x_3}{R^5}(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{R}) (\boldsymbol{Q}\boldsymbol{\cdot}\boldsymbol{R})\boldsymbol{e}_3. \end{align}

Here, $\boldsymbol {R}=\boldsymbol {x}- \boldsymbol {y}^{IM}$, and the superscript ${IM}$ stands for the mirror image of a vector or a point with respect to the wall. Equations (3.4) and (3.5) are complemented by standard FS contributions

(3.6)\begin{equation} {-}8{\rm \pi} W_j\boldsymbol{G}^{FS}_j(\boldsymbol{x};\boldsymbol{y})=\frac{\boldsymbol{W}}{r}+ \frac{(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{r})\boldsymbol{r}}{r^3}, \end{equation}

and

(3.7)\begin{equation} \frac{4{\rm \pi}}{3}Q_iW_j\boldsymbol{\tau}^{FS}_{ij}(\boldsymbol{x};\boldsymbol{y}) = \frac{(\boldsymbol{W}\boldsymbol{\cdot}\boldsymbol{r})(\boldsymbol{Q}\boldsymbol{\cdot}\boldsymbol{r})\boldsymbol{r}}{r^5}, \end{equation}

with $\boldsymbol {r}=\boldsymbol {x}-\boldsymbol {y}$. After some algebra, the combination of (3.5) and (3.7) can be transformed to match (7) of Griggs et al. (2009) (except for the typo in the second line of their equation: the closing bracket ] must be moved to the next line). Computationally, however, our (3.5) is slightly more economical.

3.2. Desingularization of boundary integrals

3.2.1. Single-layer desingularization

Let, for brevity, $\tilde {f}(\boldsymbol {x})$ be the expression in the square brackets of (3.2). The true singularity of the integrand in (3.2), when $\boldsymbol {x}=\boldsymbol {y}$, comes from $\boldsymbol {G}^{FS}$ and is eliminated in a standard way, by subtracting $\tilde {f}(\kern0.09em \boldsymbol {y})$ from $\tilde {f}(\boldsymbol {x})$. The wall-correction part (3.4) is finite for all $\boldsymbol {x}, \boldsymbol {y}\in S$, but it can be quite large, ${\sim }1/R$, when $\boldsymbol {x}$ is close to the mirror image of $\boldsymbol {y}$; this is the case for a drop moving almost in contact with the wall, of most interest in the present study. However, such near singularity only comes into play for mesh nodes $\boldsymbol {y}$ in the vicinity of the wall. Let $\boldsymbol {x}^\ast \in S$ be the mesh node on $S$ closest to $\boldsymbol {y}^{IM}$ (figure 1). By the continuity equation for $\boldsymbol {G}^{FS}$ and $\boldsymbol {G}^{C}$, the integral (3.2) can be calculated as

(3.8)\begin{align} & \int_S \tilde{f}(\boldsymbol{x})n_j(\boldsymbol{x})\boldsymbol{G}_j(\boldsymbol{x};\boldsymbol{y})\,\mbox{d}S_x= \int_S [\tilde{f}(\boldsymbol{x})-\tilde{f}(\kern0.09em \boldsymbol{y})] n_j(\boldsymbol{x})\boldsymbol{G}_j^{FS}(\boldsymbol{x};\boldsymbol{y})\,\mbox{d}S_x \nonumber\\ &\quad +\int_S[\tilde{f}(\boldsymbol{x})-\varTheta(\kern0.09em \boldsymbol{y})\tilde{f}(\boldsymbol{x}^\ast)] n_j(\boldsymbol{x};\boldsymbol{y})\boldsymbol{G}_j^{C}(\boldsymbol{x};\boldsymbol{y})\,\mbox{d}S_x. \end{align}

Here, $\varTheta (\kern0.09em \boldsymbol {y})$ is a barrier function, which is close to 1 for $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|\ll 1$ and, by definition, set to zero for $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|> h_o$, where $h_o$ is a moderately small threshold. In the present work, $h_o$ was fixed at 0.25, and one of the possible forms for $\varTheta$ was used

(3.9)\begin{equation} \varTheta(\kern0.09em \boldsymbol{y})= 1-(\|\boldsymbol{y}^{IM}-\boldsymbol{x}^\ast\|/ h_o)^4. \end{equation}

When $\varTheta \approx 1$, the subtracted term in the last integral (3.8) effectively cancels the near singularity of $\boldsymbol {G}^{C}$; for $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|> h_o$, such subtraction would not be beneficial. Also, small support for $\varTheta$ considerably improves the efficiency of our multipole-accelerated calculation of the desingularized form (3.8).

3.2.2. Free-space double-layer desingularization

The FS contribution to the integral (3.1) is desingularized in a standard way. The explicit form (3.7) yields

(3.10)\begin{align} \int_S Q_i(\boldsymbol{x})\boldsymbol{\tau}^{FS}_{ij} (\boldsymbol{x};\boldsymbol{y})n_j(\boldsymbol{x})\,\mbox{d}S_x= \frac{3}{4{\rm \pi}}\int_S\frac{\{[\boldsymbol{Q}(\boldsymbol{x})-\boldsymbol{Q}(\kern0.09em \boldsymbol{y})]\boldsymbol{\cdot}\boldsymbol{r}\} [\boldsymbol{r}\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})]\boldsymbol{r}}{r^5} \mbox{d}S_x + \frac{1}{2}\boldsymbol{Q}(\kern0.09em \boldsymbol{y}). \end{align}

Here, the subtraction of $\boldsymbol {Q}(\kern0.09em \boldsymbol {y})$ fully suppressed the singularity and makes the right-hand side integrand finite, since $\boldsymbol {r}\boldsymbol {\cdot }\boldsymbol {n}(\boldsymbol {x})=O(r^2)$ for any close pair of points $\boldsymbol {x}$ and $\boldsymbol {y}$ on a smooth surface.

3.2.3. Wall-correction double-layer desingularization

It is non-trivial to fully desingularize the integrand in the wall-correction contribution to the double-layer integral (3.1). One can, of course, follow the logic of (3.8) and subtract $Q_i(\boldsymbol {x}^\ast )$ from $Q_i(\boldsymbol {x})$ in such an integral (without an added-back term), but, for a drop in near contact with the wall, such a subtraction only reduces the integrand near singularity from $O(1/R^2)$ to $O(1/R)$. The reason is that $\boldsymbol {R}\boldsymbol {\cdot }\boldsymbol {n}(\boldsymbol {x})$ is still $O(R)$ for points $\boldsymbol {x}$ and $\boldsymbol {y}$ (on $S$) close to each other and the wall. The lack of full desingularization could be offset in the present work to a large extent by extreme resolutions (possible due to multipole acceleration) to simulate the steady-state drop velocity, but such simulations underperformed, with much longer relaxation to the steady state, excessive number of BI iterations per time step and usually quite inaccurate thin-film profiles in the most difficult runs; without full desingularization, some high-resolution runs still crashed with drop–wall overlap, before reaching the steady state.

It is relevant to note a similar, long-standing issue with full desingularization of the FS double-layer integrals (3.10) in other problems (e.g. close interaction between drops, or between drops and solid particles in unbounded space or a periodic box), where the observation point $\boldsymbol {y}$ can be slightly outside the integration surface. Again, using, for subtraction, $\boldsymbol {Q}$ at the surface point nearest to $\boldsymbol {y}$ only reduces the integrand singularity to $O(1/r)$. To the author's knowledge, neither of the works aimed at desingularized BI equations (e.g. Bazhlekov, Anderson & Meijer 2004; Klaseboer, Sun & Chan 2012; Farutin, Biben & Misbah 2014) could fully eliminate near-singular behaviour of the FS integrals (3.10), when $\boldsymbol {y}$ is off the integration surface. Zinchenko & Davis (2002) constructed the subtracted quantity $\boldsymbol {Q}^\ast$ differently from the solution of a special variational problem, which improved the spectral properties of the discretized double layer and allowed for highly concentrated emulsion flow simulations with moderate viscosity contrast, still without full desingularization. This technique, however, was found by Zinchenko & Davis (2005) to lose advantage with the increase of resolutions and was also not powerful enough for the present study. In contrast, in the high-order near-singularity subtraction scheme (Zinchenko & Davis 2006, 2008), $\boldsymbol {Q}({\boldsymbol {x}})$ is locally approximated near $\boldsymbol {x}^\ast$ as a constant plus a linear function of the coordinates in the tangential plane, and both terms are subtracted from $\boldsymbol {Q}(\boldsymbol {x})$ to fully eliminate the near singularity of the integrand. Their technique (applied to drop squeezing through constrictions) was only suitable for integration over a solid particle surface of canonic shape (sphere, spheroid, torus), when the added-back integral allows for analytical treatment. To generalize such an approach to other shapes, Gissinger, Zinchenko & Davis (2021) opted for direct numerical evaluation of the added-back integral on a superfine auxiliary mesh, but their method is still for integration over fixed solid surfaces and could not be used for a deformable drop surface in dynamical simulations.

The present work offers, for the first time, a recipe for full removal of near-singular behaviour of the double-layer integrand, suitable for any smooth, closed integration surface and any type of the fundamental stresslet $\boldsymbol {\tau }(\boldsymbol {x};\boldsymbol {y})$ (free space or wall corrected). In the present context of the wall-bounded geometry, the general identity derived in Appendix A applies to any pair of points $\boldsymbol {y}, \boldsymbol {x}^\ast \in S$ and any vector quantity $\boldsymbol {Q}^\ast$

(3.11) \begin{align} & \int_S \boldsymbol{Q}(\boldsymbol{x})\boldsymbol{\cdot}(\boldsymbol{\tau}^k)^{C}(\boldsymbol{x};\boldsymbol{y}) \boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})\,\mbox{d}S_x \nonumber\\ &\quad =\int_S \{\boldsymbol{Q}(\boldsymbol{x})-\boldsymbol{Q}^\ast{-}{\boldsymbol{\mathcal L}}(\boldsymbol{x}-\boldsymbol{x}^\ast) +[(\boldsymbol{x}-\boldsymbol{x}^\ast)\boldsymbol{\cdot}\boldsymbol{n}^\ast]\,\boldsymbol{f}^\ast\} \boldsymbol{\cdot}(\boldsymbol{\tau}^k)^{C}(\boldsymbol{x};\boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})\,\mbox{d}S_x \nonumber\\ &\qquad +\int_S[\,\boldsymbol{f}(\boldsymbol{x})-(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{n}^\ast)\boldsymbol{f}^\ast{-} (\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{f}^\ast)\boldsymbol{n}^\ast]\boldsymbol{\cdot}(\boldsymbol{G}^k)^{C}(\boldsymbol{x};\boldsymbol{y})\,\mbox{d}S_x \nonumber\\ &\qquad -{\mathcal{L}}^i_i\int_S\left\{[n_k-n_k^\ast(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{n}^\ast)] \varphi^k(\boldsymbol{x};\boldsymbol{y}) + n_k^\ast[(\boldsymbol{x}-\boldsymbol{x}^\ast)\boldsymbol{\cdot}\boldsymbol{n}^\ast] \frac{\partial \varphi^k(\boldsymbol{x};\boldsymbol{y})}{\partial n}\right\}\mbox{d}S_x , \end{align}

for $k=1,2,3$ (no summation over $k$ in the last line). Here, $(\boldsymbol {\tau }^k)^{C}$ is the wall-correction part of the second-rank tensor $\{\tau _{ij}^k\}$ for a fixed $k$; likewise, $(\boldsymbol {G}^k)^{C}$ is the wall-correction part of the Green vector $\boldsymbol {G}^k$. Further, ${\boldsymbol{\mathcal L}}(\boldsymbol {x}-\boldsymbol {x}^\ast )$ is a linear vector field

(3.12)\begin{equation} {\mathcal{L}}_i= {\mathcal{L}}_i^j (x_j-x^\ast_j), \end{equation}

with a constant matrix of coefficients ${\mathcal {L}}_i^j$ and the trace ${\mathcal {L}}_i^i$. The only limitation is ${\mathcal {L}}_i^j n_j(\boldsymbol {x}^\ast )=0$, i.e. (3.12) only depends on the projection of $\boldsymbol {x}-\boldsymbol {x}^\ast$ on the tangential plane at $\boldsymbol {x}^\ast$. The related vector field $\boldsymbol {f}(\boldsymbol {x})$ is given in coordinates as

(3.13)\begin{equation} f_j(\boldsymbol{x})= ({\mathcal{L}}_j^i + {\mathcal{L}}_i^j)n_i(\boldsymbol{x}). \end{equation}

For brevity, $\boldsymbol {n}$ in (3.11) is used for $\boldsymbol {n}(\boldsymbol {x})$, while the asterisk denotes the values of $\boldsymbol {n}$ and $\boldsymbol {f}$ at $\boldsymbol {x}^\ast$. Finally, the auxiliary functions $\varphi ^k$, specific for the half-space geometry, are

(3.14)\begin{equation} \varphi^k(\boldsymbol{x};\boldsymbol{y})={-}\frac{1}{4{\rm \pi}} \left(\frac{1}{R}\mp \frac{2y_3R_3}{R^3}\right), \end{equation}

where the upper sign is for $k=1$ and $2$, and lower sign is for $k=3$; the normal derivative in the last line of (3.11) is with respect to $\boldsymbol {x}$.

For full desingularization of (3.11), $\boldsymbol {x}^\ast$ is chosen again as the mesh node nearest to $\boldsymbol {y}^{IM}$ (figure 1), and $\boldsymbol {Q}^\ast$ is set to $\varTheta (\kern0.09em \boldsymbol {y})Q(\boldsymbol {x}^\ast )$, while the construction of ${\boldsymbol{\mathcal L}}(\boldsymbol {x}-\boldsymbol {x}^\ast )$ follows Zinchenko & Davis (2006). Namely, in the intrinsic coordinate system $(x'_1, x'_2, x'_3)$ centred at $\boldsymbol {x}^\ast$ and with the $x'_3$-axis along $\boldsymbol {n}(\boldsymbol {x}^\ast )$, ${\boldsymbol{\mathcal L}}$ is first sought as a linear function of $x'_1$ and $x'_2$ for the least-squares fit to $\boldsymbol {Q}(\boldsymbol {x})-\boldsymbol {Q}(\boldsymbol {x}^\ast )$ in the set of mesh nodes (mesh triangle vertices) directly connected to $\boldsymbol {x}^\ast$, with subsequent transformation (3.12) to global coordinates. This way, the whole expression in the braces of the second line of (3.11) is $O(\|\boldsymbol {x}-\boldsymbol {x}^\ast \|^2)+O(\|\boldsymbol {x}^\ast -\boldsymbol {y}^{IM}\|^4)$ for $\boldsymbol {x}\approx \boldsymbol {x}^\ast$ and effectively cancels the near singularity of $(\boldsymbol {\tau }^k)^{C}$. Due to $\boldsymbol {n}^\ast \boldsymbol {\cdot } \boldsymbol {f}^\ast =0$, the integrand in the third line of (3.11) is also non-singular, since the $O(1/R)$ near singularity of $(\boldsymbol {G}^k)^{C}$ is cancelled by the term in brackets; so is the last integrand in (3.11).

With high surface resolutions, full desingularization (3.11) only makes sense for $\boldsymbol {y}$ very close to the wall, and so a 2-tiered scheme is used in the algorithm. For $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|>h_o$, the integral (3.11) is handled in its original, left-hand side form. For $h_1<\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|< h_o$, only a leading-order subtraction of $\varTheta (\kern0.09em \boldsymbol {y})Q(\boldsymbol {x}^\ast )$ (akin to (3.8)) is used, ignoring all other terms/integrals associated with ${\boldsymbol{\mathcal L}}$; these terms/integrals are additionally included only for $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|< h_1$. Here, $h_1(\kern0.09em \boldsymbol {y})=(1\unicode{x2013}1.2)\varDelta$, where $\varDelta$ is the average length of the mesh edges (see below) emanating from node $\boldsymbol {y}$. Desingularization (3.11) can be organized without slowing down the BI iterations (§ 3.4), but small values of $h_1$, still beneficial, help to greatly reduce the CPU cost of the pre-iterative part associated with (3.11).

3.3. Drop surface discretization and dynamic mesh control

In the present work, requiring long-time simulations for a $\lambda \gg 1$ drop in very close contact with the wall, it was highly beneficial to use unstructured drop surface triangulations, with a fixed number of mesh nodes and fixed mesh topology in each run to the steady state. That way, surface interpolations and other undesirable, sudden mesh changes (which could disrupt iterative solution of the BI problem in the presence of strong drop–wall interaction through a thin lubrication film) are avoided altogether. Accordingly, the problem symmetry about the $x_1=0$ plane (drawn through the drop centroid) was not exploited in the solution (but naturally achieved with high resolutions), since our triangular mesh nodes do not possess such symmetry. Forcing this symmetry (only to achieve less than 2-fold gain in the code speed) would require invoking surface interpolations at each time step. However, three distinct meshing methods had to be used depending on the Bond number, as briefly outlined below; additional details are given in Appendix B.

3.3.1. Small $B$: projective meshing

An adaptation (and some simplification) of the projective mesh approach, originally developed by Zinchenko & Davis (2005) for glancing, near-contact interactions of two slightly deformable drops in free space, was found to be most suitable here to handle a slightly deformable drop near the wall. At each time moment $t$, the projection centre $\boldsymbol {O}(t)$ is placed between the drop surface centroid $\boldsymbol {x}^c$ and the wall (depending on $x_3^c$). Let $\varOmega (t)$ be a unit sphere centred at $\boldsymbol {O}(t)$. A nearly uniform, unstructured mesh with a prescribed number $N_\triangle$ of triangular elements is first constructed on $\varOmega (0)$ by now standard methods, starting from either a regular icosahedron or dodecahedron followed by a series of refinements (Zinchenko et al. 1997). The possibilities include $N_\triangle = 34\,560$, 46 080, 61 440, 77 760, 138 240, 245 760, 327 680 (abbreviated below as 35 K, etc.) used in the present simulations; the maximum-to-minimum mesh edge ratio for each of these triangulations is within 1.19. For a drop starting from spherical in close contact with the wall, the initial drop surface triangulation is simply a projection of the $\varOmega (0)$ mesh from $\boldsymbol {O}(0)$ (figure 2a). As the drops deforms and moves, so does the mesh on $S(t)$, and the projection centre $\boldsymbol {O}(t)$ evolves as well. The mesh nodes $\boldsymbol {x}^j$ are projected back from $S(t)$ onto the current $\varOmega (t)$ (figure 2b) to form a ‘parametric mesh’ on $\varOmega (t)$. This mapping is used for non-iterative calculation of the normal vector $\boldsymbol {n}(\boldsymbol {x}^j)$ and local curvature $k(\boldsymbol {x}^j)$ in the mesh nodes on $S$. The parametric, nearly uniform mesh on $\varOmega (t)$ also facilitates accurate calculation of the regularized BIs on $S(t)$, by transforming this operation to integration over $\varOmega (t)$. Compared with Zinchenko & Davis (2005), such integration is improved here by handling the corresponding mesh triangles on $\varOmega (t)$ as geodesic without the need for additional mesh nodes. By reassigning contributions from such mesh triangles on $\varOmega (t)$ to mesh vertices $\boldsymbol {x}^j$ on $S$, all smooth function integrations are economically performed as

(3.15)\begin{equation} \int_S \varphi(\boldsymbol{x})\,\mbox{d}S\approx \sum_{\boldsymbol{x}^j\in S} \varphi(\boldsymbol{x}^j)\triangle S_j, \end{equation}

where $\triangle S_j$ can still be called a mesh area associated with node $\boldsymbol {x}^j$. The BI solution for the interfacial velocity $\boldsymbol {u}$ only imposes a constraint $\boldsymbol {V}_j\boldsymbol {\cdot }\boldsymbol {n}(\boldsymbol {x}^j)=\boldsymbol {u}(\boldsymbol {x^j})\boldsymbol {\cdot }\boldsymbol {n}(\boldsymbol {x}^j)$ on the mesh node velocities $\boldsymbol {V}_j=\mbox {d}\kern0.07em\boldsymbol {x}^j/\mbox {d}t$. Remarkably, it is still possible to add tangential motion of mesh nodes on $S$ to make the parametric mesh stationary (in the reference frame moving with $\boldsymbol {O}(t)$) for the entire simulation, from $t=0$ to steady state. The rule for placing the projection centre $\boldsymbol {O}(t)$ controls near-contact mesh adaptivity (Appendix B).

Figure 2. A 2-D sketch of gap-adaptive, projective drop surface meshing near the wall. (a) For a drop starting from spherical, the initial meshing is obtained by radial projection of a uniform unit sphere meshing from the centre $\boldsymbol {O}$ placed between the wall and the drop centroid by the rule (3.15). (b) For $t>0$, migration of mesh nodes on $S(t)$, compatible with the BI solution, allows for a parametric mesh on $\varOmega (t)$ to remain stationary in the reference frame moving with $\boldsymbol {O}(t)$. A small size of the unit sphere $\varOmega$ relative to the drop is chosen just for presentation.

Figure 3(a,b) gives an example of projective meshing in a simulation with $\lambda =30$, $\theta =7.5^\circ$, $B=0.125$. A small number of elements $N_\triangle =8640$ was intentionally chosen just to facilitate mesh visualization. Due to insufficient resolution, this run could only approach an accurate, global steady-state drop shape shown in (a), but stalled shortly thereafter with drop–wall overlap and non-convergence of BI iterations, well before the steady state for drop velocity could be reached. Nevertheless, this example demonstrates the high quality of projecting meshing, with almost equilateral mesh triangles and smooth mesh transition from the near-contact zone to the outer region.

Figure 3. A simulation with $\lambda =30$, $\theta =7.5^\circ$, $B=0.125$ and $N_\triangle =8640$ using dynamic projective meshing. (a) Side view along the $x_1$-axis. (b) Diagonal view from the drop bottom showing the dimple and the near-contact mesh adaptivity (with maximum-to-minimum mesh edge ratio $r_{max}\approx 4$). Only global steady-state drop shape (not the drop speed) could be approached in this low-resolution run.

3.3.2. Intermediate $B$: $n_3$-based passive mesh stabilization

Projective meshing, although a highly accurate method, is obviously limited to small-to-moderate deformations. With this method, when the drop becomes significantly pancaked and the size of the near-contact spot comparable to the drop size, the top portion of the drop can unnecessarily receive higher resolution than the outer part of the spot. Another difficulty is the loss of smoothness of the projective mapping onto the unit sphere near the rim of the spot in this case. An alternative, robust approach for moderate $B$ should still be adaptive, with higher concentration of mesh nodes in near contact with the wall. Following the version of a passive mesh stabilization technique from Zinchenko & Davis (2013), the node velocities $\boldsymbol {V}_i= \mbox {d}\kern0.07em\boldsymbol {x}^i/\mbox {d}t$ to be used in the drop shape update are required, at each time step, to minimize a form of the ‘kinetic mesh energy function’ under the normal velocity constraints dictated by the BI solution. This function, quadratic in $\boldsymbol {V}_1$, $\boldsymbol {V}_2\cdots$, contains the rates of change of deviations between the mesh edges $\|\boldsymbol {x}_{ij}\|$ and their local target values $h_{ij}$, and also has a term to resist mesh triangle quality deterioration. The novelty of the present case is in how to construct $h_{ij}$. Curvature adaptation of Zinchenko & Davis (2013) is not the right criterion here, since local curvature in the near-contact spot is overall much smaller than on the rest of the drop. Nor was it possible to develop a stable and accurate meshing scheme with fixed topology for the present simulations based on the local drop–wall clearance, since such schemes are overly sensitive to local features which do not necessarily have a significant global effect. Note instead that, in the coordinate system of figure 1, $n_3\approx -1$ only in the spot area of the drop surface. Accordingly, with an empirical coefficient $\alpha$ set close to 1, a simple form $h(\boldsymbol {x})\sim [1+\alpha n_3(\boldsymbol {x})]^{1/2}$ for the target mesh size $h$ near $\boldsymbol {x}$ provided robust meshing of the entire surface, adaptive to the near-contact spot but not requiring knowledge of the spot size or surface clearance. The degree of adaptivity is controlled by $\alpha$ (with default 0.7). Details of this technique are described in Appendix B.

Figure 4(a,b) gives an example of $n_3$-based meshing in a simulation with ${\lambda =30}$, ${\theta =7.5^\circ}$, ${B=1}$ and intentionally small number of elements $N_\triangle =8640$ just to facilitate mesh visualization, with required adaptivity to the near-contact spot, nearly equilateral mesh triangles and smooth mesh transition from the spot to the outer region. This run could only achieve the steady-state drop shape with reasonable accuracy, but crashed shortly thereafter due to drop–wall overlap, still far from the steady state for the drop velocity.

Figure 4. A simulation with $\lambda =30$, $\theta =7.5^\circ$, $B=1$ and $N_\triangle =8640$ using $n_3$-based passive mesh stabilization. (a) Side view along the $x_1$-axis. (b) Diagonal view from the drop bottom showing the dimple and the near-contact mesh adaptivity (with $r_{max}\approx 3)$. Only global steady-state drop shape (not the drop speed) could be approached in this low-resolution run.

3.3.3. Large $B$: non-adaptive meshing

Surprisingly, for strongly pancaked drops, the best choice was non-adaptive meshing dynamically controlled by passive stabilization, even at small tilt angles when a large drop portion is in very close contact with the wall. A simpler form (Appendix B) of the kinetic mesh energy function is used, which only prevents the mesh edges from becoming irregular and maintains the quality of mesh triangles. Figure 5(a,b) demonstrates the intended non-adaptive meshing in dynamical simulation with $\lambda =30$, $\theta =7.5^\circ$ and $B=5$; the mesh triangles are seen to be almost equilateral in the spot area (b) and on the rest of the drop (a). Again, a small number $N_\triangle =8640$ of triangular elements was chosen just to facilitate mesh visualization and quickly approach, with reasonable accuracy, the steady-state drop shape for these $\lambda$, $\theta$ and $B$. The run crashed shortly thereafter with drop–wall overlap well before the steady-state drop velocity could be reached; the steady-state thin-film profile could not be simulated either with such resolution (in particular, dimple imperfection (bulging) can be seen in figure 5b). Successful simulations with $\lambda \gg 1$ and small tilt angles require much higher $N_\triangle$ (especially, for strongly pancaked drops $B\gg 1$), which was made possible by multipole acceleration described below.

Figure 5. A simulation with $\lambda =30$, $\theta =7.5^\circ$, $B=5$ and $N_\triangle =8640$ using non-adaptive passive mesh stabilization. (a) Side view along the $x_1$-axis. (b) Diagonal view from the drop bottom showing the dimple and the quality of mesh triangles. Only global steady-state drop shape (not the drop speed) could be approached in this low-resolution run.

3.4. Multipole acceleration

Since the present simulations required long-time simulations with extreme resolutions (up to $N_\triangle \sim 3\times 10^5$ triangular elements), direct node-to-node summations for the BI (3.8), (3.10) or (3.11) were severely prohibited, and a highly efficient multipole-accelerated scheme was developed instead. Starting from the year 2000, multipole acceleration has been extensively developed and used in various multidrop–multiparticle BI simulations (see Zinchenko & Davis (2000), Zinchenko & Davis (2002), Zinchenko & Davis (2008), Zinchenko & Davis (2013), Zinchenko & Davis (2017) and other works by these authors). As discussed elsewhere (e.g. Zinchenko & Davis 2008), those techniques are distinct from the hydrodynamical version of the fast multipole method (FMM) (e.g. Sangani & Mo 1996) in many respects. In particular, instead of a hierarchy of space decompositions by Cartesian grids (characteristic of general FMM), a physically motivated surface partition into compact patches is used (with better behaviour of multipole expansions), combined with a broad use of rotational transformations for high-order expansions/re-expansions and economical truncation bounds. Compared with those previous works, the present, new multipole-accelerated scheme is most in line with the algorithm of Zinchenko & Davis (2005) developed for near-contact, free-space interaction of two slightly deformable drops with high resolution, but it is necessarily more complex due to complexities of the half-space Green function and related fundamental stresslet, and due to high-order near-singularity subtraction in the double-layer integral.

The discrete form of the BI equation (3.1) is obtained by approximating the non-singular integrals on the right-hand side of (3.8), (3.10) and (3.11) by the rule (3.15), using the area elements $\triangle S_j$ associated with nodes $\boldsymbol {x}^j$. Below, we briefly outline how the multipole-accelerated scheme works for the discretized double-layer integral in (3.1).

The mesh nodes are grouped into a large number $M$ of non-overlapping, compact patches (figure 6) ${\mathcal {B}}_1,{\mathcal {B}}_2,\ldots {\mathcal {B}}_M$, as detailed in Zinchenko & Davis (2005). This partition is made at $t=0$ only and while the nodes still form a uniform mesh (i.e. before projection on the initial drop surface in the operations of § 3.3 and Appendix B), and it is assumed unchanged in the subsequent drop motion. That way, each patch retains a constant set of ${\approx }N_\triangle /(2M)$ mesh nodes; typically, 200–400 nodes per patch were optimal for multipole acceleration. If mesh nodes adapt to the near-contact area, so do the patches.

Figure 6. A schematic for multipole-accelerated BI calculations using an extended system of mesh nodes (small dots), patches (bounded by hexagonal cells) and their mirror images in the lower half-space $x_3<0$. For patches ${\mathcal {B}}_\alpha$ and ${\mathcal {B}}_\beta$ well separated from ${\mathcal {B}}_\delta$, their ${\mathcal {B}}_\alpha \to {\mathcal {B}}_\delta$ and ${\mathcal {B}}_\beta \to {\mathcal {B}}_\delta$ contributions are calculated by singular-to-regular re-expansion of Lamb's series. The ${\mathcal {B}}_\gamma \to {\mathcal {B}}_\delta$ contribution is computed by pointwise Lamb's singular series for nodes $\boldsymbol {y}\in {\mathcal {B}}_\delta$ well outside the shell ${\mathcal {D}}_\gamma$, and by direct node-to-node summation otherwise.

Let ${\mathcal {B}}(\kern0.09em \boldsymbol {y})$ be the patch containing mesh node $\boldsymbol {y}$. The contribution of nodes $\boldsymbol {x}^j\in {\mathcal {B}}(\kern0.09em \boldsymbol {y})$ to the right-hand side of (3.10) is calculated by direct summation

(3.16)\begin{equation} \boldsymbol{\varPsi}_{SI}(\kern0.09em \boldsymbol{y})=\frac{3}{4{\rm \pi}}\sum_{\boldsymbol{x}^j\in {\mathcal{B}}(\kern0.09em \boldsymbol{y})} \frac{\{[\boldsymbol{Q}(\boldsymbol{x}^j)-\boldsymbol{Q}(\kern0.09em \boldsymbol{y})]\boldsymbol{\cdot}\boldsymbol{r}_j\} [\boldsymbol{r}_j\boldsymbol{\cdot}\boldsymbol{W}(\boldsymbol{x}^j)]\boldsymbol{r}_j} {r_j^5},\quad \boldsymbol{W}(\boldsymbol{x}^j)=\boldsymbol{n}(\boldsymbol{x}^j)\triangle S_j \end{equation}

(where $\boldsymbol {r}_j=\boldsymbol {x}^j-\boldsymbol {y}$, and $\boldsymbol {x}^j=\boldsymbol {y}$ is excluded from summation). Handling of other FS contributions to (3.10) is assisted by the subtraction matrices

(3.17)\begin{equation} \boldsymbol{\varPi}_1(\kern0.09em \boldsymbol{y})={-}\frac{3}{4{\rm \pi}}\sum_{\boldsymbol{x}^j\in S\setminus {\mathcal{B}}(\kern0.09em \boldsymbol{y})} \frac{[\boldsymbol{r}_j\boldsymbol{\cdot}\boldsymbol{W}(\boldsymbol{x}^j)]\boldsymbol{r}_j\boldsymbol{r}_j}{r_j^5}, \end{equation}

independent of $\boldsymbol {Q}$ and calculated for all $\boldsymbol {y}\in S$ before the BI iterations (see below). Additional matrices

(3.18)\begin{equation} \boldsymbol{\varPi}_2(\kern0.09em \boldsymbol{y}) ={-}\sum_{\boldsymbol{x}^j\in S} \boldsymbol{\tau}^{C}(\boldsymbol{x}^j;\boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{W}(\boldsymbol{x}^j), \end{equation}

precalculated only for $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|< h_o$, help with handling the discretized wall-correction contribution (3.11). Regarding other parts of (3.11), we note that the best local linear fit ${\boldsymbol{\mathcal L}}(\boldsymbol {x}-\boldsymbol {x}^\ast )$ to $\Delta \boldsymbol {Q}(\boldsymbol {x})= \boldsymbol {Q}(\boldsymbol {x}) - \boldsymbol {Q}(\boldsymbol {x}^\ast )$ can be written (Zinchenko & Davis 2006) as

(3.19)\begin{equation} {\boldsymbol{\mathcal L}}(\boldsymbol{x}-\boldsymbol{x}^\ast)= \sum_{\boldsymbol{x}^j\in{\mathcal{A}}(\boldsymbol{x}^\ast)} (a_j x'_1 + b_j x'_2)\Delta \boldsymbol{Q}(\boldsymbol{x}^j). \end{equation}

Here, ${\mathcal {A}}(\boldsymbol {x}^\ast )$ is the set of mesh nodes directly connected to $\boldsymbol {x}^\ast$; coefficients $a_j$ and $b_j$ depend on the local surface geometry near $\boldsymbol {x}^\ast$ and the choice of the intrinsic coordinates $(x'_1, x'_2)$ in the tangential plane at $\boldsymbol {x}^\ast$ (see § 3.2.1). Let $\varGamma ^k_l(\kern0.09em \boldsymbol {y})$ and $\varGamma ^k_{l+3}(\kern0.09em \boldsymbol {y})$ be the total ${\mathcal {L}}$-related contributions to (3.11) when ${\boldsymbol{\mathcal L}}(\boldsymbol {x}-\boldsymbol {x}^\ast )$ has the special form $x'_1\boldsymbol {e}_l$ or $x'_2\boldsymbol {e}_l$, respectively. These 18 scalars are also precalculated for all necessary nodes $\boldsymbol {y}$ before the iterations.

The discretized form of the double layer (3.1) to compute at every iteration becomes

(3.20)\begin{align} & \int_S\boldsymbol{Q}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\tau}(\boldsymbol{x};\boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})\mbox{d}S_x \approx \frac{1}{2}\boldsymbol{Q}(\kern0.09em \boldsymbol{y}) + \boldsymbol{\varPsi}_{SI}(\kern0.09em \boldsymbol{y}) + \boldsymbol{Q}(\kern0.09em \boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{\varPi}_1(\kern0.09em \boldsymbol{y}) \nonumber\\ &\quad +\varTheta(\kern0.09em \boldsymbol{y})\boldsymbol{Q}(\boldsymbol{x}^\ast)\boldsymbol{\cdot}\boldsymbol{\varPi}_2(\kern0.09em \boldsymbol{y}) + \varTheta_2(\kern0.09em \boldsymbol{y})\sum_{\boldsymbol{x}^j\in{\mathcal{A}}(\boldsymbol{x}^\ast)} [a_j\boldsymbol{\varGamma}_l(\kern0.09em \boldsymbol{y}) + b_j\boldsymbol{\varGamma}_{l+3}(\kern0.09em \boldsymbol{y})] [Q_l(\boldsymbol{x}^j)-Q_l(\boldsymbol{x}^\ast)] \nonumber\\ &\quad + \sum_{\boldsymbol{x}^j\in S\setminus {\mathcal{B}}(\kern0.09em \boldsymbol{y})} \boldsymbol{Q}(\boldsymbol{x}^j) \boldsymbol{\cdot} \boldsymbol{\tau}^{FS}(\boldsymbol{x}^j;\boldsymbol{y})\boldsymbol{\cdot} \boldsymbol{W}(\boldsymbol{x^j})+ \sum_{\boldsymbol{x}^j\in S} \boldsymbol{Q}(\boldsymbol{x}^j) \boldsymbol{\cdot} \boldsymbol{\tau}^{C}(\boldsymbol{x}^j;\boldsymbol{y})\boldsymbol{\cdot} \boldsymbol{W}(\boldsymbol{x^j}), \end{align}

with $\boldsymbol {\varGamma }_l=(\varGamma ^1_l,\varGamma ^2_l,\varGamma ^3_l)$; the barrier $\varTheta _2(\kern0.09em \boldsymbol {y})=1$ for $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|< h_1$ and zero otherwise.

The form (3.20) unloads the iterative part of the algorithm, but at the non-trivial expense of precalculating matrices $\boldsymbol {\varPi _1}$, $\boldsymbol {\varPi _2}$ and $\boldsymbol {\varGamma }$. Before addressing this issue, it is instructive to consider the multipole-accelerated scheme for the last two sums in (3.20). To invoke logical analogy with the algorithm of Zinchenko & Davis (2005) for two drops in free space, it is useful to view the last sum as the contribution from the mirror images $(\boldsymbol {x}^j)^{IM}$ of nodes $\boldsymbol {x}^j$. Accordingly, the node and patch system are extended to the lower half-space: for every node $\boldsymbol {x}^j$ and patch ${\mathcal {B}}_\gamma$ in the upper half-space there are the corresponding node $\boldsymbol {x}^{j+M}= (\boldsymbol {x}^j)^{IM}$ and patch ${\mathcal {B}}_{\gamma +M}={\mathcal {B}}_{\gamma }^{IM}$ in the lower half-space. Assuming this extension, a minimal spherical shell ${\mathcal {D}}_\gamma$ with centre $\boldsymbol {x}_\gamma ^o$ and radius $d_\gamma ^o$ is constructed (with sufficient accuracy) around each patch ${\mathcal {B}}_\gamma$. The total contribution to (3.20) from any patch ${\mathcal {B}}_\gamma$ ($\gamma \leq 2M$) is a Stokes flow as a function of $\boldsymbol {y}$, regular outside the shell ${\mathcal {D}}_\gamma$ and expandable there into Lamb's singular series about the shell centre $\boldsymbol {x}_\gamma ^o$. So, one can write either of the two expansions for a patch ${\mathcal {B}}$ lying in the upper half-space

(3.21)\begin{align} & \sum_{\boldsymbol{x}^j\in {\mathcal{B}}} \boldsymbol{Q}(\boldsymbol{x}^j) \boldsymbol{\cdot}\left\{\begin{array}{c} \boldsymbol{\tau}^{FS}(\boldsymbol{x}^j;\boldsymbol{y})\\ \boldsymbol{\tau}^{C}(\boldsymbol{x}^j;\boldsymbol{y}) \end{array}\right\}\boldsymbol{\cdot}\boldsymbol{W}(\boldsymbol{x}^j) =\sum_{\nu=0}^{\infty} \left[\boldsymbol{\nabla}\boldsymbol{\times}(\chi_{-(\nu+1)}\boldsymbol{Y}_{\gamma}) +\boldsymbol{\nabla}\varPhi_{-(\nu+1)} \vphantom{\frac{(\nu-2)Y_{\gamma}^2\boldsymbol{\nabla}p_{-(\nu+1)}} {2\nu(2\nu-1)}}\right.\nonumber\\ &\quad + \left.\frac{(\nu+1)p_{-(\nu+1)}\boldsymbol{Y}_{\gamma}}{\nu(2\nu-1)} - \frac{(\nu-2)Y_{\gamma}^2\boldsymbol{\nabla}p_{-(\nu+1)}} {2\nu(2\nu-1)}\right] \end{align}

(omitting the terms with $p_{-(\nu +1)}$ and $\chi _{-(\nu +1)}$ when $\nu =0$). Here, the expansion vector $\boldsymbol {Y}_{\gamma }=\boldsymbol {y}-\boldsymbol {x}_\gamma ^o$ and the negative-order solid harmonics are

(3.22)\begin{equation} \left.\begin{gathered} p_{-(\nu+1)}(\boldsymbol{Y}_{\gamma})= \sum_{\mu={-}\nu}^\nu A_{-(\nu+1),\mu}^{\gamma} \left(\frac{d_\gamma^o}{Y_{\gamma}}\right)^{\nu+1} {\mathcal{Y}}_{\nu,\mu}(\boldsymbol{Y}_{\gamma}), \\ \varPhi_{-(\nu+1)}(\boldsymbol{Y}_{\gamma})= \sum_{\mu={-}\nu}^\nu B_{-(\nu+1),\mu}^{\,\gamma} \left(\frac{d_\gamma^o}{Y_{\gamma}}\right)^{\nu+1} {\mathcal{Y}}_{\nu,\mu}(\boldsymbol{Y}_{\gamma}), \\ \chi_{-(\nu+1)}(\boldsymbol{Y}_{\gamma})= \sum_{\mu={-}\nu}^\nu C_{-(\nu+1),\mu}^{\gamma} \left(\frac{d_\gamma^o}{Y_{\gamma}}\right)^{\nu+1} {\mathcal{Y}}_{\nu,\mu}(\boldsymbol{Y}_{\gamma}), \end{gathered}\right\} \end{equation}

with complex coefficients $A_{-(\nu +1),\mu }^{\gamma }$, $B_{-(\nu +1),\mu }^{\gamma }$ and $C_{-(\nu +1),\mu }^{\gamma }$; ${\mathcal {Y}}_{\nu,\mu }(\boldsymbol {r})$ is the standard, normalized spherical harmonic. For the choice $\boldsymbol {\tau }^{FS}$ in (3.21), $\gamma \leq M$ is the one with ${\mathcal {B}}_\gamma ={\mathcal {B}}$; for $\boldsymbol {\tau }^{C}$ in (3.21), $\gamma >M$ corresponds to the mirror image of ${\mathcal {B}}$. For $\boldsymbol {\tau }^{FS}$, an efficient algorithm (Zinchenko & Davis 2008) was used to generate Lamb's $A$-, $B$- and $C$-coefficients in (3.22) to an arbitrary order; for $\boldsymbol {\tau }^{C}$, it was a novel, far more involved task undertaken in the present work (see Appendix C for details).

Based on Lamb's singular series generated for every patch ${\mathcal {B}}_\gamma$ to a sufficiently high order, the total of the last two sums in (3.20) is computed as follows. For every patch ${\mathcal {B}}_\gamma (\gamma \leq 2M)$ sufficiently separated from a chosen patch ${\mathcal {B}}_\delta (\delta \leq M)$ (i.e. with enough clearance between ${\mathcal {D}}_\gamma$ and ${\mathcal {D}}_\delta$), Lamb's singular series (3.21) is re-expanded near ${\mathcal {D}}_\delta$ into Lamb's regular series

(3.23)\begin{equation} \sum_{n=1}^{\infty}\left[\boldsymbol{\nabla}\boldsymbol{\times}(\chi_n\boldsymbol{Y}_{\delta}) +\boldsymbol{\nabla}\varPhi_n -\frac{np_n\boldsymbol{Y}_{\delta}}{(n+1)(2n+3)} + \frac{(n+3)Y_{\delta}^2\boldsymbol{\nabla}p_n}{2(n+1)(2n+3)}\right], \end{equation}

with positive-order solid harmonics $p_n(\boldsymbol {Y}_{\delta })$, $\varPhi _n(\boldsymbol {Y}_{\delta })$ and $\chi _n(\boldsymbol {Y}_{\delta })$. This re-expansion employs rotational transformations of spherical harmonics by Wigner functions (Zinchenko 1994; Zinchenko & Davis 2000), which makes this operation much less cumbersome to implement and reduces the operation cost from ${\sim }(\max (\nu,n))^4$ to ${\sim }(\max (\nu,n))^3$ compared with straightforward re-expansion without rotations. The cumulative regular series (3.23) from all the above patches ${\mathcal {B}}_\gamma$ is then calculated pointwise for all $\boldsymbol {y}\in {\mathcal {B}}_\delta$.

Regarding remaining patches ${\mathcal {B}}_\gamma$ (not sufficiently separated from ${\mathcal {B}}_\delta$), their contributions to the last two sums (3.20) for $\boldsymbol {y}\in {\mathcal {B}}_\delta$ are computed individually by either pointwise calculation of their Lamb's singular series (if $\boldsymbol {y}$ is sufficiently outside the shell ${\mathcal {D}}_\gamma$) or (in very rare cases) by direct node-to-node summation using the explicit forms (3.5) or (3.7).

Since the cumbersome elements $\boldsymbol {\varGamma }_l(\kern0.09em \boldsymbol {y})$ and $\boldsymbol {\varGamma }_{l+3}(\kern0.09em \boldsymbol {y})$ in (3.20) are required only for $\boldsymbol {y}$ very close to the wall, it was sufficient to precalculate them before the iterations in the simplest manner by direct node-to-node summation, with a small effect on the total computational cost. In contrast, tensor $\boldsymbol {\varPi }_2(\kern0.09em \boldsymbol {y})$ (active in a much broader range $\|\boldsymbol {y}^{IM}-\boldsymbol {x}^\ast \|< h_o$) and tensor $\boldsymbol {\varPi }_1(\kern0.09em \boldsymbol {y})$ (necessary on the entire surface $S$) had to be precalculated by multipole-accelerated schemes. Such schemes mimic the parts of the already described scheme for the last two sums in (3.20), just applied three times for constant vectors $\boldsymbol {Q}(\boldsymbol {x}^j)=\boldsymbol {e}_k$ with $k=1,2,3$.

Calculation of the inhomogeneous term (3.2) $\boldsymbol {F}(\kern0.09em \boldsymbol {y})$ also required multipole acceleration done along the same lines. Discretizing the non-singular right-hand side integrals (3.8) by the rule (3.15), their contributions are then split into parts coming (i) from $\tilde f(\boldsymbol {x})$ (in both integrals) (ii) from $\tilde f(\kern0.09em \boldsymbol {y})$ and (iii) from $\varTheta (\kern0.09em \boldsymbol {y})\tilde {f}(\boldsymbol {x}^\ast )$ in the brackets. Part (i) logically parallels the total contribution from the last two integrals in (3.20), while parts (ii) and (iii) are calculated in the same logical manner as tensors $\boldsymbol {\varPi }_1(\kern0.09em \boldsymbol {y})$ and $\boldsymbol {\varPi }_2(\kern0.09em \boldsymbol {y})$, respectively. All these tasks employ Lamb's singular series (3.21), this time for single-layer patch contributions generated to an arbitrary order, as outlined in Appendix C.

A crucial feature in all these multipole-accelerated operations is a broad spectrum of economical truncation bounds for multipole expansions/re-expansions. Based on plausible estimations of the convergence rate for Lamb's series, the algorithm decides (in a rational, nearly optimal way) how many terms to retain in singular-to-regular re-expansions (depending primarily on the clearance between the spherical shells around patches), and how to truncate singular Lamb's series in pointwise calculations (depending mostly on the clearance from an observation point to the spherical shell around a patch). The construction of truncation bounds controlled by a single, intuitive precision parameter $\varepsilon$ mostly follows that for two nearly touching drops in free space (Zinchenko & Davis 2005), with technical details in Zinchenko & Davis (2000). Here, again, the logical analogy between the drop mirror image in the present problem and the second drop in Zinchenko & Davis (2005) is used; a few adjustments specific to the present problem are listed in Appendix D. To enhance performance, the algorithm also has a threshold $k_o$ (set to 20) of the order of multipoles: if, for a given $\varepsilon$, particular multipole operations require higher orders, then direct node-to-node summations are invoked instead. There is some freedom in construction of truncation bounds, but the multipole-accelerated solution converges, as $\varepsilon \to 0$, to the same (much slower) non-multipole solution, regardless of $k_o$. This convergence is demonstrated in Appendix D, together with performance benchmarks.

4. Numerical analysis

4.1. Validation test: a spherical drop near a vertical wall

Although the present work is largely devoted to the most challenging case of small tilt angles $\theta$, the opposite case of a drop in very close proximity to a vertical wall ($\theta =90^\circ$) serves as a non-trivial check for correctness of the present multipole-accelerated BI code. In the Stokes regime and with infinite surface tension, a spherical drop, initially placed at an arbitrary clearance $\delta _o$ from the wall, would continue settling at a constant speed and without lateral migration; the drop–wall interaction in this case is conveniently characterized by the non-dimensional correction factor $\varDelta (\lambda,\delta _0)$ to the Hadamar–Rybchinsky drag force on an isolated drop. It is of interest to explore how the present solution can approach this limiting case and predict $\varDelta (\lambda,\delta _0)$. A deformable, initially spherical drop in this situation, however, would drift away from the wall for viscosity ratios $\lambda >3.12$ (Magnaudet, Takagi & Legendre 2003). Since this drift disappears in the limit $B\to 0$, it was legitimate in our tests for small $B\neq 0$ to exclude this drift altogether by forcing the drop surface centroid to stay at a distance $1+\delta _0$ from the wall. That way, the long-time, steady-state (non-dimensional) drop velocity $U_{st}$ along the wall could be unambiguously attained (figure 7) for selected $\lambda =10$ and 300. It is this velocity (not the initial velocity at $t=0$ when the drop was strictly spherical) and the related drag coefficient $\varDelta = 1/U_{st}$ that must be compared with those for a drop with infinite surface tension. The reason is that, in the BI formulation (3.2), the fully developed, small curvature variation $\Delta k$ (due to the normal stress balance on the interface) is divided by the Bond number to produce an $O(1)$ effect as $B\to 0$. In contrast, the infinite surface-tension formulation (see below) is only based on the tangential stress balance, and replaces the normal stress balance by postulating the spherical shape. In this respect, comparison with the case of a non-deformable drop is also a sensitive check on correctness of drop shape perturbation in our BI simulations. The difference between $U(t=0)$ and $U_{st}$ is modest for $\lambda =300$ (figure 7b), but is more significant (almost 1.4 times) for $\lambda =10$, $\delta _o=0.005$ and $B=0.007$ (figure 7a).

Figure 7. Transient speed of a slightly deformable drop settling along a vertical wall by the multipole-accelerated BI algorithm for (a) $\lambda =10$ and (b) $\lambda =300$. The drop starts from spherical, and its surface centroid is constrained to stay at a distance $1+\delta _0$ from the wall for the entire simulation.

Since small drop–wall clearances are of primary interest here, it would be attractive to use for comparison a highly accurate, $B=0$ bispherical (bipolar) coordinate solution for a spherical drop moving parallel to the wall, but such a solution could not be found in the literature (only the bubble case $\lambda =0$ for $\delta _o\geq 0.0453$ was handled, Meyyappan & Subramanian 1987). For this reason, the bispherical coordinate code of Zinchenko (1980) for two generally unequal spherical drops with arbitrary viscosities and instantaneous velocities normal to the line of centres had to be used at a size ratio of $10^{-6}$ and viscosity of the larger drop $10^7$ times larger than the medium viscosity to accurately represent the case of one spherical drop near a solid wall. However, to overcome ill conditioning for such extreme parameters, the code of Zinchenko (1980) had to be run in the quadruple precision mode.

In table 1, as $B\to 0$, our results by the multipole-accelerated BI code (with moderately adaptive, projective meshing and up to $N_\triangle =78\,\textrm {K}$ elements to guarantee the accuracy shown) are in excellent agreement with the $B=0$ bispherical coordinate solution; this limit is approached faster for smaller $\lambda$ and larger $\delta _o$. For $\lambda =300$ and $B=0$, our values of $\varDelta$ in table 1 are also consistent with the near-contact asymptotics

(4.1)\begin{equation} \varDelta(\infty, \delta_o)= \frac{8}{15}L+0.954 -\frac{4}{3} \frac{(L/10-0.193)^2} {(2L/5+0.371)},\quad L={-}\ln(\delta_o), \end{equation}

for a solid sphere moving with free rotation (which is the extension of the results from O'Neill & Stewartson 1967). A different validation (for small and large deformations) of the present multipole-accelerated BI code will be demonstrated in Appendix D, including code performance benchmarks.

Table 1. The drag force correction factor $\varDelta$ for a drop moving parallel to a plane wall; the $B\neq 0$ results are from the present multipole-accelerated, BI code, and the $B=0$ values are obtained by the bispherical coordinate code of Zinchenko (1980).

4.2. The effect of high-order double-layer desingularization

The benefits of full double-layer BI desingularization, achieved through (3.11), are most pronounced when it is applied to simulations with modest resolutions. In figure 8 for $\lambda =30$, $\theta =15^\circ$ and $B=0.125$, all four simulations employed projective meshing with the same pattern of near-contact adaptivity (as described in § 3.3.1 and § B.1) and started from a spherical drop shape with the initial drop–wall clearance $\delta _o=0.01$ or $0.02$. The curves 1 and 2 are for high resolutions of $N_\triangle =46$ and 61 K, respectively, with full desingularization. Reliable steady state was achieved in these two runs for both the drop–wall clearance $\delta _{min}(t)$ (figure 8a) and the transient drop velocity $U(t)$ (figure 8b); the latter is graphically indistinguishable between the two runs. Numerical convergence for $\delta _{min}$ up to steady state (with $\delta _{min}^{st}\approx 0.0061{-}0.0062$) is also good, considering how sensitive this small quantity is.

Figure 8. (a) Drop–wall clearance $\delta _{min}(t)$ and (b) transient drop speed in the simulations with $\lambda =30$, $\theta =15^\circ$ and $B=0.125$. Lines 1: $N_\triangle =46\,\textrm {K}$; 2: $N_\triangle =61\,\textrm {K}$; 3 and 4: $N_\triangle =8640$. The run for lines 3 was performed without high-order double-layer desingularization, leading to a crash at $t\approx 10$ with drop–wall overlap.

In contrast, the low-resolution run ($N_\triangle =8640$) with only leading-order subtraction in (3.11) (curves 3) crashed early due to drop–wall overlap (figure 8a), with subsequent divergence of BI iterations not allowing the simulation to proceed to steady state for the drop velocity (figure 8b). However, activating high-order desingularization for the same low resolution $N_\triangle =8640$ gives much better results (curves 4), with steady state achieved for the clearance (figure 8a) and drop velocity (figure 8b). In this run, roughly 12 % of the mesh nodes participated in the high-order desingularization (3.11) at large times. Although $\delta _{min}^{st}$ for this run is only half of the correct value, $U_{st}$ is in fair agreement with the correct result from the high-resolution simulations in figure 8(b). Note that such an improvement is achieved in a fundamental way through full BI desingularization, without attempts to impose an artificial repulsive barrier to prevent drop–wall overlap. Full desingularization is still beneficial for much higher resolutions, which is demonstrated in the next subsection.

4.3. Drop velocity: convergence analysis

Below, for the difficult combination $\lambda =60$, $\theta =7.5^\circ$, we explore the effects of discretization and solution strategies, primarily on the steady-state drop velocity $U_{st}$. It is far from obvious which of the three discretization strategies: (i) projective meshing of § 3.3.1, labelled as PR in what follows, (ii) $n_3$-based meshing (N3) of § 3.3.2 or (iii) non-adaptive meshing (NA) of § 3.3.3 to choose in each case, and the answer depends on the Bond number. In all the cases, acceptable results required very high surface resolutions only possible in dynamical simulations owing to multipole acceleration.

For $B=0.125$, with the drop exhibiting small deformation all the way to steady state, the drop speed trajectory $U(t)$ is presented in figure 9(a) for four simulations. Runs 1 and 2 employ projective meshing and started from a spherical drop shape with $\delta _o=0.01$. The steady-state meshing pattern for this $B$ is well represented by the low-resolution run in figure 3, but successful simulations of the steady-state drop speed and thin-film profile required much finer resolutions than $N_\triangle =8640$ could provide. Indeed, for runs 1 and 2 in figure 9(a) with $N_\triangle =61$ and 78 K, respectively, $U_{st}$ still differs by 2 %. To increase the resolution further, a one-time surface remeshing to $N_\triangle =138\,\textrm {K}$ was done at the end of run 2 preserving the mesh adaptivity pattern (by rebuilding the parametric mesh, as described in Appendix D), to provide the initial conditions for run 3 in figure 9(a), until a new steady-state drop velocity was reached. The difference between the $U_{st}$ values for 78 and 138 K resolutions is 1 %, indicating sufficient convergence (an additional run PR-246K done in a similar manner, but not included in figure 9(a), changed the steady-state drop speed from PR-138K by less than 0.17 %); the converged steady-state drop–wall clearance $\delta _{min}^{st}$ is 0.0034.

Figure 9. Transient drop speed for $\lambda =60$, $\theta =7.5^\circ$, with (a) $B=0.125$, (b) $B=1$, (c) $B=2$ and (d) $B=5$, using various discretization schemes (PR, N3, NA) combined with various resolutions $N_\triangle$. In (b), run 5 (dashed line) is the repeat of run 3, but without the default high-order double-layer desingularization. In each panel, the inset shows the steady-state drop–wall configuration (side view along $x_1$) from the highest-resolution run.

To test the alternative, $n_3$-based discretization scheme of § 3.3.2, the run N3-61K (line 4) was performed with $\alpha$ therein set to 0.84, to have approximately the same near-contact adaptivity as for the other runs in figure 9(a) (judging by the same maximum-to-minimum mesh edge ratio of four at this $B$). However, the steady-state drop speed prediction by N3-61K is only half as accurate as by PR-61K. Presumably, superiority of PR discretization is due to the higher accuracy of surface integration (after desingularization) though a uniform parametric mesh on the auxiliary unit sphere, while in the N3 scheme such integration has to be done on a strongly non-uniform native mesh on the drop surface. Additional tests confirmed that projective meshing is the method of choice for all Bond numbers $B\leq 0.5$ at small tilts $\theta$.

For higher $B\ge 1$, however, projective meshing is not a robust approach, and the choice was only between the $n_3$-based (with default $\alpha =0.7$) and non-adaptive surface discretizations. At $B=1$ and this $\theta$, the steady-state adaptivity pattern for $n_3$-based meshing is well represented by the lower-resolution run in figure 4. For this $B$, an excellent agreement is observed in figure 9(b) between the runs N3-61K and N3-78K all the way to steady state, with only 0.75 % difference in the steady-state drop velocities (the additional runs N3-138K and N3-246K, not included in figure 9(b), gave indistinguishable steady-state drop speed, to within 0.1 %). The convergent steady-state drop–wall clearance is 0.0066. Compared with the case $B=0.125$ (figure 9a), the numerical convergence for $B=1$ is faster, but it still requires high surface resolutions only feasible with multipole acceleration. For the N3-46K run (line 3 in figure 9b), the steady-state drop speed is only accurate to 6 %. Although the lubrication effects are now spread on a wider near-contact area than for $B=0.125$, near-contact mesh adaptivity is still essential for accurate results. For example, run 4 in figure 9(b) with non-adaptive mesh and $N_\triangle =46$ K is not nearly as good as N3-46K. Finally, the dashed line 5 in figure 9(b) demonstrates, again, the significance of full double-layer desingularization (3.11) for small tilt angles $\theta$ even in high-resolution simulations. This run was a repeat of N3-46K, only the default full desingularization was replaced by the simple leading-order subtraction (i.e. $h_1$ set to zero); the result was a crash with drop–wall overlap and quite incorrect drop speed not reaching the steady state. Note also that the pronounced shoulder of the curves in figure 9(b) (even upon numerical convergence) is due to the artificial initial condition (spherical shape at the small gap of 0.01) and drop deformability, making the surface clearance $\delta _{min}(t)$ non-monotonic, so that it grows in a small time interval after the initial decrease. Once $\delta _{min}(t)$ resumes its regular decline, the shoulder of the curves for $U(t)$ disappears.

For $B=2$ in figure 9(c), only the simulations N3-78K and N3-138K show good convergence, with $U_{st}$ differing by 1.2 % between the two runs (further increase of $U_{st}$ in an additional run N3-246K was only 0.5 %); the steady-state drop–wall gap is 0.0064. The run N3-61K underestimates the steady-state drop speed by approximately 5 %. The cruder resolution $N_\triangle =46$ K could not give satisfactory results at all, especially with non-adaptive meshing (line 5 for NA-46K), and so $n_3$-based meshing is still a preferable discretization strategy for this $B=2$.

It was particularly difficult to obtain the correct steady-state drop speed for a strongly pancaked drop, $B=5$ (figure 9d), with imperative ultrahigh surface resolutions; also, the partially deflated form of the BI equation (3.1) (instead of fully deflated) was mostly beneficial in this case to reach numerical convergence. The steady-state pattern of non-adaptive meshing for this $B$ is well represented by the lower-resolution run in figure 5. In figure 9(d), runs NA-246K and NA-328K both reliably reached the steady state, with excellent convergence all the way (the difference in $U_{st}$ between the two runs being only 0.5 %); the steady-state drop–wall clearance is 0.0046 for NA-246K and 0.0049 for NA-328K. However, run NA-138K (curve 3) could not reach the steady state due to discretization errors and had to be stopped (it would significantly underpredict $U_{st}$ anyway). Run NA-46K (line 4) had to be stopped even sooner due to poor accuracy. Since a strongly pancaked drop on a gently inclined wall is characterized by a large near-contact spot, one could expect a meshing adaptive to that region to perform better. However, the opposite is observed, and run N3-46K (line 5) is even less successful than NA-46K. Among possible explanations for this phenomenon, it can be noted that, for a given $N_\triangle$, mesh concentration in a large near-contact area significantly depletes the node density on the rest of the drop and makes it less uniform there, thus negatively affecting the accuracy of global surface integration. Such numerical integration errors are amplified at $\lambda \gg 1$, because of the proximity of $(\lambda -1)/(\lambda +1)$ to the spectrum of the BI equation (3.1). Also, for a strongly pancaked drop, the near-contact stresses are distributed on a larger area (than for small or moderate $B$), which makes them smaller, and it may be even less important to use very fine mesh in the spot area at the expense of the outer and transition regions. Thus, non-adaptive meshing with an ultrahigh number $N_\triangle$ of elements appears to be the only successful way in this case.

5. Physical results

Figure 10(a–d), the main result of this effort, presents the non-dimensional steady-state drop speed $U$ (as defined in § 2) for small-to-moderate tilt angles $7.5^\circ \leq \theta \leq 30^\circ$ in the wide ranges of the Bond number $0.0625\leq B \leq 5$ and viscosity ratios $1\leq \lambda \leq 300$. For each data point, the relative error does not exceed 0.5 % (and is often much smaller). To this end, super-resolutions with up to $N_\triangle =328$ K had to be used for the most difficult combinations of the small tilt $\theta =7.5^\circ$ with high viscosity ratios $\lambda \ge 60$. However, the case $\theta =7.5^\circ$, $\lambda =300$, $B=5$ still could not be studied, requiring even higher resolution to avoid non-convergence of BI iterations on the way to steady state. (This difficulty did not stem from the surface curvature/normal vector calculations by the default best paraboloid-spline method employed here for all $B\geq 1$. Using instead the high-order normal vector/curvature algorithm from Appendix B of Zinchenko & Davis (2006) for $\theta =7.5^\circ$, $\lambda =300$, $B=5$ and $N_\triangle =328$ K, the run stalled at almost the same time moment, with non-convergence of BI iterations. It remains to be seen in future work if improved schemes for surface integration could handle this case better.) On the other hand, for $\lambda =O(1)$ and moderate tilt angles $\theta \geq 22.5^\circ$, the resolution $N_\triangle =35$ K was always sufficient for high accuracy of $U$. For convenience, all the data from figure 10(a–d) are also presented in the Supplementary Material in tabular form available at https://doi.org/10.1017/jfm.2024.998.

Figure 10. Steady-state drop speed for tilt angles (a) $\theta =7.5^\circ$, (b) $\theta =15^\circ$, (c) $\theta =22.5^\circ$ and (d) $\theta =30^\circ$. In each panel, curves (top to bottom) correspond to viscosity ratios $\lambda =1,3,10,30,60$ and $300$, respectively.

The early efforts of Griggs et al. (2008, 2009) used a much simpler BI algorithm without multipole acceleration and high-order, double-layer near-singularity subtraction, and they focused on the range of tilt angles $\theta \geq 30^\circ$ and viscosity ratios $\lambda = O(1)$ to calculate the steady-state drop speed. While their resolution (up to $N_\triangle \sim 8640$) was quite acceptable in that generic range (especially in the simpler, most favourable case $\lambda =1$, see table 2), a few attempts therein to consider just slightly smaller tilts $\theta$ and/or slightly higher $\lambda$ quickly led to loss of accuracy. For example, at $\theta =15^\circ$, $B=3$ and $\lambda =1$, their solution underestimates $U$ by just 5 % (compared with the present code), but the error reaches $-$12 % for $\lambda =10$ at the same $\theta$ and $B$. In contrast, for $\lambda =20$ (the highest viscosity ratio in Griggs et al. 2008, 2009), $\theta =15^\circ$ and $B=4$, their solution overpredicted the drop speed by 15 %. A similar loss of accuracy is observed in the results of Griggs et al. (2009) for $\theta =7.5^\circ$, $\lambda =8.5$ and $0.25\leq B\leq 2$. Moreover, just slightly more extreme conditions could not be simulated at all with such algorithm and resolutions $N_\triangle \sim 8640$ due to non-convergence of BI iterations on the way to steady state. For example, simulations with $\lambda =60$, $\theta =30^\circ$ and $B=5$ failed; so did the runs with $\lambda =60$, $\theta =22.5^\circ$ and all $B\geq 2$. Simulations with $\lambda =300$ and $\theta =30^\circ$ could not succeed for any $B$ for the same reason (let alone more challenging, small tilt angles $\theta$). Thus, the present solution in figure 10(a–d) not only greatly improves the accuracy for generic values of $\theta$, $\lambda$ and $B$, but also greatly expands the parameter range owing to superhigh resolutions, feasible through multipole acceleration combined with the high-order near-singularity subtraction in the double-layer BI.

Table 2. Non-dimensional steady-state drop speed for $\theta =30^\circ$ and $\lambda =1$. The results from figure 5(c) of Griggs et al. (2008) are taken with the factor of two to account for our scale (2.2).

Inclusion of $\sin \theta$ in the velocity scale (2.2) considerably mitigates the dependence of $U$ on $\theta$ (for small $B$), but $U$ is still a sharp function of the Bond number at small tilt angles (figure 10a), especially with high viscosity ratios and near the spherical drop limit. For $\theta =7.5^\circ$, $U$ is a decreasing function of $B$ in the whole range $B\geq 0.0625$ for all $\lambda \geq 1$, which is confirmed by some simulations (with $\lambda =1$, 60 or 300) extended to $B=0.0625$ (figure 10a). At a larger tilt $\theta =15^\circ$ (figure 10b), the same behaviour is observed, except for $\lambda =1$ where the trend is just starting to reverse and $U(B)$ reaches a maximum at $B\approx 0.15$; high accuracy of our simulations (especially at $\lambda =1$) gives confidence that this trend reversal is not a numerical effect. For $\theta =22.5^\circ$ (figure 10c) and $\lambda =1$, this maximum of $U(B)$ becomes more noticeable; also, the trend reversal is seen to emerge at $\lambda =60$ while $U(B)$ remains monotonic to $B=0.0625$ for the higher viscosity ratio of 300. Finally, at the tilt angle of $30^\circ$ (figure 10d) and $\lambda =1$, the maximum of $U(B)$ is most pronounced, shifting to a higher Bond number $B\approx 0.5$; the $\lambda =60$ curve also shows a (very weak) maximum at $B\approx 0.1$, but the trend reversal is still not seen to $B=0.0625$ for $\lambda =300$. Note that Griggs et al. (2009) also observed weakly non-monotonic behaviour of $U(B)$ at $\theta =30^\circ$ and $\lambda \approx 1$, in both their experiments and simulations (table 2).

In their analytical study, Hodges et al. (2004) identified 11 asymptotic regimes for a non-wetting drop sedimenting on a gently inclined wall, to give leading-order estimates for the steady-state drop velocity and thin-film thickness depending on the relations between $\theta$, $\lambda$ and $B$. Their analysis naturally assumes logarithmically very broad variation of these parameters; the boundaries on different asymptotic regions are very tight, unless the tilt $\theta$ can be made extremely small (which was not possible in the present simulations). For strongly pancaked drops $B=5$ at $\theta =7.5^\circ$ and moderate viscosity ratios $\lambda =10{-}30$, an approximate agreement between our precise drop-speed results and the asymptotic theory predictions (region $\mbox {PII}_1$ in Hodges et al. 2004) was observed, although the difference could reach 1.5 times. With this exception, even though the present solution greatly expands the parameter space compared with Griggs et al. (2008, 2009), it was still hardly possible to make quantitative connections with the asymptotic theory of Hodges et al. (2004). For $B\ll 1$, one source of such difficulties is obvious from our figure 10. According to Hodges et al. (2004), as $B\to 0$ for fixed $\theta$ and $\lambda$, our $U(B)$ must approach zero, albeit slowly like ${\sim }|\ln (\theta ^2 B)|^{-1}$ (note the difference in the definitions of the non-dimensional drop speed $U$ between the present work and their work). Even at $\theta =30^\circ$ and $\lambda =1$, it would require extremely small $B$ (well past the reversal trend for $U(B)$ in figure 10d) to observe a substantial decrease of the drop speed, let alone smaller tilt angles and/or higher viscosity ratios. Such Bond numbers are computationally not feasible and would correspond to extremely small drop–wall clearance (see below), where surface roughness and other small-scale effects can come into play in practice.

Qualitatively, our trend reversal for $U(B)$ can be explained by the theory of Hodges et al. (2004). Namely, before reaching zero as $B\to 0$, $U(B)$ must first go through either $\mbox {FIII}_1$ or $\mbox {FIII}_2$ asymptotic regions of the map in their figure 12, with the scalings $U(B)\sim B^{-1/4}\theta ^{1/2}$ in both regions. Quantitatively, however, the asymptotics of Hodges et al. (2004) for these regions do not describe the well-tested results in our figure 10 at small $B$. One reason may be the neglect (in theory) of rim-edge energy dissipation, which scales (for $B\ll 1$) like $1.4B^{-1/5}(Ca)^{2/15}$ relative to the total dissipation elsewhere in the rim (Hodges et al. 2004); here, $Ca=\mu _e U_{dim}/\sigma$ is the capillary number based on the dimensional drop speed $U_{dim}$; this ratio must be small for the asymptotic theory to be applicable. However, based on our drop-speed results (figure 10), this ratio is approximately unity for all small $B$. (By a similar analysis, pancaked drops $B\sim 5$ are only slightly more favourable for the neglect of rim-edge dissipation.)

It was possible, however, to make comparisons with the experimental work of Rahman & Waghmare (2018) who measured the steady-state drop speed on a tilted glass plane submerged in a different liquid with viscosity ratio $\lambda =0.1$, $10$ or $110$ at very small Bond numbers. No sticking threshold was reported, and so their data presumably do not suffer from the contact angle hysteresis and correspond to the present model of a non-wetting drop. For $\lambda \geq 10$, $\theta \leq 15^\circ$ and $B\leq 0.04$, the dimensional drop speed $U_{dim}$ was fitted in Rahman & Waghmare (2018) by a semi-empirical relation

(5.1)\begin{equation} U_{dim}\sim \frac{B\sigma\sin\theta}{\mu_d B^{3/2}+6{\rm \pi}\mu_e}\end{equation}

(although it is difficult to judge how accurately this fit represents their data since their dimensional speed covers a very broad range represented logarithmically). Using our velocity scale (2.2) to make (5.1) non-dimensional, comparison was made with our precise results in figure 10 for all pairs with $\theta \leq 15^\circ$, $\lambda \geq 10$ and $B\leq 0.125$ (table 3). In most cases, there is only a modest difference, although, in a few instances (with $\theta =7.5^\circ$ and large $\lambda \geq 30$), (5.1) overpredicts our $U$ by 38 %–42 %. This equation does not capture the non-monotonic behaviour of our $U(B)$ in figure 10 and quickly becomes inadequate with the increase in the Bond number. Again, the lack of simulation data for $B\ll 0.0625$ did not allow for more representative comparisons with the experiments of Rahman & Waghmare (2018).

Table 3. Comparison of the non-dimensional steady-state drop speed $U$ from the present simulations with the semi-empirical formula of Rahman & Waghmare (2018). For each set of $\theta$, $B$ and $\lambda$, the upper value is from our figure 10; the lower value is from (5.1) scaled on (2.2).

Another difficulty with connecting our work to Hodges et al. (2004) is that their 3-D thin-film thickness is characterized by a single parameter $\epsilon$, while we pursue a more complete description of the dimple observed in the present simulations (as well as in those of Griggs et al. 2008, 2009); this dimple very slowly disappears as $B\to 0$. Instead of the topographic maps (Griggs et al. 2008, 2009) for the film profile, we use here a more compact dimple characterization in steady state by three parameters, $\delta _{min}$, $\delta _{av}$ and $\delta _{max}$. The first one (already used in § 4.2) is the global minimum of the mesh node-to-wall distances, while $\delta _{av}$ and $\delta _{max}$ are, respectively, the surface average and maximum of such distances calculated in the dimple region only. This region is unambiguously defined as the set of drop-mesh nodes where at least one of the principal surface curvatures $k_1$ or $k_2$ is negative. Obviously, such definition makes sense for small tilt angles, when the drop is in near contact with the wall; note that $\delta _{min}$ is reached very slightly outside the rim of the so-defined dimple.

These three metrics, $\delta _{min}$, $\delta _{av}$ and $\delta _{max}$ (all well convergent with respect to resolutions) are presented in figure 11 for $\theta =7.5^\circ$ (a–c) and $15^\circ$ (d–f). In each panel, red squares are for $\lambda =1$, green circles for $\lambda =60$ and black diamonds for $\lambda =300$. For small $B$, each quantity in figure 11(a–c) approximately obeys the scaling ${\sim }B^{1/2}\sin \theta$, in qualitative agreement with the asymptotic analysis of Hodges et al. (2004), although this scaling is likely not the ultimate one for $B\to 0$, since their theoretical estimate also includes a logarithmically weak factor of $O(|\ln (\theta ^2 B)|^{-1})$ (without the after-log term). However, the above scaling fits our $B\ll 1$ results in figure 11(a–c) with numerical factors almost independent of $\lambda =1{-}300$, but sensitive to the choice of the metric (namely, the factor is within 0.07–0.076 for $\delta _{min}$, and 0.13–0.14 for $\delta _{max}$). The asymptotic theory does not discern such thin-film details sufficiently for comparison with our simulations (in particular, $\delta _{min}$ is not resolved in the theory).

Figure 11. Steady-state thin-film metrics (a,d) $\delta _{min}$, (b,e) $\delta _{av}$ and (c, f) $\delta _{max}$ vs Bond number for (a–c) $\theta =7.5^\circ$ and (d–f) $\theta =15^\circ$. Red squares: $\lambda =1$, green circles: $\lambda =60$, black diamonds: $\lambda =300$.

For a pancaked drop $B\geq O(1)$, the geometry of the lubrication space becomes strongly sensitive to the viscosity ratio. While the film thickness saturates with the increase in the Bond number to $B\approx 1{-}2$ for a highly viscous drop $\lambda =60$, it continues to grow for $\lambda =1$, at least in the range $B\leq 5$, thus making the film much thicker for a homoviscous drop compared with a drop with high viscosity ratio. This observation is confirmed for all the metrics in figure 11(a–c).

For a larger tilt $\theta =15^\circ$ and high viscosity ratios, only $\delta _{min}$ saturates with the increase in the Bond number (figure 11d), while the two other thin-film metrics continue to grow monotonically with $B$ (figure 11e, f); the film is still thicker for a homoviscous drop than for a drop with $\lambda \geq 60$. Comparison between (a–c) and (d–f) shows that a linear scaling of the film thickness parameters with the tilt angle $\theta$ or $\sin \theta$ roughly holds for $\delta _{av}$ and $\delta _{max}$ in case of highly viscous drops, and for all three metrics in case of homoviscous drops. However, significant deviations from this scaling are observed for $\delta _{min}$ in case of highly viscous pancaked drops.

The asymptotic lubrication analysis of Hodges et al. (2004) for $\theta \ll 1$ predicts, to the leading order, uniform film thickness $\delta$ from the back to the front of the film, i.e. in the notations of our figure 1, $\delta$ must be only a function of $x_1$. To verify this prediction, $\delta$ was computed for $\theta =7.5^\circ$ on the central cross-section of the dimple region (i.e. in the $(x_2, x_3)$-plane drawn through the drop centroid $\boldsymbol {x}^c$) and plotted vs $x_2 -x_2^c$. (Some lack of smoothness in figure 12 may be due to presentation, with the simplest, zero-order interpolation from the 3-D mesh of triangles to the central plane.) For a pancaked drop $B=2$ (figure 12a), $\delta$ is indeed practically constant in the dimple cross-section, except for the narrow range of $x_2-x_2^c$ at the tail where the film is substantially thinner. This behaviour is observed for $\lambda =1$, 60 and especially for 300. Note that the global minimum of $\delta$ for each of these $\lambda$ is not approached in figure 12(a); it is significantly smaller (see figures 11(a) and 13(a)) and reached well outside the central plane $(x_2, x_3)$, so that $\delta$ remains a strong function of $x_1$. For a slightly deformable drop $B=0.0625$ (figure 12b), however, even film uniformity along the drop motion direction cannot be assumed as an approximation; here, a large portion of the film cross-section (near the tail) is much thinner than the rest of the film in the $(x_2, x_3)$-plane and does not approach a constant thickness as $x_2$ grows.

Figure 12. Profile of the steady-state thin-film thickness in the central cross-section of the dimple region for $\theta =7.5^\circ$. Panels show (a) $B=2$; (b) $B=0.0625$. Lines 1: $\lambda =1$, 2: $\lambda =60$, 3: $\lambda =300$.

Figure 13. Contour plots of the steady-state thin-film thickness $\delta (x_1,x_2)$ in the near-contact spot for $\theta =7.5^\circ$ and $\lambda =300$. (a) $B=2$; (b) $B=0.0625$. View from the drop bottom. In (a), the outer boundary shows the entire drop projection on the $(x_1,x_2)$-plane. The minimum clearance $\delta _{min}$ is 0.0049 for (a) and 0.00248 for (b).

The lubrication space geometry in the entire near-contact spot is also much different between the $B=2$ and $B=0.0625$ cases, as demonstrated in figure 13 for $\lambda =300$. While the contour lines for drop–wall clearance in figure 13(a) tend to be parallel to the $x_2$-axis, there is no such region at $B=0.0625$ in figure 13(b). On the other hand, the distribution of $\delta$ in the trailing part of the spot resembles axisymmetric for $B=0.0625$ – a feature absent for $B=2$. These observations additionally explain particular difficulties of comparison with the asymptotic theories for $B\ll 1$, requiring much smaller $\theta$ and $B$ than possible in rigorous simulations.

Finally, all the results in figures 11 and 12 indicate that a non-uniform 3-D film profile with finite thickness is approached, albeit slowly, for every fixed $\theta$ and $B$, as $\lambda \to \infty$. This last observation is in line with the arguments of Hodges et al. (2004) that the steady-state, high viscosity ratio limit must correspond to pure drop sliding as a rigid body, since tumbling is prohibited by the drop non-sphericity. It can be added to their arguments that the lack of fore–aft symmetry of the drop shape developing in this limit explains how such mode of the drop motion is possible at all at zero Reynolds number. A solid sphere moving parallel to the wall would only experience a lateral hydrodynamic force, not capable to counter the gravity effect and maintain a steady-state film thickness. More generally, lateral translation of a non-spherical solid particle with fore–aft symmetry would not produce a normal hydrodynamic force either (as follows from the reciprocal theorem) necessary for the particle to glide, and so a loss of such symmetry is crucial to maintain a steady-state regime.

Besides sliding, Hodges et al. (2004) also classify other possible steady-state kinematics for a drop motion down an inclined wall. One of them is tank treading for a pancaked drop, where the circulation fluid velocity inside the drop is comparable in magnitude to the drop speed, and is uniform over the near-contact spot. The other regime is rolling of a nearly spherical drop, with rigid-body rotational motion comparable in magnitude to the drop speed. It is of interest to see if the present simulations for $\lambda \gg 1$ and small tilt angles fall, at least qualitatively, into any of these regimes.

To this end, the fluid circulation velocity $\boldsymbol {u}-\boldsymbol {U}$ was computed along the central cross-section of the whole drop surface (i.e. in the $(x_2, x_3)$-plane drawn through the drop centroid) and presented in figure 14(a–l) for different combinations of $\lambda$, $\theta$ and $B$. For each combination, a reference vector (placed inside the drop) represents the steady-state drop velocity $\boldsymbol {U}$. For each field vector, its length relative to the length of the reference vector gives the local circulation intensity ${\mathcal {C}}= \|\boldsymbol {u}-\boldsymbol {U}\|/U$ at the starting surface point of the field vector. The maximum and minimum values of ${\mathcal {C}}$ over the entire cross-section contour are also listed for each set of $\lambda$, $\theta$ and $B$; the maximum is reached in the near-contact spot close to the trailing edge, while the minimum is observed just in front of the spot. The global drop shape in figure 14 is a strong function of $B$, but is almost insensitive to $\lambda$ and $\theta$, since it is close (for small tilts) to the hydrostatic equilibrium shape on a horizontal plane unaffected at all by the viscosity ratio (e.g. Hodges et al. 2004).

Figure 14. Steady-state circulation velocity $\boldsymbol {u}-\boldsymbol {U}$ along the central cross-section of the drop surface. Length of each field vector on a contour relative to the length of the reference vector (inside this contour) gives local circulation intensity ${\mathcal {C}}= \|\boldsymbol {u}-\boldsymbol {U}\|/U$.

However, $\lambda$ and $\theta$ have a strong effect on the circulation intensity ${\mathcal {C}}$. For a pancake drop $B=2$ and tilt $\theta =15^\circ$, the sliding regime is quickly approached as $\lambda$ is increased from 60 to 300: in figure 14(b), the circulation intensity is within 5 % on the whole drop contour, which is an almost perfect sliding regime, while the $\lambda =60$ motion (a) is better described as tank treading (with large slip at the wall). Decreasing the Bond number to 0.5 (c,d) enhances the circulation, especially for $\lambda =300$, but the kinematics in (d) is still close to sliding, judging by the small values of ${\mathcal {C}}$. For an even smaller $B=0.0625$ in (e, f), only the vicinity of the near-contact spot is shown; the rest of the drop is nearly spherical, with nearly constant ${\mathcal {C}}\approx {\mathcal {C}}_{min}$. A marked difference between (e, f) is much more uniform circulation for $\lambda =300$, and so the kinematics in ( f) is best described as rolling with sliding; the circulation here does not exceed 27 %.

At the smaller tilt $\theta =7.5^\circ$ (g–l), convergence to the pure sliding regime $\lambda =\infty$ is much slower and would require much larger values of the viscosity ratio; in (h) for $\lambda =300$ and $B=2$, the circulation intensity still reaches 11 % (cf with (b)). Again, for the smallest $B=0.0625$, increasing $\lambda$ from 60 (k) to 300 results in a much more uniform circulation, so that the kinematics in (l) resembles rolling with sliding; the value of ${\mathcal {C}}_{max}=0.35$ is still away from unity. It appears that for liquid–liquid systems in a purely hydrodynamical formulation, with a smooth solid wall and in the absence of singular adhesive forces, the thin-film lubrication can never provide for perfect rolling (i.e. ${\mathcal {C}}\approx 1$).

It is not obvious from figure 10 how our results could be used to accurately predict the drop speed for an arbitrary viscosity ratio from 1 to 300, or in a wider range $\lambda >300$. In this respect, it was found helpful (mostly for moderate-to-large Bond numbers) to plot $U(\lambda +1)/(\lambda +2/3)$ vs $\lambda ^{-1/3}$ (figure 15). For $B=1$, $2$ and $5$, and all tilt angles from $7.5^\circ$ to $30^\circ$, a smooth dependence on $\lambda ^{-1/3}$ is observed, allowing for accurate interpolation to an arbitrary $\lambda \in (1,300)$. Moreover, at the small tilt of $7.5^\circ$, the results are almost linear with $\lambda ^{-1/3}$ to $\lambda =300$, allowing for reasonably credible predictions of the drop speed to $\lambda \approx 1000$ by extrapolation. Theoretical explanation of such, almost linear, behaviour in a wide range of $\lambda$ is lacking and it is, obviously, not the ultimate behaviour for $\lambda \to \infty$. Indeed, for fixed $B$ and $\theta$, the finite-$\lambda$ correction to the drop speed must be asymptotically $O(\lambda ^{-1})$ (since it would become a problem of regular perturbation). This change of behaviour is seen most vividly for the moderate tilt of $30^\circ$ (figure 15d), where the convergence to the sliding limit $\lambda =\infty$ is much faster for $B=5$ than for $B=1$ due to larger drop–wall separation in the former case. Again, a qualitative connection is seen with the work of Hodges et al. (2004): from their figure 12, the sliding regime requires $\lambda \gg (B^{1/2}\sin ^2 \theta )^{-1}$, i.e. smaller viscosity ratios as the Bond number and/or the tilt angle are increased.

Figure 15. Steady-state values of $U(\lambda +1)/(\lambda +2/3)$ for (a) $\theta =7.5^\circ$, (b) $\theta =15^\circ$, (c) $\theta =22.5^\circ$ and (d) $\theta =30^\circ$. Curves (top to bottom) are for $B=1$, $2$ and $5$, respectively.

It would be interesting, albeit challenging, to have special simulation or semi-analytical methods for the sliding limit $\lambda =\infty$. One difficulty is that the relevant solid-body shape with a dimple and wide near-contact spot is unknown a priori and must be a part of the solution. Also, the available asymptotic approach to handle close hydrodynamical interaction of a non-spherical solid particle with a wall or another particle (Cox 1974; Claeys & Brady 1989; Staben, Zinchenko & Davis 2006) is based on lubrication analysis at the single point of near contact and would not be applicable in the present case.

6. Conclusions

Gravity-driven, low-Reynolds-number, steady-state motion of a viscous, immiscible deformable drop on an inclined plane wall in another liquid has been studied from first principles, by rigorous 3-D BI simulations. In the absence of singular adhesive forces and wall surface roughness, such a drop is unable to make contact with the wall, and remains separated by a lubricating film due to drop motion. The focus has been on small-to-moderate wall inclination angles $\theta$ and high drop-to-medium viscosity ratios $\lambda$. In these conditions, the drop stays extremely close to the wall with strong hydrodynamical interaction, and the successful solution has met a much greater challenge than for the generic case (Griggs et al. 2008, 2009) of moderate $\lambda$ and $\theta$. The BI formulation is based on the half-space Green function and related fundamental stresslet to reduce the problem to a BI equation for the interfacial velocity on the drop surface only. The wall-correction contribution to this equation leads to near-singular integrands for a drop in near contact with the wall. A novel high-order near-singularity subtraction in the double-layer BI offers (for the first time) full desingularization of the BI equation for any drop–wall clearance, and greatly reduces the numerical trend for drop–wall overlap. Still, the proximity of $(\lambda -1)/(\lambda +1)$ (for $\lambda \gg 1$) to the spectrum of the BI equation necessitates extreme drop surface resolutions to avoid non-convergence of BI iterations on the way to steady state, eliminate overlaps altogether and obtain the convergent steady-state drop speed with respect to discretization. Such superhigh resolutions (with up to $3\times 10^5$ boundary elements) were made possible in dynamical simulations through highly efficient multipole acceleration; the techniques for such acceleration are far more involved in some aspects than in the case of the free-space Green function and stresslet. Fixed topology triangulations were used to discretize the drop surface and track its evolution to steady state with different adaptive/non-adaptive versions required, depending mostly on the Bond number range.

Using this novel algorithm, the non-dimensional steady-state drop speed $U$ (scaled on the Hadamar–Rybchinsky settling speed (2.2) in reduced gravity $g\sin \theta$) was studied systematically, and with high accuracy, for $\theta$ down to $7.5^\circ$, $\lambda$ up to 300 and the Bond number $0.0625\leq B\leq 5$ covering the range from nearly spherical to strongly pancaked drops. The results are consistent with the semi-empirical formula of Rahman & Waghmare (2018) roughly describing their experimental data for small $B$ and $\lambda =10{-}100$. At $\theta =7.5^\circ$, our $U$ is a monotonically decreasing function of $B$ for all $\lambda =1{-}300$; the dependence on $B$ is strong for high viscosity ratios. In contrast, at $\theta =30^\circ$, $U(B)$ becomes noticeably non-monotonic for $\lambda =1$, but not for high $\lambda$. Comparison with the leading-order asymptotic theory of Hodges et al. (2004) could not be made here, since their drop-speed behaviour $U(B)\to 0$ presumably requires much smaller Bond numbers than it is possible to simulate. Our results for the thickness of the dimpled lubrication zone ($\delta _{min}$, $\delta _{av}$ and $\delta _{max}$), however, all approximately obey the scaling ${\sim }B^{1/2}\sin \theta$ for small $B$, in qualitative agreement with the asymptotic analysis of Hodges et al. (2004) (although their additional, inverse logarithmic factor could not be verified). For small $B$, the above thin-film metrics are practically unaffected by $\lambda$ from 1 to 300. In contrast, for pancaked drops $B\geq O(1)$ and small inclination angles, the lubrication space geometry becomes strongly sensitive to the viscosity ratio. While the film thickness saturates with the increase in the Bond number to $B\approx 1{-}2$ for highly viscous drops, it continues to grow for $\lambda =1$, at least in the range $B\leq 5$, thus making the film much thicker than for $\lambda \gg 1$ drops.

At $\theta \ll 1$, the asymptotic analysis of Hodges et al. (2004) predicted almost uniform film thickness in the drop motion direction (but strong non-uniformity in the orthogonal direction). The present simulations confirm this prediction for a pancaked drop ($B=2$) at all viscosity ratios $\lambda =1{-}300$, but not for a nearly spherical drop ($B=0.0625$). For this reason, quantitatively, the regime $B\to 0$ of the asymptotic theory may require much smaller inclination angles than possible to simulate.

The steady-state fluid circulation velocity $\boldsymbol {u}-\boldsymbol {U}$ on the central cross-section of the drop surface was also monitored to observe changes in the flow kinematics. Increasing $\lambda$ from 60 to 300 for a pancaked drop $B=2$ leads to a sliding regime (with small circulation); this approach to sliding is much faster for the $15^\circ$ than for $7.5^\circ$ tilt angle. For a slightly deformed drop ($B=0.0625$) at $\lambda =300$ and $\theta =7.5^\circ$, the circulation is close to uniform, and the mode of motion is best described as rolling, but with a large slip ($|\boldsymbol {u}-\boldsymbol {U}|\ll U)$. Neither perfect tank treading, nor perfect rolling (with the drop fluid locally pinned to the wall in both cases) can be approached for liquid–liquid systems in the pure hydrodynamical formulation considered herein; those modes of motion could only result with the help from additional mechanisms (e.g. strong adhesive forces or wall surface roughness).

Finally, it was shown in this work how our drop-speed results obtained for a limited set of $\lambda \in [1,300]$ can be very accurately interpolated to an arbitrary viscosity ratio in that range, at least for $B\geq 1$; also, for small tilt angles, a reasonably credible extrapolation to $\lambda \approx 1000$ can be made. The present results can stimulate further experimental work, as well as the development of advanced, high-order asymptotic theories for better comparison with simulations.

The present methodology can be applied to a variety of other problems, e.g. motion of a surfactant-covered drop or a vesicle on an inclined wall. The ability to handle extremely small drop–wall clearances through multipole acceleration with ultrahigh surface resolutions opens a way to directly include short-range colloidal forces (Van der Waals attraction and electrostatic repulsion from the Derjaguin–Landau–Verwey–Overbeek theory) in 3-D BI simulations (although first-principle simulations with rough walls still present a formidable challenge). The full double-layer desingularization tool found in this work is even more universal, since it is also applicable to near-singular integrals with the free-space fundamental stresslet. Hence, numerous problems of near-contact drop–drop, or drop–particle interaction in free space or a periodic box with extreme viscosity ratio can be solved using this tool (thus breaking the tradition that such 3-D BI problems have been addressed so far for $\lambda =O(1)$ only).