Hostname: page-component-cd9895bd7-q99xh Total loading time: 0 Render date: 2024-12-28T16:18:36.662Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  26 November 2024

Alexander Z. Zinchenko*
Affiliation:
Department of Chemical and Biological Engineering, University of Colorado, Boulder, CO 80309-0424, USA
*
Email address for correspondence: [email protected]

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.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
© The Author(s), 2024. Published by Cambridge University Press.

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 Reference Dussan and Chow1983; Dussan Reference Dussan1985). 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é Reference Richard and 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é Reference Aussillous and 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 (Reference Mahadevan and Pomeau1999) 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 (Reference Mahadevan and Pomeau1999) was recently derived (Schnitzer, Davis & Yariv Reference Schnitzer, Davis and Yariv2020) 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 (Reference Mahadevan and Pomeau1999) and Schnitzer et al. (Reference Schnitzer, Davis and Yariv2020) 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é Reference Richard and Quéré1999; Aussillous & Quéré Reference Aussillous and Quéré2001; Aussillous Reference Aussillous2002) for a drop (a mixture of water and glycerol) on a super-hydrophobic wall in the air (see also Quéré (Reference 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. Reference Del Castillo, Ohnishi, White, Carnie and Horn2011), surface roughness (Quéré Reference 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 Reference Aussillous2002).

Hodges, Jensen & Rallison (Reference Hodges, Jensen and Rallison2004) 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é (Reference Richard and 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. (Reference Hodges, Jensen and Rallison2004), 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. (Reference Hodges, Jensen and Rallison2004), this behaviour is in stark contrast to the experiments of Richard & Quéré (Reference Richard and 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é (Reference Richard and 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. (Reference Hodges, Jensen and Rallison2004) 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é Reference Aussillous and Quéré2002; Griggs, Zinchenko & Davis Reference Griggs, Zinchenko and Davis2009; Rahman & Waghmare Reference Rahman and Waghmare2018). 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 (Reference Griggs, Zinchenko and Davis2008) and Griggs et al. (Reference Griggs, Zinchenko and Davis2009) 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 (Reference Blake1971) 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009), 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009) (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. (Reference Hodges, Jensen and Rallison2004) 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009) 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009), 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. (Reference Hodges, Jensen and Rallison2004); 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009), 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 (Reference Blake1971) and reviewed in Pozrikidis (Reference Pozrikidis1992), 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 Reference Rallison and Acrivos1978; Pozrikidis Reference Pozrikidis1992), 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 (Reference Zinchenko, Rother and Davis1997). 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 Reference Kim and Karilla1991; Pozrikidis Reference Pozrikidis1992) 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 (Reference Blake1971) and Pozrikidis (Reference Pozrikidis1992) (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. (Reference Griggs, Zinchenko and Davis2009) (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 Reference Bazhlekov, Anderson and Meijer2004; Klaseboer, Sun & Chan Reference Klaseboer, Sun and Chan2012; Farutin, Biben & Misbah Reference Farutin, Biben and Misbah2014) could fully eliminate near-singular behaviour of the FS integrals (3.10), when $\boldsymbol {y}$ is off the integration surface. Zinchenko & Davis (Reference Zinchenko and Davis2002) 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 (Reference Zinchenko and Davis2005) 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 Reference Zinchenko and Davis2006, Reference Zinchenko and Davis2008), $\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 (Reference Gissinger, Zinchenko and Davis2021) 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 (Reference Zinchenko and Davis2006). 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 (Reference Zinchenko and Davis2005) 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. Reference Zinchenko, Rother and Davis1997). 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 (Reference Zinchenko and Davis2005), 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 (Reference Zinchenko and Davis2013), 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 (Reference Zinchenko and Davis2013) 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 (Reference Zinchenko and Davis2000), Zinchenko & Davis (Reference Zinchenko and Davis2002), Zinchenko & Davis (Reference Zinchenko and Davis2008), Zinchenko & Davis (Reference Zinchenko and Davis2013), Zinchenko & Davis (Reference Zinchenko and Davis2017) and other works by these authors). As discussed elsewhere (e.g. Zinchenko & Davis Reference Zinchenko and Davis2008), those techniques are distinct from the hydrodynamical version of the fast multipole method (FMM) (e.g. Sangani & Mo Reference Sangani and Mo1996) 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 (Reference Zinchenko and Davis2005) 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 (Reference Zinchenko and Davis2005). 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 Reference Zinchenko and Davis2006) 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 (Reference Zinchenko and Davis2005) 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 Reference Zinchenko and Davis2008) 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 Reference Zinchenko1994; Zinchenko & Davis Reference Zinchenko and Davis2000), 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 Reference Zinchenko and Davis2005), with technical details in Zinchenko & Davis (Reference Zinchenko and Davis2000). Here, again, the logical analogy between the drop mirror image in the present problem and the second drop in Zinchenko & Davis (Reference Zinchenko and Davis2005) 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 Reference Magnaudet, Takagi and Legendre2003). 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 Reference Meyyappan and Subramanian1987). For this reason, the bispherical coordinate code of Zinchenko (Reference Zinchenko1980) 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 (Reference Zinchenko1980) 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 Reference O'Neill and Stewartson1967). 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 (Reference Zinchenko1980).

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(ad), 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 (Reference Zinchenko and Davis2006) 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(ad) 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009) 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. Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009), $\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. (Reference Griggs, Zinchenko and Davis2009) 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(ad) 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. (Reference Griggs, Zinchenko and Davis2008) 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. (Reference Griggs, Zinchenko and Davis2009) 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. (Reference Hodges, Jensen and Rallison2004) 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. Reference Hodges, Jensen and Rallison2004) 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. (Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009), it was still hardly possible to make quantitative connections with the asymptotic theory of Hodges et al. (Reference Hodges, Jensen and Rallison2004). For $B\ll 1$, one source of such difficulties is obvious from our figure 10. According to Hodges et al. (Reference Hodges, Jensen and Rallison2004), 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. (Reference Hodges, Jensen and Rallison2004). 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. (Reference Hodges, Jensen and Rallison2004) 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. Reference Hodges, Jensen and Rallison2004); 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 (Reference Rahman and Waghmare2018) 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 (Reference Rahman and Waghmare2018) 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 (Reference Rahman and Waghmare2018).

Table 3. Comparison of the non-dimensional steady-state drop speed $U$ from the present simulations with the semi-empirical formula of Rahman & Waghmare (Reference Rahman and Waghmare2018). 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. (Reference Hodges, Jensen and Rallison2004) 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. Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009); this dimple very slowly disappears as $B\to 0$. Instead of the topographic maps (Griggs et al. Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009) 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$ (ac) and $15^\circ$ (df). 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(ac) approximately obeys the scaling ${\sim }B^{1/2}\sin \theta$, in qualitative agreement with the asymptotic analysis of Hodges et al. (Reference Hodges, Jensen and Rallison2004), 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(ac) 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 (cf) $\delta _{max}$ vs Bond number for (ac) $\theta =7.5^\circ$ and (df) $\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(ac).

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 11ef); the film is still thicker for a homoviscous drop than for a drop with $\lambda \geq 60$. Comparison between (ac) and (df) 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. (Reference Hodges, Jensen and Rallison2004) 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. (Reference Hodges, Jensen and Rallison2004) 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. (Reference Hodges, Jensen and Rallison2004) 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(al) 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. Reference Hodges, Jensen and Rallison2004).

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 (ef), 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 (ef) 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$ (gl), 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. (Reference Hodges, Jensen and Rallison2004): 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 Reference Cox1974; Claeys & Brady Reference Claeys and Brady1989; Staben, Zinchenko & Davis Reference Staben, Zinchenko and Davis2006) 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. Reference Griggs, Zinchenko and Davis2008, Reference Griggs, Zinchenko and Davis2009) 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 (Reference Rahman and Waghmare2018) 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. (Reference Hodges, Jensen and Rallison2004) 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. (Reference Hodges, Jensen and Rallison2004) (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. (Reference Hodges, Jensen and Rallison2004) 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).

Supplementary material

Supplementary material is available at https://doi.org/10.1017/jfm.2024.998.

Acknowledgements

This work utilized, in part, the RMACC Alpine supercomputer. Alpine is jointly funded by the University of Colorado Boulder, the University of Colorado Anschutz, Colorado State University and the National Science Foundation (award 2201538).

Declaration of interests

The author reports no conflict of interest.

Appendix A. Details of full double-layer desingularization

Let $g(\boldsymbol {x})$ be a smooth scalar field on a smooth closed surface $S$ with the outward unit normal $\boldsymbol {n}(\boldsymbol {x})$, and $\varphi (\boldsymbol {x})$ be a harmonic function regular inside and around $S$. Then, for any $\boldsymbol {x}^\ast \in S$

(A1)\begin{equation} \int_S g(\boldsymbol{x})\varphi(\boldsymbol{x})\,\mbox{d}S= \int_S[g(\boldsymbol{x})-g^\ast(\boldsymbol{n}^\ast\boldsymbol{\cdot}\boldsymbol{n})]\varphi(\boldsymbol{x})\,\mbox{d}S + \int_S g^\ast[(\boldsymbol{x}-\boldsymbol{x}^\ast)\boldsymbol{\cdot}\boldsymbol{n}^\ast] \frac{\partial\varphi}{\partial n}\,\mbox{d}S. \end{equation}

Here and below, the asterisk denotes the values at $\boldsymbol {x}=\boldsymbol {x}^\ast$; where the argument of $\boldsymbol {n}$, $g$ (or other quantities) is omitted for brevity, it is assumed to be the integration point $\boldsymbol {x}$. To derive (A1), Green's formula is applied to the harmonic functions $w(\boldsymbol {x})=(\boldsymbol {x}-\boldsymbol {x}^\ast )\boldsymbol {\cdot }\boldsymbol {n}^\ast$ and $\varphi (\boldsymbol {x})$

(A2)\begin{equation} \int_S\varphi\frac{\partial w}{\partial n}\mbox{d}S= \int_S\varphi(\boldsymbol{n}^\ast\boldsymbol{\cdot}\boldsymbol{n})\,\mbox{d}S= \int_S w\frac{\partial\varphi}{\partial n}\mbox{d}S= \int_S\left[(\boldsymbol{x}-\boldsymbol{x}^\ast)\boldsymbol{\cdot}\boldsymbol{n}^\ast\right] \frac{\partial\varphi}{\partial n}\mbox{d}S, \end{equation}

which immediately gives (A1). Note that, if $\varphi (\boldsymbol {x})$ is near singular at $\boldsymbol {x}\approx \boldsymbol {x}^\ast$, this near singularity is suppressed in the right-hand side integrals (A1), since $g-g^\ast (\boldsymbol {n}^\ast \boldsymbol {\cdot }\boldsymbol {n})= O(\|\boldsymbol {x}-\boldsymbol {x}^\ast \|)$ and $(\boldsymbol {x}-\boldsymbol {x}^\ast )\boldsymbol {\cdot }\boldsymbol {n}^\ast = O(\|\boldsymbol {x}-\boldsymbol {x}^\ast \|^2)$ for $\boldsymbol {x}\to \boldsymbol {x}^\ast$.

Another identity, necessary to derive (3.11), holds for an arbitrary $\boldsymbol {x}^\ast \in S$, arbitrary vector field $\boldsymbol {f}(\boldsymbol {x})= \boldsymbol {f}_\parallel (\boldsymbol {x}) +f_n(\boldsymbol {x})\boldsymbol {n}(\boldsymbol {x})$ on $S$, decomposed into the tangential ($\boldsymbol {f}_\parallel$) and normal ($\,f_n= \boldsymbol {f}\boldsymbol {\cdot }\boldsymbol {n}$) components, and arbitrary Stokes velocity field $\boldsymbol {G}(\boldsymbol {x})$ regular inside and around $S$

(A3)\begin{align} \int_S\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{G}(\boldsymbol{x})\,\mbox{d}S &= \int_S[\,\boldsymbol{f}-(\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{n}^\ast)\boldsymbol{f}_\parallel^\ast{-} (\boldsymbol{f}_\parallel^\ast\boldsymbol{\cdot}\boldsymbol{n})\boldsymbol{n}^\ast{-}f_n^\ast\boldsymbol{n}]\boldsymbol{\cdot} \boldsymbol{G}(\boldsymbol{x})\,\mbox{d}S \nonumber\\ &\quad +\int_S[(\boldsymbol{x}-\boldsymbol{x}^\ast)\boldsymbol{\cdot}\boldsymbol{n}^\ast] \,\boldsymbol{f}_\parallel^\ast\boldsymbol{\cdot}\boldsymbol{\tau}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{n}\,\mbox{d}S. \end{align}

Here, $\boldsymbol {\tau }$ is the stress tensor for flow $\boldsymbol {G}(\boldsymbol {x})$ (assuming unit viscosity)

(A4)\begin{equation} \tau_{ij}(\boldsymbol{x})={-}p(\boldsymbol{x})\delta_{ij} + \boldsymbol{\nabla}_i G_j(\boldsymbol{x})+\boldsymbol{\nabla}_j G_i(\boldsymbol{x}), \end{equation}

with the corresponding pressure field $p(\boldsymbol {x})$ and Cartesian derivatives $\boldsymbol {\nabla }_i=\partial /\partial x_i$; relation (A3) was derived in Zinchenko & Davis (Reference Zinchenko and Davis2017) irrespective of the form of $\boldsymbol {G}(\boldsymbol {x})$.

Using a linear vector field ${\boldsymbol{\mathcal L}}(\boldsymbol {x}-\boldsymbol {x}^\ast )$ (3.12), Gauss’ theorem and zero divergence of $\boldsymbol {\tau }$, we have

(A5)\begin{equation} \int_S{\mathcal{L}}_i\tau_{ij}n_j\mbox{d}S= {\mathcal{L}}_i^j\int_V ({-}p\delta_{ij}+\boldsymbol{\nabla}_i G_j+ \boldsymbol{\nabla}_j G_i)\,\mbox{d}V. \end{equation}

A part of the integral (A5) over the drop volume $V$ can be transformed back to a surface integral

(A6)\begin{equation} {\mathcal{L}}_i^j\int_V (\boldsymbol{\nabla}_i G_j+ \boldsymbol{\nabla}_j G_i)\,\mbox{d}V= ({\mathcal{L}}_i^j+ {\mathcal{L}}_j^i) \int_S n_i G_j\,\mbox{d}S, \end{equation}

which has the left-hand side form (A3) and can be replaced by the right-hand side, if the vector field $\boldsymbol {f}(\boldsymbol {x})$ is defined by (3.13). Note that $f_n^\ast =0$ is dictated by our construction of ${\boldsymbol{\mathcal L}}(\boldsymbol {x}-\boldsymbol {x}^\ast )$.

To apply all these prerequisites to the derivation of (3.11), $\boldsymbol {G}(\boldsymbol {x})$ and $\boldsymbol {\tau }(\boldsymbol {x})$ in the above relations must be replaced with the wall-correction parts $(\boldsymbol {G}^k)^{C}(\boldsymbol {x};\boldsymbol {y})$ and $(\boldsymbol {\tau }^k)^{C}(\boldsymbol {x};\boldsymbol {y})$ of the Green vector $\boldsymbol {G}^k$ and related fundamental stresslet $\tau ^k_{ij}$, respectively, for a fixed $k=1,2,3$. The corresponding wall-correction pressure part $(p^k)^{C}(\boldsymbol {x};\boldsymbol {y})$, to replace $p(\boldsymbol {x})$ in (A5), can be written as

(A7)\begin{equation} (p^k)^{C}(\boldsymbol{x};\boldsymbol{y})= \boldsymbol{\nabla}_k\varphi^k(\boldsymbol{x};\boldsymbol{y}), \end{equation}

with the auxiliary function (3.14), using the explicit form Pozrikidis (Reference Pozrikidis1992) and Blake (Reference Blake1971) for $(p^k)^{C}$. Note there is no summation here over $k$. Accordingly, the remaining part of the volume integral (A5) can be also transformed into a surface integral and handled by (A1) (with $g(\boldsymbol {x})$ and $\varphi (\boldsymbol {x})$ replaced by $n_k(\boldsymbol {x})$ and $\varphi ^k(\boldsymbol {x};\boldsymbol {y})$, respectively), since $\varphi ^k$ is a harmonic function of $\boldsymbol {x}$. As a result, all the ${\mathcal {L}}$-related contributions to the right-hand side of (3.11) total to zero, which completes the proof of this identity.

It is useful to note that the above novel approach also resolves the long-standing issue of full desingularization of double-layer BIs over an arbitrary deformable surface with the free-space fundamental stresslet, when the observation point $\boldsymbol {y}$ is slightly outside this surface. Such issues arise, e.g. in challenging BI problems for near-contact drop–drop or drop–particle interactions in free space or a periodic box. The only changes, compared with (3.11), are in the choice of $\boldsymbol {x}^\ast$ to be the nearest mesh point to $\boldsymbol {y}$, and in a simpler form for $\varphi ^k(\boldsymbol {x};\boldsymbol {y})$, namely, $\varphi ^k=1/(4{\rm \pi} r)$ for all $k=1,2,3$ (with our normalization of the free-space Green function).

Appendix B. Details of meshing algorithms and discretization

B.1. Projective meshing

At each time moment $t$, the projection centre $\boldsymbol {O}(t)=(x_1^c,x_2^c,O_3)$ is placed between the drop surface centroid $\boldsymbol {x}^c$ and the wall by the rule $O_3=x_3^c -c/x_3^c$. The empirical factor $c$ was chosen to be 0.18, 0.28, 0.40 and 0.5 for $B=0.5$, $0.25$, 0$.125$ and $0.0625$, respectively, regardless of $\lambda$, tilt angle $\theta$ and the resolution. Of the two schemes detailed in Zinchenko & Davis (Reference Zinchenko and Davis2005) for non-iterative calculation of the normal vector $\boldsymbol {n}$ and local curvature $k$ in the mesh nodes on $S$, the fourth-order scheme for $\boldsymbol {n}(\boldsymbol {x}^j)$ and $k(\boldsymbol {x}^j)$ is chosen herein, requiring two layers of mesh nodes around $\boldsymbol {x}^j$. (Note the typos in their (3.17): there must be $(\partial R/\partial x)^2$ and $(\partial R/\partial y)^2$ in lines 1 and 3, respectively, instead of $\partial R/\partial x$ and $\partial R/\partial y$.) Computation of regularized BIs on $S(t)$ is reduced to integration over $\varOmega (t)$, using $\mbox {d}S= \rho ^3\mbox {d}\varOmega /(\boldsymbol {\rho }\boldsymbol {\cdot }\boldsymbol {n})$ with $\boldsymbol {\rho }= \boldsymbol {x}-\boldsymbol {O}$. Contribution of each geodesic mesh triangle $\triangle _{sph}\in \varOmega (t)$ is calculated as outlined in Appendix A of Zinchenko & Davis (Reference Zinchenko and Davis2013) (high-order version) based on the integrand values in the vertices of the corresponding triangle $\triangle \in S(t)$ and three more mesh nodes $\boldsymbol {x}^j$ around $\triangle$. The node velocities $\boldsymbol {V}_j$ to update the mesh nodes $\boldsymbol {x}^j$ are calculated as (Zinchenko & Davis Reference Zinchenko and Davis2005)

(B1)\begin{equation} \boldsymbol{V}_j= \boldsymbol{\dot O}+\frac{\boldsymbol{u}(\boldsymbol{x^j})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x}^j) - \boldsymbol{\dot O}\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x}^j)} {(\boldsymbol{x}^j -\boldsymbol{O})\boldsymbol{\cdot} \boldsymbol{n}(\boldsymbol{x}_j)} (\boldsymbol{x}_j -\boldsymbol{O}), \end{equation}

to satisfy the constraint imposed by the BI solution and also keep all the vectors $\boldsymbol {\rho }_j/\|\rho _j\|$ of the parametric mesh stationary, i.e. $\mbox {d}\boldsymbol {\rho }_j/\mbox {d}t\parallel \boldsymbol {\rho }_j$.

The projection centre velocity $\boldsymbol {\dot O}$ can be expressed in terms of the drop centroid velocity $\boldsymbol {\dot x}^c$. For the purposes of mesh control (B1), it was sufficient to regard the mesh triangles on $S$ as flat, which gives $\boldsymbol {\dot x}^c$ as a relatively simple function of all nodes $\boldsymbol {x}^j$ and their time derivatives $\boldsymbol {V}_j$, and so (B1) can be satisfied precisely by a few iterations. However, as in Zinchenko & Davis (Reference Zinchenko and Davis2005), $\boldsymbol {\dot x}^c$ is almost insensitive to the tangential components of $\boldsymbol {V}_j$, and it was sufficient to take these components from the preceding time step, in addition to the normal components provided at the current step by the BI solution, for very accurate calculation of $\boldsymbol {\dot x}^c$ and $\boldsymbol {\dot O}$ in (B1) without iterations. That way, in the reference frame moving with the projection centre $\boldsymbol {O}$, all the nodes $\boldsymbol {\rho }_j/\|\rho _j\|$ of the parametric mesh remained stationary to about five digits for the entire simulation, from $t=0$ to steady state. Note how projective meshing does not require surface interpolations whatsoever.

The above values of $c$ were found empirically to provide moderate mesh adaptivity to the near-contact zone, with maximum-to-minimum mesh edge ratio (over the entire drop surface) $r_{max}=2.9$, $3.3$, $4.0$ and $4.7$ for $B=0.5$, $0.25$, $0.125$ and $0.0625$, respectively, at the steady state. Although optimal adaptivity is unknown and would require more experimenting, overly aggressive mesh node concentration in the film region (e.g. with $r_{max}=12$) was clearly disadvantageous, making it more difficult to reach convergence for the drop steady-state velocity with respect to $N_\triangle$; much larger $r_{max}$ could even lead to a crash. Also, it was not a robust approach in the present work to base mesh adaptivity on the local value of the drop–wall surface clearance.

B.2. The $n_3$-based passive mesh stabilization

The ‘kinetic mesh energy function’ is defined as

(B2)\begin{equation} F=\sum_{\boldsymbol{x}_{ij}}\left[\frac{\mbox{d}}{\mbox{d}t} \left (\frac{\boldsymbol{x}_{ij}^2}{h_{ij}^2} + \frac{h_{ij}^2}{\boldsymbol{x}_{ij}^2}\right)\right]^2 + 0.4\sum_\triangle\frac{1}{C_\triangle^2} \left(\frac{\mbox{d}C_\triangle}{\mbox{d}t}\right)^2. \end{equation}

In the second term of (B2) (which resists mesh triangle quality deterioration), summation is over all mesh triangles $\triangle$, where $C_\triangle =S_\triangle /(a^2 + b^2 + c^2)$ is the compactness of triangle $\triangle$ with flat area $S_\triangle$ and sides $a,b,c$. The precise value of the coefficient before the second sum (B2) is unimportant, it should be just $O(1)$ for intermediate $B$.

In the first term of (B2), summation is over all mesh edges $\boldsymbol {x}_{ij}= \boldsymbol {x}^j-\boldsymbol {x}^i$ between directly connected nodes, and this term works to prevent significant deviation of $\|\boldsymbol {x}_{ij}\|$ from local target values $h_{ij}$. With an empirical coefficient $\alpha$ close to 1

(B3)\begin{equation} h_{ij}^2= \tfrac{1}{2}(h_i^2 + h_j^2),\quad h_i^2= K[1+\alpha n_3(\boldsymbol{x}^i)]. \end{equation}

The factor $K$ is calculated at each time step as

(B4)\begin{equation} K= \frac{4}{\sqrt{3}N_\triangle}\sum_{\boldsymbol{x}^j\in S} \frac{\Delta S_j}{1+\alpha n_3(\boldsymbol{x}^j)}, \end{equation}

as if all mesh triangles were ideally equilateral with side $h$ slowly varying along the surface. Now, the area element $\Delta S_j$ associated with node $\boldsymbol {x}^j$, is simply 1/3 of the sum of flat mesh triangle areas sharing the node $\boldsymbol {x}^j$, corresponding to the standard trapezoidal integration rule. These area elements are also used in numerical evaluation (3.15) of the desingularized BIs to solve the BI problem; the local curvature and normal vector are calculated in the mesh nodes by the best paraboloid-spline method of Zinchenko & Davis (Reference Zinchenko and Davis2000).

Regardless of the form (B3) for the target values $h_{ij}$, a great simplification (Zinchenko & Davis Reference Zinchenko and Davis2013) can be made to the kinetic energy (B2), still allowing for sufficient, adaptive mesh control. Namely, $h_{ij}$ is regarded as time independent in the differentiation (B2), since spatial variation of the internode distances $\|\boldsymbol {x}_{ij}\|$ is far more important. With this simplification, (B2) can be handily written as a quadratic function of all the velocities $\boldsymbol {V}_i$ and conditionally minimized by conjugate-gradient iterations (as in the simplest, original passive mesh stabilization version of Zinchenko et al. (Reference Zinchenko, Rother and Davis1997)).

Assuming that the drop starts from spherical, the initial distribution of mesh nodes is made compatible with the form (B3), although not exactly. This form implies that the mesh size on the top of the spherical drop must be $p$ times larger than on the bottom, with $p=[(1+\alpha )/(1-\alpha )]^{1/2}$. This minimal requirement is satisfied if the initial mesh on the drop is obtained by projection from a uniform mesh of triangles on a unit sphere (§ 3.3.1) centred at $\boldsymbol {O}$, provided that the projection centre $\boldsymbol {O}$ is placed closer to the wall than the drop centre by $(p-1)/(p+1)$.

Again, to avoid overadaptive meshing, the default $\alpha$ coefficient is set to 0.7. With this choice, the steady-state mesh edge ratio $r_{max}$ was about 3.0 and 3.5 in the $B=1$ and $B=2$ simulations, respectively.

B.3. Non-adaptive meshing

The form of Zinchenko & Davis (Reference Zinchenko and Davis2002) (Appendix A therein) is used for the kinetic mesh energy function

(B5)\begin{equation} F=\sum_{\boldsymbol{x}_{ij}}\frac{1}{\|\boldsymbol{x}_{ij}\|^4} \left[\frac{\mbox{d}}{\mbox{d}t}\|\boldsymbol{x}_{ij}\|^2\right]^2 + 8\sum_\triangle\frac{1}{C_\triangle^2} \left(\frac{\mbox{d}C_\triangle}{\mbox{d}t}\right)^2. \end{equation}

To within a factor, the second term in (B5) is the same as in (B2) and helps to maintain the quality of mesh triangles. The curvature, normal vectors and area elements $\triangle S_j$ are calculated in the mesh nodes the same way as for $n_3$-based meshing (§ B2).

Appendix C. Lamb's singular series for patch contributions

Unlike in the rest of the paper, this section assumes, for more compact algebra, Pozrikidis’ (Reference Pozrikidis1992) normalization of the Green function and related fundamental pressure. Accordingly, all Lamb's singular series generated as described in this section, will need an additional multiplication by $-1/(8{\rm \pi} )$ to be used in the algorithm of § 3.

C.1. Double-layer wall-correction contributions

Technically, the most complicated element is generation of Lamb's series (3.21) for a double-layer patch contribution coming from the wall-correction part $\boldsymbol {\tau }^{C}$ of the fundamental stresslet. In the systematic approach below, allowing one to generate, without approximations, an arbitrary number of terms of the expansion (3.21), the first few, general steps almost parallel the scheme developed by Zinchenko & Davis (Reference Zinchenko and Davis2008) for free-space stresslet contributions, but then the algebra becomes far more involved in the present case, due to the cumbersome form of $\boldsymbol {\tau }^C(\boldsymbol {x};\boldsymbol {y})$ (see (3.5)). Additional care must be taken from the beginning, however, because $\boldsymbol {\tau }^C$ depends on both $\boldsymbol {x}$ and $\boldsymbol {y}$ (not only on $\boldsymbol {x}-\boldsymbol {y}$) and is lacking symmetry properties inherent in the free-space case. To simplify the notations for the following derivations, let

(C1)\begin{equation} \boldsymbol{\varPsi}(\kern0.09em \boldsymbol{y})=\sum_{\boldsymbol{x}\in\mathcal{B}} Q_k(\boldsymbol{x})W_s(\boldsymbol{x}) \boldsymbol{\tau}_{ks}(\boldsymbol{x};\boldsymbol{y}), \end{equation}

where the summation is over all mesh nodes $\boldsymbol {x}$ within a patch $\mathcal {B}$ (lying in the upper half-space $x_3>0$), and only the wall-correction part is retained in $\boldsymbol {\tau }_{ks}= (\tau _{ks}^1,\tau _{ks}^2,\tau _{ks}^3)$. Accordingly,

(C2)\begin{equation} \boldsymbol{\tau}_{ks}(\boldsymbol{x};\boldsymbol{y}) ={-} {\boldsymbol{\mathcal P}}(\boldsymbol{x};\boldsymbol{y})\delta_{ks} + \boldsymbol{\nabla}_s\boldsymbol{G}_k(\boldsymbol{x};\boldsymbol{y}) + \boldsymbol{\nabla}_k\boldsymbol{G}_s(\boldsymbol{x};\boldsymbol{y}), \end{equation}

where the Green vectors $\boldsymbol {G}_k= (G_k^1, G_k^2,G_k^3)$ and the corresponding vector of pressures ${\boldsymbol{\mathcal P}}$ only include parts coming from the wall correction. Unless otherwise stated, the Cartesian derivatives $\boldsymbol {\nabla }_k= \partial /\partial x_k$ are with respect to $\boldsymbol {x}$.

The Stokes equations $\nabla ^2 \boldsymbol {G}_k=\boldsymbol {\nabla }_k {\boldsymbol{\mathcal P}}$ and harmonicity of ${\boldsymbol{\mathcal P}}(\boldsymbol {x};\boldsymbol {y})$ (with respect to $\boldsymbol {x}$) imply $\nabla ^2 \boldsymbol {\tau }_{ks}(\boldsymbol {x};\boldsymbol {y})= 2\boldsymbol {\nabla }_k\boldsymbol {\nabla }_s{\boldsymbol{\mathcal P}}(\boldsymbol {x};\boldsymbol {y})$ and harmonicity of auxiliary vector fields

(C3)\begin{equation} \boldsymbol{t}_{ks}(\boldsymbol{x};\boldsymbol{y})=\boldsymbol{\tau}_{ks}(\boldsymbol{x};\boldsymbol{y}) - \tfrac{1}{4}[(x-x^o)_k \tilde{{\boldsymbol{\mathcal P}}}_s(\boldsymbol{x};\boldsymbol{y})+ (x-x^o)_s \tilde{{\boldsymbol{\mathcal P}}}_k(\boldsymbol{x};\boldsymbol{y})], \end{equation}

where the vectors $\tilde {{\boldsymbol{\mathcal P}}}_k(\boldsymbol {x};\boldsymbol {y})$ are defined as $2\boldsymbol {\nabla }_k{\boldsymbol{\mathcal P}}(\boldsymbol {x};\boldsymbol {y})$, and $\boldsymbol {x}^o$ is chosen as the centre of the minimal spherical shell around the block $\mathcal {B}$.

Next, a special form (Zinchenko & Davis Reference Zinchenko and Davis2000) of the Taylor series for harmonic functions is used

(C4)\begin{equation} f(\boldsymbol{x})= \sum_{\nu=0}^\infty \sum_{\mu={-}\nu}^\nu \partial_{\nu,\mu}f(\boldsymbol{x})|_{\boldsymbol{x}=\boldsymbol{x}^o} Z_{\nu,\mu}(\boldsymbol{x}-\boldsymbol{x}^o). \end{equation}

Here, the differential operator $\partial _{\nu,\mu }$ is defined for $0\leq |\mu |\leq \nu$ as

(C5)\begin{equation} \partial_{\nu,\mu}= (\boldsymbol{\nabla}_1 -\mbox{i}\boldsymbol{\nabla}_2)^\mu \boldsymbol{\nabla}_3^{\nu-|\mu|}, \end{equation}

(with $(\boldsymbol {\nabla }_1 -\mbox {i}\boldsymbol {\nabla }_2)^\mu = (-1)^{\mu } (\boldsymbol {\nabla }_1 +\mbox {i}\boldsymbol {\nabla }_2)^{-\mu }$ when $\mu <0$) $\mbox {i}=\sqrt {-1}$, and positive-order solid harmonics $Z_{\nu,\mu }(\boldsymbol {r})$ are related to standard normalized spherical harmonics ${\mathcal {Y}}_{\nu,\mu }(\boldsymbol {r})$ as

(C6)\begin{equation} Z_{\nu,\mu}(\boldsymbol{r})=\frac{2{\rm \pi}^{1/2}r^\nu{\mathcal{Y}}_{\nu,\mu}(\boldsymbol{r})} {[(2\nu+1)(\nu-\mu)!(\nu+\mu)!]^{1/2}}. \end{equation}

Relation (C4) can be shown to be equivalent to (27) of Sangani & Mo (Reference Sangani and Mo1996).

Substituting $\boldsymbol {\tau }_{ks}(\boldsymbol {x};\boldsymbol {y})$ from (C3) into (C1) and applying (C4) to the harmonic functions $\boldsymbol {t}_{ks}(\boldsymbol {x};\boldsymbol {y})$, $\tilde {{\boldsymbol{\mathcal P}}}_k(\boldsymbol {x};\boldsymbol {y})$ and $\tilde {{\boldsymbol{\mathcal P}}}_s(\boldsymbol {x};\boldsymbol {y})$ yields

(C7)\begin{equation} \boldsymbol{\varPsi}(\kern0.09em \boldsymbol{y})=\sum_{\nu=0}^\infty \sum_{\mu={-}\nu}^\nu [\tilde{E}_{\nu,\mu,k,s}\partial_{\nu,\mu}\boldsymbol{t}_{ks}(\boldsymbol{x};\boldsymbol{y}) +\tilde{D}_{\nu,\mu,k}\partial_{\nu,\mu}\tilde{{\boldsymbol{\mathcal P}}}_k (\boldsymbol{x};\boldsymbol{y})]_{\boldsymbol{x}=\boldsymbol{x}^0}. \end{equation}

Here, the patch double-layer moments are

(C8)\begin{equation} \tilde{E}_{\nu,\mu,k,s}=\sum_{\boldsymbol{x}\in\mathcal{B}}W_{(s}Q_{k)}(\boldsymbol{x}) Z_{\nu,\mu}(\boldsymbol{x}-\boldsymbol{x}^o), \end{equation}

and

(C9)\begin{equation} \tilde{D}_{\nu,\mu,k}=\frac{1}{4}\sum_{\boldsymbol{x}\in\mathcal{B}} (\boldsymbol{x}-\boldsymbol{x}^o)\boldsymbol{\cdot} [Q_k(\boldsymbol{x})\boldsymbol{W}(\boldsymbol{x})+W_k(\boldsymbol{x})\boldsymbol{Q}(\boldsymbol{x})] Z_{\nu,\mu}(\boldsymbol{x}-\boldsymbol{x}^o), \end{equation}

with symmetrization $W_{(s}Q_{k)}= (W_s Q_k+W_k Q_s)/2$ in (C8).

The next step is to introduce the auxiliary vectors

(C10)\begin{equation} \boldsymbol{g}_k(\boldsymbol{x};\boldsymbol{y})=\boldsymbol{G}_k(\boldsymbol{x};\boldsymbol{y}) - \tfrac{1}{2}(x-x^o)_k{\boldsymbol{\mathcal P}}(\boldsymbol{x};\boldsymbol{y})\end{equation}

(still only keeping the parts of $\boldsymbol {G}$ and ${\boldsymbol{\mathcal P}}$ coming from the wall correction). These vectors, as well as ${\boldsymbol{\mathcal P}}(\boldsymbol {x};\boldsymbol {y})$, are harmonic functions of $\boldsymbol {x}$. Also, as follows from (C2) and (C10), $\boldsymbol {t}_{ks}= \boldsymbol {\nabla }_k \boldsymbol {g}_s + \boldsymbol {\nabla }_s \boldsymbol {g}_k$. For these reasons, the high-order derivatives in (C7) can be expressed in terms of $\partial _{\nu +1,\mu '}\boldsymbol {g}_k$, $\partial _{\nu +1,\mu '}\boldsymbol {g}_s$ and $\partial _{\nu +1,\mu '}{\boldsymbol{\mathcal P}}_k$ (with $\mu '=\mu -1$, $\mu$ or $\mu +1$) by the same algebra as in Zinchenko & Davis (Reference Zinchenko and Davis2008, § 4.2 therein) to arrive at a more tractable form

(C11)\begin{equation} \boldsymbol{\varPsi}(\kern0.09em \boldsymbol{y})=\sum_{\nu=1}^\infty \sum_{\mu={-}\nu}^\nu [E_{\nu,\mu,k}\partial_{\nu,\mu}\boldsymbol{g}_{k}(\boldsymbol{x};\boldsymbol{y})+ D_{\nu,\mu}\partial_{\nu,\mu}{\boldsymbol{\mathcal P}}(\boldsymbol{x};\boldsymbol{y})]_{\boldsymbol{x}=\boldsymbol{x}^o}. \end{equation}

The new $E$- and $D$-coefficients are calculated as

(C12)\begin{align} E_{\nu,\mu,k}=\tilde{E}_{\nu-1,\mu-1,k,1} + \mbox{i} \tilde{E}_{\nu-1,\mu-1,k,2} - \tilde{E}_{\nu-1,\mu+1,k,1}+ \mbox{i} \tilde{E}_{\nu-1,\mu+1,k,2}+ 2\tilde{E}_{\nu-1,\mu,k,3}, \end{align}

and

(C13)\begin{equation} D_{\nu,\mu}=\tilde{D}_{\nu-1,\mu-1,1} + \mbox{i} \tilde{D}_{\nu-1,\mu-1,2} - \tilde{D}_{\nu-1,\mu+1,1}+ \mbox{i} \tilde{D}_{\nu-1,\mu+1,2}+ 2\tilde{D}_{\nu-1,\mu,3}\end{equation}

(assuming that any coefficients $\tilde {E}_{\nu ',\mu ',k,s}$ or $\tilde {D}_{\nu ',\mu ',k}$ with $|\mu '|>\nu '$ appearing on the right-hand side of (C12) or (C13) are set to zero).

The explicit form of $\boldsymbol {g}_{k}$ (not as complicated as for $\boldsymbol {\tau }_{ks}$) can be derived from our definition (C10) and the relations from Blake (Reference Blake1971) or Pozrikidis (Reference Pozrikidis1992) for $\boldsymbol {G}$ and ${\boldsymbol{\mathcal P}}$

(C14)\begin{equation} \boldsymbol{g}_k(\boldsymbol{x};\boldsymbol{y})=2y_3\left[\boldsymbol{\nabla}_3\frac{1}{R}\boldsymbol{e}_k - y_3\boldsymbol{\nabla}\boldsymbol{\nabla}_k\frac{1}{R} -\delta_{3k}\boldsymbol{\nabla}\frac{1}{R}+ R^o_k\boldsymbol{\nabla}\boldsymbol{\nabla}_3\frac{1}{R}\right]^{IM} -\frac{1}{R}\boldsymbol{e}_k + R^o_k\boldsymbol{\nabla}\frac{1}{R}. \end{equation}

Here, $R=\|\boldsymbol {x}- \boldsymbol {y}^{IM}\|$, $\boldsymbol {R}^o=\boldsymbol {x}^o-\boldsymbol {y}^{IM}$ and the superscript ${IM}$ stands for the mirror image of a vector or a point with respect to the wall. Applying the operator $\partial _{\nu,\mu }$ to (C14), one can obtain

(C15)\begin{align} & E_{\nu,\mu,k}\partial_{\nu,\mu}\boldsymbol{g}_{k}(\boldsymbol{x};\boldsymbol{y})= 2y_3\left[\left(\boldsymbol{\nabla}_3\partial_{\nu,\mu}\frac{1}{R}\right) \boldsymbol{E}_{\nu,\mu} - y_3\boldsymbol{\nabla}\left(\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{\nabla} \partial_{\nu,\mu}\frac{1}{R}\right)-E_{\nu,\mu,3}\boldsymbol{\nabla}\partial_{\nu,\mu}\frac{1}{R} \right.\nonumber\\ &\quad +\left.(\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R}^o)\boldsymbol{\nabla}\left(\boldsymbol{\nabla}_3 \partial_{\nu,\mu}\frac{1}{R}\right)\right]^{IM} - \left(\partial_{\nu,\mu}\frac{1}{R}\right)\boldsymbol{E}_{\nu,\mu} + (\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R}^o)\boldsymbol{\nabla}\partial_{\nu,\mu}\frac{1}{R}, \end{align}

where, for brevity, $\boldsymbol {E}_{\nu,\mu }= (E_{\nu,\mu,1},E_{\nu,\mu,2}, E_{\nu,\mu,3})$; the expression (C15) must be used at $\boldsymbol {x}=\boldsymbol {x}^o$ to give the first contribution $\boldsymbol {\varPsi }_{\nu,\mu }^{g}(\kern0.09em \boldsymbol {y})$ in the brackets of (C11).

For the wall-correction part ${\boldsymbol{\mathcal P}}= 4y_3[\boldsymbol {\nabla }\boldsymbol {\nabla }_3({1}/{R})]^{IM}+ 2\boldsymbol {\nabla }({1}/{R})$ of the vector of pressures, we have a simpler expression

(C16)\begin{equation} D_{\nu,\mu}\partial_{\nu,\mu}{\boldsymbol{\mathcal P}}(\boldsymbol{x};\boldsymbol{y})= \boldsymbol{\varPsi}_{\nu,\mu}^{\mathcal{P}}(\kern0.09em \boldsymbol{y})= D_{\nu,\mu}\left[4y_3\left(\boldsymbol{\nabla}\boldsymbol{\nabla}_3\partial_{\nu,\mu} \frac{1}{R}\right)^{IM} + 2\boldsymbol{\nabla}\partial_{\nu,\mu}\frac{1}{R}\right], \end{equation}

again to be used at $\boldsymbol {x}=\boldsymbol {x}^o$.

To convert (C11) into Lamb's singular series (3.21) for the Stokes flow $\boldsymbol {\varPsi }(\kern0.09em \boldsymbol {y})$ and to generate the harmonics $p_{-(n+1)}(\boldsymbol {Y})$, $\varPhi _{-(n+1)}(\boldsymbol {Y})$ and $\chi _{-(n+1)}(\boldsymbol {Y})$, we use general relations from Happel & Brenner (Reference Happel and Brenner1986), in particular

(C17)\begin{equation} \boldsymbol{\varPsi}(\kern0.09em \boldsymbol{y})\boldsymbol{\cdot}\boldsymbol{Y}= \sum_{n=0}^\infty \left[\frac{n+1}{2(2n-1)} Y^2 p_{-(n+1)} - (n+1)\varPhi_{-(n+1)}\right]. \end{equation}

Instead of $\boldsymbol {Y}_{\gamma }$ defined in § 3.4, a simpler notation $\boldsymbol {Y}$ is sufficient here for the expansion vector $\boldsymbol {y}-(\boldsymbol {x}^o)^{IM}$.

Using $\boldsymbol {Y}^{IM}= -\boldsymbol {R}^o$, one can obtain for the contribution from (C15) to $\boldsymbol {\varPsi }(\kern0.09em \boldsymbol {y})\boldsymbol {\cdot }\boldsymbol {Y}$

(C18)\begin{align} & E_{\nu,\mu,k}\boldsymbol{Y}\boldsymbol{\cdot}\partial_{\nu,\mu}\boldsymbol{g}_{k}(\boldsymbol{x};\boldsymbol{y})={-}2y_3\left[(\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R}) \left(1+R\frac{\partial}{\partial R}\right) \left(\boldsymbol{\nabla}_3\partial_{\nu,\mu}\frac{1}{R}\right) \right.\nonumber\\ &\quad -\left.y_3 R\frac{\partial}{\partial R}\left(\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{\nabla} \partial_{\nu,\mu}\frac{1}{R}\right)-E_{\nu,\mu,3}R\frac{\partial}{\partial R} \partial_{\nu,\mu}\frac{1}{R}\right] \nonumber\\ &\quad +(\boldsymbol{E}^{IM}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R})\partial_{\nu,\mu}\frac{1}{R} +(\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R})\left({-}R\frac{\partial}{\partial R} \partial_{\nu,\mu}\frac{1}{R} + 2R_3\boldsymbol{\nabla}_3\partial_{\nu,\mu}\frac{1}{R}\right). \end{align}

Here, $R\partial /\partial R=\boldsymbol {R}\boldsymbol {\cdot }\boldsymbol {\nabla }$ contains only differentiation along $\boldsymbol {R}= \boldsymbol {x}-\boldsymbol {y}^{IM}$, which must be set to $\boldsymbol {R}^o= \boldsymbol {x}^o-\boldsymbol {y}^{IM}$ in the final result. Since $\partial _{\nu,\mu }(1/R)$ and $\boldsymbol {\nabla }\partial _{\nu,\mu }(1/R)$ are homogeneous functions of $\boldsymbol {R}$ of degrees $-(\nu +1)$ and $-(\nu +2)$, respectively, the terms with $R\partial /\partial R$ are simplified by the Euler theorem. Also, using $y_3= R_3-x^o_3$, the right-hand side of (C18) can be transformed after some algebra into

(C19)\begin{align} & (\nu+2) \left\{{-}2y_3[E_{\nu,\mu,1}(R_3\boldsymbol{\nabla}_1-R_1\boldsymbol{\nabla}_3) + E_{\nu,\mu,2}(R_3\boldsymbol{\nabla}_2-R_2\boldsymbol{\nabla}_3)-x_3^o\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{\nabla}] \partial_{\nu,\mu}\frac{1}{R} \right.\nonumber\\ &\quad + \left.(\boldsymbol{E}^{IM}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R}) \partial_{\nu,\mu}\frac{1}{R}\right\} + 2x_3^o[(\nu+1)E_{\nu,\mu,3}+ (\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}\boldsymbol{R})\boldsymbol{\nabla}_3] \partial_{\nu,\mu}\frac{1}{R}. \end{align}

In a similar, but simpler manner, one can account for the pressure contribution in (C11) to $\varPsi (\kern0.09em \boldsymbol {y})\boldsymbol {\cdot }\boldsymbol {Y}$. Using (C16), we have at $\boldsymbol {x}=\boldsymbol {x}^o$

(C20)\begin{equation} D_{\nu,\mu}\boldsymbol{Y}\boldsymbol{\cdot}\partial_{\nu,\mu}{\boldsymbol{\mathcal P}}(\boldsymbol{x};\boldsymbol{y})= D_{\nu,\mu}[4(\nu+3)R_3\boldsymbol{\nabla}_3 -4(\nu+2)x_3^o\boldsymbol{\nabla}_3 +2(\nu+1)] \partial_{\nu,\mu}\frac{1}{R}. \end{equation}

The following steps employ the recurrent relations between the fundamental derivatives $\partial _{\nu,\mu }(1/R)$ with respect to $\boldsymbol {R}$:

(C21ac)$$\begin{gather} (\boldsymbol{\nabla}_1 +\mbox{i}\boldsymbol{\nabla}_2)\partial_{\nu,\mu}={-}\partial_{\nu+1,\mu-1},\quad (\boldsymbol{\nabla}_1 -\mbox{i}\boldsymbol{\nabla}_2)\partial_{\nu,\mu}= \partial_{\nu+1,\mu+1},\quad \boldsymbol{\nabla}_3 \partial_{\nu,\mu}=\partial_{\nu+1,\mu}, \end{gather}$$
(C22)$$\begin{gather}(R_1 +\mbox{i}R_2)\partial_{\nu,\mu}= \frac{R^2}{2\nu+1} \partial_{\nu+1,\mu-1} -\frac{(\nu+\mu-1)(\nu+\mu)}{2\nu+1}\partial_{\nu-1,\mu-1}, \end{gather}$$
(C23)$$\begin{gather}(R_1 -\mbox{i}R_2)\partial_{\nu,\mu}={-}\frac{R^2}{2\nu+1}\partial_{\nu+1,\mu+1} +\frac{(\nu-\mu-1)(\nu-\mu)}{2\nu+1}\partial_{\nu-1,\mu+1}, \end{gather}$$
(C24)$$\begin{gather}R_3\partial_{\nu,\mu}={-}\frac{(\nu-\mu)(\nu+\mu)}{2\nu+1}\partial_{\nu-1,\mu}- \frac{R^2}{2\nu+1} \partial_{\nu+1,\mu} \end{gather}$$
(C25)$$\begin{gather}(R_2\boldsymbol{\nabla}_3 - R_3\boldsymbol{\nabla}_2)\partial_{\nu,\mu}=\frac{\mbox{i}}{2} [(\nu+\mu)\partial_{\nu,\mu-1} + (\nu-\mu)\partial_{\nu,\mu+1}], \end{gather}$$
(C26)$$\begin{gather}(R_3\boldsymbol{\nabla}_1 - R_1\boldsymbol{\nabla}_3)\partial_{\nu,\mu}=\frac{1}{2} [(\nu+\mu)\partial_{\nu,\mu-1} - (\nu-\mu)\partial_{\nu,\mu+1}], \end{gather}$$
(C27)$$\begin{gather}(R_1\boldsymbol{\nabla}_2 - R_2\boldsymbol{\nabla}_1)\partial_{\nu,\mu}={-}\mbox{i}\mu\partial_{\nu,\mu} \end{gather}$$

(omitting, for brevity, the argument $1/R$ of $\partial _{\nu,\mu })$. Of course, (C21ac)–(C27) are equivalent to recurrent formulae for standard negative-order solid harmonics owing to Maxwell's relation

(C28)\begin{equation} \partial_{\nu,\mu}\frac{1}{R}= ({-}1)^{\nu-\mu} \left[\frac{4{\rm \pi}(\nu-\mu)!(\nu+\mu)!}{2\nu+1}\right]^{1/2} \frac{{\mathcal{Y}}_{\nu,-\mu}(\boldsymbol{R})}{R^{\nu+1}}, \end{equation}

but using (temporarily) $\partial _{\nu,\mu }(1/R)$ instead of such harmonics makes algebra simpler and more attractive.

Using $y_3=R_3^o-x_3^o$ and the recurrent relations (C21ac)–(C27) in (C19) and (C20) allows one to write cumulative contributions from (C19) and (C20) to $\boldsymbol {\varPsi }(\kern0.09em \boldsymbol {y})\boldsymbol {\cdot }\boldsymbol {Y}$ in the form (C17) and thereby generate $p_{-(n+1)}$ and $\varPhi _{-(n+1)}$, first as combinations of $\partial _{n,m}(1/R)$. It remains to use (C28) and ${\mathcal {Y}}_{\nu,\mu }(\boldsymbol {R})/R^{\nu +1}= (-1)^\mu {\mathcal {Y}}_{\nu,\mu }(\boldsymbol {Y})/Y^{\nu +1}$ to complete calculation of the coefficients $A_{-(\nu +1),\mu }$ and $B_{-(\nu +1),\mu }$ in Lamb's singular series (3.21) for a double-layer patch contribution coming from the wall-correction part of the fundamental stresslet.

A different technique is used to generate the harmonics $\chi _{-(n+1)}$ in Lamb's series (3.21) based on the general identity (Happel & Brenner Reference Happel and Brenner1986)

(C29)\begin{equation} [\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\varPsi}(\kern0.09em \boldsymbol{y})]\boldsymbol{\cdot} \boldsymbol{Y}= \sum_{n=1}^\infty n(n+1) \chi_{-(n+1)}. \end{equation}

First, (C15) is recast in terms of differentiations with respect to the observation point. Assuming that all further derivatives till the end of this subsection are now with respect to $\boldsymbol {y}$, we have

(C30)\begin{align} \boldsymbol{\varPsi}_{\nu,\mu}^{g}(\kern0.09em \boldsymbol{y}) &= 2y_3[(\boldsymbol{\nabla}_3 f)\boldsymbol{E}_{\nu,\mu}^{IM} - y_3\boldsymbol{\nabla}(\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{\nabla}f) +E_{\nu,\mu,3}\boldsymbol{\nabla}f+ (\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{Y})\boldsymbol{\nabla}(\boldsymbol{\nabla}_3 f)] \nonumber\\ &\quad -f\boldsymbol{E}_{\nu,\mu} + (\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{Y}) \boldsymbol{\nabla}f -2(\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{Y}) (\boldsymbol{\nabla}_3 f) \boldsymbol{e}_3, \end{align}

where, for brevity, $f=(-1)^\mu \partial _{\nu,\mu }(1/Y)$. After some algebra, we have from (C30)

(C31)\begin{align} \boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\varPsi}_{\nu,\mu}^{g}(\kern0.09em \boldsymbol{y}) &= 4\boldsymbol{e}_3\boldsymbol{\times}[(\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{Y}) \boldsymbol{\nabla}(\boldsymbol{\nabla}_3 f)+(\boldsymbol{\nabla}_3 f) \boldsymbol{E}_{\nu,\mu}^{IM} - y_3\boldsymbol{\nabla} (\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{\nabla}f)] \nonumber\\ &\quad +2\boldsymbol{E}_{\nu,\mu}\boldsymbol{\times}\boldsymbol{\nabla}f, \end{align}

so that

(C32)\begin{align} & \boldsymbol{Y}\boldsymbol{\cdot}[\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\varPsi}_{\nu,\mu}^{g}(\kern0.09em \boldsymbol{y})]= 4[(\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{Y}) (Y_2\boldsymbol{\nabla}_1-Y_1\boldsymbol{\nabla}_2)+Y_2 E_{\nu,\mu,1}-Y_1 E_{\nu,\mu,2}]\boldsymbol{\nabla}_3f \nonumber\\ &\quad +[4Y_3(Y_1\boldsymbol{\nabla}_2-Y_2\boldsymbol{\nabla}_1) +4x_3^o(Y_2\boldsymbol{\nabla}_1-Y_1\boldsymbol{\nabla}_2)] (\boldsymbol{E}_{\nu,\mu}^{IM}\boldsymbol{\cdot}\boldsymbol{\nabla}f) + 2\boldsymbol{E}_{\nu,\mu}\boldsymbol{\cdot}(\boldsymbol{\nabla}f\boldsymbol{\times}\boldsymbol{Y}), \end{align}

where the split $y_3= Y_3-x_3^o$ was used. Matching the form (C32) to the expansion (C29) is not obvious and requires further detailed derivation. Using repeatedly the recurrent formulae (C21ac)–(C27) (with $\boldsymbol {Y}$ instead of $\boldsymbol {R}$), (C32) can be transformed to

(C33)\begin{align} \boldsymbol{Y}\boldsymbol{\cdot} [\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\varPsi}_{\nu,\mu}^{g}(\kern0.09em \boldsymbol{y})] &= ({-}1)^\mu\{(\nu-\mu)(2\mu+1)(\mbox{i}E_{\nu,\mu,1}-E_{\nu,\mu,2})\partial_{\nu,\mu+1} \nonumber\\ &\quad -(\nu+\mu)(2\mu-1)(\mbox{i}E_{\nu,\mu,1}+E_{\nu,\mu,2})\partial_{\nu,\mu-1} +2\mbox{i}\mu E_{\nu,\mu,3}\partial_{\nu,\mu} \nonumber\\ &\quad + 2x_3^o[(\mu+1)(\mbox{i}E_{\nu,\mu,1}-E_{\nu,\mu,2})\partial_{\nu+1,\mu+1} \nonumber\\ &\quad -(\mu-1)(\mbox{i}E_{\nu,\mu,1}+E_{\nu,\mu,2})\partial_{\nu+1,\mu-1} -2\mbox{i}\mu E_{\nu,\mu,3}\partial_{\nu+1,\mu}]\}, \end{align}

where the argument $1/Y$ of the differential operator (C5) has been omitted for brevity. The terms with $Y^2\partial _{\nu +2, \mu \pm 1}$ and $Y^2\partial _{\nu +2, \mu }$, appearing at the intermediate stage of derivations, cancel out in (C33), making it, indeed, a harmonic function.

The contribution of $\boldsymbol {\varPsi }_{\nu,\mu }^{\mathcal {P}}(\kern0.09em \boldsymbol {y})$ to the left-hand side of (C29) is easy to handle. Recasting (C16) in terms of differentiations with respect to $\boldsymbol {y}$, then taking the curl and using the recurrent relations (C21ac) and (C27), we have

(C34)\begin{equation} \boldsymbol{Y}\boldsymbol{\cdot}[\boldsymbol{\nabla}\boldsymbol{\times}\boldsymbol{\varPsi}_{\nu,\mu}^{\mathcal{P}}(\kern0.09em \boldsymbol{y})]= 8({-}1)^{\mu+1}\mbox{i}\mu D_{\nu,\mu}\partial_{\nu+1,\mu}\frac{1}{Y}. \end{equation}

Summing up (C33) and (C34), rearranging indices and comparing the result with (C29) yields the solid harmonics $\chi _{-(n+1)}(\boldsymbol {Y})$, first as combinations of $\partial _{n,m}(1/Y)$. Finally, (C28) allows us to complete calculation of the coefficients $C_{-(\nu +1),\mu }$ in Lamb's singular series (3.21) for a wall-induced double-layer patch contribution.

Apart from the validations (§ 4 and Appendix D) of the whole algorithm, correctness of our generated coefficients $A_{-(\nu +1),\mu }$, $B_{-(\nu +1),\mu }$ and $C_{-(\nu +1),\mu }$ was confirmed for randomly selected patches $\mathcal {B}$ in the upper half-space and observation points $\boldsymbol {y}$ just slightly outside the minimal spherical shell around the mirror image of $\mathcal {B}$. For randomly set $Q_k(\boldsymbol {x})$ and $W_s(\boldsymbol {x})$, calculation of (C1) by Lamb's singular series (3.27) could be made arbitrarily close to the exact result by direct point to summations with increasing the number of terms in Lamb's series; the agreement to 9 or more digits was observed.

C.2. Double-layer free-space contributions

Efficient generation of Lamb's singular series for free-space patch contributions (3.21) based on the patch double-layer moments was developed earlier (Zinchenko & Davis Reference Zinchenko and Davis2008) and later used in many multipole-accelerated BI simulations for deformable drops (e.g. Zinchenko & Davis Reference Zinchenko and Davis2013). Still, it is helpful to note how the above methodology for wall-correction contributions is easily adapted for the free-space contributions. If $\boldsymbol {\tau }_{ks}(\boldsymbol {x};\boldsymbol {y})$ in (C1) is the free-space part $-6r_kr_s\boldsymbol {r}/r^5$ of the fundamental stresslet (with $\boldsymbol {r}=\boldsymbol {x}-\boldsymbol {y}$), then (C2) to (C13) still hold, but the expressions for the vector of pressures and auxiliary functions (C10) are now much simpler

(C35a,b)\begin{equation} \boldsymbol{g}_k (\boldsymbol{x};\boldsymbol{y})=\frac{1}{r}\boldsymbol{e}_k - r^o_k\boldsymbol{\nabla}^x\frac{1}{r},\quad {\boldsymbol{\mathcal P}}(\boldsymbol{x};\boldsymbol{y})={-}2\boldsymbol{\nabla}^x\frac{1}{r}. \end{equation}

Here, $\boldsymbol {r}^o= \boldsymbol {x}^o-\boldsymbol {y}$, and the expansion vector $\boldsymbol {Y}$ is now $-\boldsymbol {r}^o$. The same logic as for wall-correction contributions, but with much less effort, applies to convert the free-space version of (C11) into Lamb's singular series (3.21).

C.3. Single-layer patch contributions

To generate Lamb's singular series (3.21) for single-layer patch contributions related to the discretized form of the right-hand side of (3.8), with either $\boldsymbol {G}^{FS}$ or $\boldsymbol {G}^{C}$, the patch moments $E_{\nu,\mu,k}$ and $D_{\nu,\mu }$ must be redefined. Instead of the operations (C12) and (C13), these coefficients are now directly calculated as

(C36a,b)\begin{align} {E}_{\nu,\mu,k}=\sum_{\boldsymbol{x}\in\mathcal{B}}W_k(\boldsymbol{x})Z_{\nu,\mu}(\boldsymbol{x}-\boldsymbol{x}^o),\quad {D}_{\nu,\mu}=\frac{1}{2}\sum_{\boldsymbol{x}\in\mathcal{B}} (\boldsymbol{x}-\boldsymbol{x}^o)\boldsymbol{\cdot} \boldsymbol{W}(\boldsymbol{x})Z_{\nu,\mu}(\boldsymbol{x}-\boldsymbol{x}^o), \end{align}

where $\boldsymbol {W}(\boldsymbol {x})$ is either $\tilde {f}(\boldsymbol {x}^j)\boldsymbol {n}(\boldsymbol {x}^j)\Delta S_j$ or $\boldsymbol {n}(\boldsymbol {x}^j)\Delta S_j$. With these definitions, the single-layer patch contributions still have the form (C11). Accordingly, the rest of § C.1, starting from (C14), applies without changes to generate Lamb's singular series (3.21) for single-layer, wall-correction patch contributions. For free-space, single-layer patch contributions, just a different form (C35a,b) of $\boldsymbol {g}_k$ and ${\boldsymbol{\mathcal P}}$ with the expansion vector $\boldsymbol {Y}=\boldsymbol {y}-\boldsymbol {x}^o$ is used.

Appendix D. Miscellaneous details of the algorithm

D.1. Time stepping

Not only the minimum mesh size limits a stable time step (as the Courant stability rule requires), but this step is also affected by elevated local surface curvatures and (in our case) the drop proximity to the wall; there are no rational and universal strategies, though, for all BI problems. Combining the semi-empirical approaches of Zinchenko & Davis (Reference Zinchenko and Davis2005, Reference Zinchenko and Davis2008), the dimensional time step was set as

(D1)\begin{equation} \Delta t=K_{\Delta t}\frac{\mu_e}{\sigma}\min(\varDelta_1, 0.7\varDelta_2), \end{equation}

where $K_{\Delta t}$ is a numerical factor

(D2)\begin{equation} \varDelta_1=\min_i\left\{\frac{\Delta x_i}{a\max[|k_1(\boldsymbol{x}^i)|,|k_2(\boldsymbol{x}^i)|]}\right\}, \end{equation}

and

(D3)\begin{equation} \varDelta_2=\min_{i,j}\|\boldsymbol{x}^i-\boldsymbol{x}^j\|,\quad \boldsymbol{x}^i\in S,\enspace \boldsymbol{x}^j\in S^{IM}. \end{equation}

In (D2), the minimum is over all mesh nodes $\boldsymbol {x}^i\in S$ with principal curvatures $k_1$ and $k_2$; the shortest distance from $\boldsymbol {x}^i$ to the nodes directly connected to it is $\Delta x_i$. In (D3), $\varDelta _2$ is the minimum node-to-node distance between $S$ and its mirror image $S^{IM}$, but excluding pairs $(i,j)$ with $\boldsymbol {x}^j$ being the nearest node to $\boldsymbol {x}^i$ (since such influences are excluded by the leading-order near-singularity subtraction in BI); the $\varDelta _2$ term rarely came into play in the present simulations due to extreme resolutions used.

The simulations used mostly $K_{\Delta t}=5{-}10$ (with typically 4–6 BI iterations per time step) to avoid sharp temporal changes. Remarkably, for high $\lambda \geq 60$, as the steady state was roughly approached, it was possible to gradually increase $K_{\Delta t}$ to much higher values ($\sim$100–200) without losing numerical stability, and accurately finalize the steady state without excessive computational expenses. This highly beneficial feature may be due to using passive mesh stabilization with its soft stability limitations.

D.2. One-time projective remeshing

For small to moderately small $B$, when projective meshing was used to discretize and track the drop surface (§ 3.3.1), there was an easy way to continue any such simulation with a new number of mesh triangles ($N_\triangle '$) and a new number of patches $M'$ (if necessary to retain efficiency of multipole acceleration) by one-time remeshing of the parent triangulation (of $N_\triangle$ elements), while inheriting near-contact adaptivity from the parent mesh. This is, for example, how curve 3 in figure 9(a) with $N_\triangle =138$ K was continued from curve 2 (with $N_\triangle =78$ K), instead of doing the 138 K resolution run from scratch. First, a new uniform, parametric mesh is generated on the unit sphere $\varOmega (t)$ centred at $\boldsymbol {O}(t)$, with $N_\triangle '$ elements and nodes partitioned into $M'$ patches. The radial distance $\rho =\|\boldsymbol {x}-\boldsymbol {O}(t)\|$ for nodes $\boldsymbol {x}\in S$ of the parent mesh, which serves mapping between $S(t)$ and $\varOmega (t)$ (figure 2b), is interpolated from the old to new parametric mesh. To this end, a local fourth-order polynomial of the two coordinates in the tangential plane at an old node $\boldsymbol {\xi }\in \varOmega (t)$ is used, best fitted to the values of $\ln \rho$ in the two layers of nodes around $\boldsymbol {\xi }$. Once $\rho$ is known on the new parametric mesh, it gives the new nodes on $S(t)$ to continue the simulation. The validity and accuracy of this approach was verified in selected cases by repeating the new-resolution simulation from scratch, only to arrive at the same steady-state results.

D.3. Economical truncation bounds

Construction of bounds for patch contributions with $\tilde {f}(\boldsymbol {x})$ on the right-hand side of (3.8), and for double-layer expansions (3.21), closely follows that of Zinchenko & Davis (Reference Zinchenko and Davis2000) for many deformable drops composed of blocks. Namely, their blocks are replaced here by an extended system of patches (figure 6), and only their near-field bounds are relevant to our problem. The strength factors $C_\gamma$ and $\tilde {C}_\gamma$ given by (3.85) and (3.92) therein are additionally increased 4 times for $x_3<0$ patches to account for complexity of $\boldsymbol {G}^{C}$ and $\boldsymbol {\tau }^{C}$. In other respects, their (3.82) to (3.92) are applied without changes to give the bounds for singular-to-regular Lamb's series re-expansion, and for pointwise calculation of Lamb's singular series. With $e_{nf}=0.1$ and $\tilde {e}_{nf}=10$ chosen here, these truncation bounds are controlled by a single intuitive precision parameter $\varepsilon$. They are calculated before the iterations (using $\boldsymbol {Q}$ from the preceding time step to estimate $\tilde {C}_\gamma$ for double-layer expansions (3.21)).

A modification of this scheme to construct truncation bounds for patch contributions associated with $\tilde {f}(\kern0.09em \boldsymbol {y})$ and $\tilde {f}(\boldsymbol {x}^\ast )$ on the right-hand side of (3.8), and with the subtraction tensors (3.17) $\boldsymbol {\varPi }_1(\kern0.09em \boldsymbol {y})$ and (3.18) $\boldsymbol {\varPi }_2(\kern0.09em \boldsymbol {y})$ mostly follows Appendix B of Zinchenko & Davis (Reference Zinchenko and Davis2005). However, with our $e_{nf}=0.1$ and $\tilde {e}_{nf}=10$, the strength factors in their (B6) were increased 10 and 40 times for $x_3>0$ and $x_3<0$ patches, respectively.

D.4. Validation and performance of multipole acceleration

As an additional check of correctness and efficiency of our multipole-accelerated BI algorithm, the inhomogeneous term (3.2) $\boldsymbol {F}(\kern0.09em \boldsymbol {y})$ and iterative solution of the BI equation (3.1), corresponding to the final drop–wall configuration with $\lambda =60$, $\theta =7.5^\circ$, $B=0.125$ and $N_\triangle =138$ K in figure 9(a), were computed differently by simple node-to-node summations in the regularized integrals (3.8), (3.10) and (3.11), without any multipole operations. Two metrics were used to quantify deviations of such a direct-summation solution $(\boldsymbol {F}_{ex},\boldsymbol {u}_{ex})$ from our multipole-accelerated solutions $(\boldsymbol {F},\boldsymbol {u})$ for the same configuration, controlled by the intuitive precision parameter $\varepsilon$. Namely, for the inhomogeneous term (3.2)

(D4)$$\begin{gather} \delta_1(\boldsymbol{F},\boldsymbol{F}_{ex})= \frac{1}{\langle\boldsymbol{F}_{ex}^2\rangle^{1/2}}\max_{\boldsymbol{x}^i\in S}\|\boldsymbol{F}(\boldsymbol{x}^i)-\boldsymbol{F}_{ex}(\boldsymbol{x}^i)\|, \end{gather}$$
(D5)$$\begin{gather}\delta_2(\boldsymbol{F},\boldsymbol{F}_{ex})= \frac{1}{\langle\boldsymbol{F}_{ex}^2\rangle^{1/2}} [\langle(\boldsymbol{F}-\boldsymbol{F}_{ex})^2\rangle]^{1/2}, \end{gather}$$

where the angular brackets stand for averaging over the whole drop surface $S$; $\delta _1$ quantifies maximum pointwise deviation between $\boldsymbol {F}$ and $\boldsymbol {F}_{ex}$ in the mesh nodes relative to root-mean-square $\boldsymbol {F}_{ex}$, while $\delta _2$ is the integral measure of such deviation. The deviations $\delta _1(\boldsymbol {u},\boldsymbol {u}_{ex})$ and $\delta _2(\boldsymbol {u},\boldsymbol {u}_{ex})$ are defined the same way through $\boldsymbol {u}$ and $\boldsymbol {u}_{ex}$; however, for test purposes below in the multipole-accelerated solution of BI equation (3.1), the inhomogeneous term $\boldsymbol {F}$ is set to $\boldsymbol {F}_{ex}$ to isolate the effect of multipole truncation errors on the double-layer integral calculation. The initial approximation for the iterative solution $\boldsymbol {u}_{ex}$ was taken from the prelast time step of the multipole-accelerated dynamical simulation (with $\varepsilon = 10^{-6}$), followed by 50 GMRES iterations to make $\boldsymbol {u}_{ex}$ fully independent of the initial approximation. For consistency, the same initial approximation and 50 iterations were also used below in multipole-accelerated calculations of $\boldsymbol {u}$ at various values of the intuitive precision parameter $\varepsilon$.

Table 4 demonstrates that $\delta _i(\boldsymbol {F},\boldsymbol {F}_{ex})$ and $\delta _i(\boldsymbol {u},\boldsymbol {u}_{ex})$ tend to zero, as $\varepsilon \to 0$, which proves convergence of our code to the direct-summation code. The CPU times in table 4 are for single-core calculations on a personal Linux workstation with 3.6 GHz clock speed. As $\varepsilon$ tightens, the multipole-accelerated solution $(\boldsymbol {F},\boldsymbol {u})$ becomes practically indistinguishable from $(\boldsymbol {F}_{ex},\boldsymbol {u}_{ex})$ and is still much faster than by the direct summations. Somewhat surprisingly, values of $\varepsilon \sim 10^{-5}{-}10^{-6}$ were always sufficient in the present dynamical simulations, as further tightening of $\varepsilon$ did not affect four significant figures of the steady-state drop speed (and could only marginally affect the minimum drop–wall clearance in rare cases); an appreciable maximum error $\delta _1(\boldsymbol {F},\boldsymbol {F}_{ex})$ appears not to have a significant effect, since the root-mean-square error $\delta _2(\boldsymbol {F},\boldsymbol {F}_{ex})$ is two orders of magnitude less. For $\varepsilon \sim 10^{-5}{-}10^{-6}$ and resolution $N_\triangle =138$ K in table 4, multipole acceleration speeds up $\boldsymbol {F}$-calculation by a factor of 23–27, and makes BI iterations 25–35 times faster, compared with the direct-summation algorithm. Precalculation of matrices (3.17), (3.18) and $\boldsymbol {\varGamma }$ before the iterations adds insignificantly to the cost of the multipole-accelerated solution.

Table 4. Convergence of the multipole-accelerated solution $(\boldsymbol {F},\boldsymbol {u})$ to the direct-summation solution $(\boldsymbol {F}_{ex},\boldsymbol {u}_{ex})$, as the precision parameter $\varepsilon \to 0$, for one steady drop configuration with $\lambda =60$, $\theta =7.5^\circ$, $B=0.125$ and $N_\triangle =138$ K. The CPU times are in seconds on a single core. The direct-summation solution requires 67 s for the inhomogeneous term $\boldsymbol {F}_{ex}$ and 64 s per iteration for $\boldsymbol {u}_{ex}$.

A similar analysis was performed for a strongly pancaked drop with a large contact spot and higher resolution $N_\triangle =246$ K, corresponding to the final, steady-state drop–wall configuration in the run $\lambda =60$, $\theta =7.5^\circ$, $B=5$ of figure 9(d). Again, tightening $\varepsilon$ can make the multipole-accelerated solution arbitrarily close to the result by the simplest node-to-node summations (table 5). Compared with the above case of small drop deformation, even more dramatic advantage of the multipole-accelerated solution is observed. For $\varepsilon \sim 10^{-5}{-}10^{-6}$ sufficient for the present dynamical simulations, our algorithm calculates the inhomogeneous term $\boldsymbol {F}$ 34–40 times faster, and performs the BI iteration 31–39 times faster than by direct node-to-node summations.

Table 5. Convergence of the multipole-accelerated solution $(\boldsymbol {F},\boldsymbol {u})$ to the direct-summation solution $(\boldsymbol {F}_{ex},\boldsymbol {u}_{ex})$, as $\varepsilon \to 0$, for one steady drop configuration with $\lambda =60$, $\theta =7.5^\circ$, $B=5$ and $N_\triangle =246$ K. The direct-summation solution requires 224 s for the inhomogeneous term $\boldsymbol {F}_{ex}$ and 197 s per iteration for $\boldsymbol {u}_{ex}$, to be compared with the CPU times (in seconds) for multipole-accelerated solution at various $\varepsilon$.

In most of the runs, the parallel version of the multipole-accelerated code (achieved through OMP programming) was used, with an additional boost in performance; for example, the speed up due to parallelization was 4.5-fold on eight cores of a personal Linux workstation with 3.6 GHz clock speed. Parallelization of the direct-summation code only gave a similar boost, so it was still far from possible to use such a code in high-resolution dynamical simulations of the present work.

Despite crucial speed up in the BI solution owing primarily to multipole acceleration and (to a lesser extent) OMP programming, the dynamical high-resolution simulations are still extensive because of a large number of time steps to reach steady state. It took 260 and 230 K steps to generate curve 3 in figure 9(a) and curve 1 in figure 9(d), respectively.

References

Aussillous, P. 2002 Les gouttes enrobées. Doctoral thesis, Université Paris VI.Google Scholar
Aussillous, P. & Quéré, D. 2001 Liquid marbles. Nature 411, 924927.CrossRefGoogle ScholarPubMed
Aussillous, P. & Quéré, D. 2002 Bubbles creeping in a viscous liquid along a slightly inclined wall. Europhys. Lett. 59, 370376.CrossRefGoogle Scholar
Bazhlekov, I.B., Anderson, P.D. & Meijer, H.E.H. 2004 Nonsingular boundary integral method for deformable drops in viscous flows. Phys. Fluids 16, 10641081.CrossRefGoogle Scholar
Blake, J.R. 1971 A note on the image system for a stokeslet in a no-slip boundary. Proc. Camb. Phil. Soc. 70, 303310.CrossRefGoogle Scholar
Del Castillo, L.A., Ohnishi, S., White, L.R., Carnie, S.L. & Horn, R.G. 2011 Effect of disjoining pressure on terminal velocity of a bubble sliding along an inclined wall. J. Colloid Interface Sci. 364, 505511.CrossRefGoogle ScholarPubMed
Claeys, I.L. & Brady, J.F. 1989 Lubrication singularities of the grand resistance tensor for 2 arbitrary particles. Physicochem. Hydrodyn. 11, 261293.Google Scholar
Cox, R.G. 1974 The motion of suspended particles almost in contact. Intl J. Multiphase Flow 1, 343371.CrossRefGoogle Scholar
Dussan, V.E.B. 1985 On the ability of drops or bubbles to stick to non-horizontal surfaces of solids. Part 2. Small drops or bubbles having contact angles of arbitrary size. J. Fluid Mech. 151, 120.CrossRefGoogle Scholar
Dussan, V.E.B. & Chow, R.T.-P. 1983 On the ability of drops or bubbles to stick to non-horizontal surfaces of solids. J. Fluid Mech. 137, 129.CrossRefGoogle Scholar
Farutin, A., Biben, T. & Misbah, C. 2014 3D numerical simulations of vesicle and inextensible capsule dynamics. J. Comput. Phys. 275, 539568.CrossRefGoogle Scholar
Gissinger, J.R., Zinchenko, A.Z. & Davis, R.H. 2021 Drop squeezing between arbitrary smooth obstacles. J. Fluid Mech. 908, A33.CrossRefGoogle Scholar
Griggs, A.J., Zinchenko, A.Z. & Davis, R.H. 2008 Gravity-driven motion of a deformable drop or bubble near an inclined plane at low Reynolds number. Intl J. Multiphase Flow 34, 408418.CrossRefGoogle Scholar
Griggs, A.J., Zinchenko, A.Z. & Davis, R.H. 2009 Creeping motion and pending breakup of drops and bubbles near an inclined wall. Phys. Fluids 21, 093303.CrossRefGoogle Scholar
Happel, J. & Brenner, H. 1986 Low Reynolds Number Hydrodynamics. Martinus Nijhoff.Google Scholar
Hodges, S.R., Jensen, O.E. & Rallison, J.M. 2004 Sliding, slipping, and rolling: the sedimentation of a viscous drop down a gently inclined plane. J. Fluid Mech. 512, 95131.CrossRefGoogle Scholar
Kim, S. & Karilla, S. 1991 Microhydrodynamics: Principles and Selected Applications. Butterworth-Heinemann.Google Scholar
Klaseboer, E., Sun, Q. & Chan, D.Y.C. 2012 Non-singular boundary integral methods for fluid mechanics applications. J. Fluid Mech. 696, 468478.CrossRefGoogle Scholar
Magnaudet, J., Takagi, S. & Legendre, D. 2003 Drag, deformation and lateral migration of a buoyant drop moving near a wall. J. Fluid Mech. 476, 115157.CrossRefGoogle Scholar
Mahadevan, L. & Pomeau, Y. 1999 Rolling drops. Phys. Fluids 11, 24492453.CrossRefGoogle Scholar
Meyyappan, M. & Subramanian, R.S. 1987 Thermocapillary migration of a gas bubble in an arbitrary direction with respect to a plane surface. J. Colloid Interface Sci. 115, 206219.CrossRefGoogle Scholar
O'Neill, M.E. & Stewartson, K. 1967 On the slow motion of a sphere parallel to a nearby plane wall. J. Fluid Mech. 27, 705724.CrossRefGoogle Scholar
Pozrikidis, C. 1992 Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press.CrossRefGoogle Scholar
Quéré, D. 2005 Non-sticking drops. Rep. Prog. Phys. 68, 24952532.CrossRefGoogle Scholar
Rahman, M.R. & Waghmare, P.R. 2018 Influence of outer medium viscosity on the motion of rolling droplets down an incline. Phys. Rev. Fluids 3, 023601.CrossRefGoogle Scholar
Rallison, J.M. & Acrivos, A. 1978 A numerical study of the deformation and burst of a viscous drop in an extensional flow. J. Fluid Mech. 89, 191200.CrossRefGoogle Scholar
Richard, D. & Quéré, D. 1999 Viscous drops rolling on a tilted non-wettable solid. Europhys. Lett. 48, 286291.CrossRefGoogle Scholar
Sangani, A.S. & Mo, G. 1996 An $O(N)$ algorithm for Stokes and Laplace interactions of particles. Phys. Fluids 8, 19902010.CrossRefGoogle Scholar
Schnitzer, O., Davis, A.M.J. & Yariv, E. 2020 Rolling of non-wetting droplets down a gently inclined plane. J. Fluid Mech. 903, A25.CrossRefGoogle Scholar
Staben, M.E., Zinchenko, A.Z. & Davis, R.H. 2006 Dynamic simulation of spheroid motion between two parallel walls in low-Reynolds-number Poiseuille flow. J. Fluid Mech. 553, 187226.CrossRefGoogle Scholar
Zinchenko, A.Z. 1980 The slow asymmetric motion of two drops in a viscous medium. J. Appl. Math. Mech. 44, 3037.CrossRefGoogle Scholar
Zinchenko, A.Z. 1994 An efficient algorithm for calculating multiparticle thermal interaction in a concentrated dispersion of spheres. J. Comput. Phys. 111, 120135.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2000 An efficient algorithm for hydrodynamical interaction of many deformable drops. J. Comput. Phys. 157, 539587.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2002 Shear flow of highly concentrated emulsions of deformable drops by numerical simulations. J. Fluid Mech. 455, 2162.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2005 A multipole-accelerated algorithm for close interaction of slightly deformable drops. J. Comput. Phys. 207, 695735.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2006 A boundary-integral study of a drop squeezing through interparticle constrictions. J. Fluid Mech. 564, 227266.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2008 Algorithm for direct numerical simulation of emulsion flow through a granular material. J. Comput. Phys. 227, 78417888.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2013 Emulsion flow through a packed bed with multiple drop breakup. J. Fluid Mech. 725, 611663.CrossRefGoogle Scholar
Zinchenko, A.Z. & Davis, R.H. 2017 General rheology of highly concentrated emulsions with insoluble surfactant. J. Fluid Mech. 816, 661704.CrossRefGoogle Scholar
Zinchenko, A.Z., Rother, M.A. & Davis, R.H. 1997 A novel boundary-integral algorithm for viscous interaction of deformable drops. Phys. Fluids 9, 14931511.CrossRefGoogle Scholar
Figure 0

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

Figure 1

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 2

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.

Figure 3

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.

Figure 4

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.

Figure 5

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.

Figure 6

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.

Figure 7

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).

Figure 8

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.

Figure 9

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.

Figure 10

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.

Figure 11

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).

Figure 12

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).

Figure 13

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

Figure 14

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 15

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).

Figure 16

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$.

Figure 17

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.

Figure 18

Table 4. Convergence of the multipole-accelerated solution $(\boldsymbol {F},\boldsymbol {u})$ to the direct-summation solution $(\boldsymbol {F}_{ex},\boldsymbol {u}_{ex})$, as the precision parameter $\varepsilon \to 0$, for one steady drop configuration with $\lambda =60$, $\theta =7.5^\circ$, $B=0.125$ and $N_\triangle =138$ K. The CPU times are in seconds on a single core. The direct-summation solution requires 67 s for the inhomogeneous term $\boldsymbol {F}_{ex}$ and 64 s per iteration for $\boldsymbol {u}_{ex}$.

Figure 19

Table 5. Convergence of the multipole-accelerated solution $(\boldsymbol {F},\boldsymbol {u})$ to the direct-summation solution $(\boldsymbol {F}_{ex},\boldsymbol {u}_{ex})$, as $\varepsilon \to 0$, for one steady drop configuration with $\lambda =60$, $\theta =7.5^\circ$, $B=5$ and $N_\triangle =246$ K. The direct-summation solution requires 224 s for the inhomogeneous term $\boldsymbol {F}_{ex}$ and 197 s per iteration for $\boldsymbol {u}_{ex}$, to be compared with the CPU times (in seconds) for multipole-accelerated solution at various $\varepsilon$.

Supplementary material: File

Zinchenko supplementary material

Zinchenko supplementary material
Download Zinchenko supplementary material(File)
File 35.6 KB