4.1. General desingularization formulae

A repository of desingularization techniques that permit stable hydrodynamic interactions between drops and solid particles is provided by Zinchenko & Davis (Reference Zinchenko and Davis2006), including drop and solid self-interactions, solid–solid interaction, solid–drop and drop–solid contributions. In certain cases (solid self-interaction, solid–solid interaction and solid–drop contributions), additional integrals appear during the subtraction procedure that must be calculated with near-analytical accuracy to facilitate tight-squeezing simulations. The complete desingularization expressions and relevant ‘high-accuracy integrals’ in each case are listed below.

For solid self-interaction ($\boldsymbol {y}\in \hat {S}_\alpha$), singularity subtraction is provided by

(4.1)\begin{align} & \int_{\hat{S}_\alpha} \boldsymbol{q}(\boldsymbol{x})\boldsymbol{\cdot}[2\boldsymbol{\tau} (\boldsymbol{r})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})+ \eta \boldsymbol{G}(\boldsymbol{r})] \,\mbox{d} S_x \nonumber\\ &\quad =\boldsymbol{q}(\boldsymbol{y}) + \int_{\hat{S}_\alpha} [\boldsymbol{q}(\boldsymbol{x})-\boldsymbol{q}(\boldsymbol{y})] \boldsymbol{\cdot}[2\boldsymbol{\tau}(\boldsymbol{r}) \boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})+\eta \boldsymbol{G}(\boldsymbol{r})] \, \mbox{d} S_x + \eta\boldsymbol{q}(\boldsymbol{y})\int_{\hat{S}_\alpha} \boldsymbol{G}(\boldsymbol{r}) \,\mbox{d} S_x, \end{align}

where the final integral on the right-hand side of (4.1) must be computed with high accuracy.

Solid–solid interactions are transformed slightly differently, but contain the same integral requiring near-analytical treatment

(4.2)\begin{align} & \int_{\hat{S}_\beta} \boldsymbol{q}(\boldsymbol{x})\boldsymbol{\cdot}[2\boldsymbol{\tau} (\boldsymbol{r})\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})+ \eta \boldsymbol{G}(\boldsymbol{r})] \,\mbox{d} S_x \nonumber\\ &\quad =\int_{\hat{S}_\beta} [\boldsymbol{q}(\boldsymbol{x})-\boldsymbol{q}(\boldsymbol{x}^{*})] \boldsymbol{\cdot}[2\boldsymbol{\tau}(\boldsymbol{r})\boldsymbol{\cdot}\boldsymbol{n} (\boldsymbol{x})+\eta \boldsymbol{G}(\boldsymbol{r})] \,\mbox{d} S_x + \eta\boldsymbol{q}(\boldsymbol{x}^{*})\int_{\hat{S}_\beta}\boldsymbol{G}(\boldsymbol{r}) \,\mbox{d} S_x, \end{align}

where $\boldsymbol {x}^{*} = \boldsymbol {x}^{*}(\boldsymbol {y},\beta )$ is the mesh node on $\hat {S}_\beta$ closest to $\boldsymbol {y}\in \hat {S}_\alpha \neq \hat {S}_\beta$.

Likewise, the single-layer part of the solid–drop contribution utilizes a simple regularization also involving $\boldsymbol {G}(\boldsymbol {r})$

(4.3)\begin{equation} \int_{\hat{S}_\beta} \boldsymbol{q}(\boldsymbol{x})\boldsymbol{\cdot} \boldsymbol{G}(\boldsymbol{r})\,\mbox{d} S_x = \int_{\hat{S}_\beta} [\boldsymbol{q}(\boldsymbol{x})-\boldsymbol{q}(\boldsymbol{x}^{*})] \boldsymbol{\cdot}\boldsymbol{G}(\boldsymbol{r})\,\mbox{d} S_x + \boldsymbol{q}(\boldsymbol{x}^{*}) \int_{\hat{S}_\beta} \boldsymbol{G}(\boldsymbol{r})\,\mbox{d} S_x. \end{equation}

Finally, the high-order near-singularity subtraction used for the solid–drop double-layer contribution was found to be especially critical for stable tight-squeezing simulations. To this end, for $\boldsymbol {y}$ on a drop surface $\tilde {S}_{\alpha }$ and its nearest mesh node $\boldsymbol {x}^{*}$ on a solid surface $\hat {S}_{\beta }$, the density function $\boldsymbol {q}(\boldsymbol {x})$ near $\boldsymbol {x}^{*}$ is approximated by $\boldsymbol {q}(\boldsymbol {x^{*}})$ plus a linear function $\mathscr {L}(\boldsymbol {x}-\boldsymbol {x}^{*})$; the coefficients of the linear form $\mathscr {L}$ are found by least-squares fitting to the values of ${\rm \Delta} \boldsymbol {q}(\boldsymbol {x}^{j})=\boldsymbol {q}(\boldsymbol {x}^{j})-\boldsymbol {q}(\boldsymbol {x}^{*})$ in the solid mesh nodes $\boldsymbol {x}^{j}$ directly connected to $\boldsymbol {x}^{*}$. Subtracting $\boldsymbol {q}(\boldsymbol {x}^{*})+\mathscr {L}(\boldsymbol {x}-\boldsymbol {x}^{*})$ from $\boldsymbol {q}(\boldsymbol {x})$ completely desingularizes the double-layer part of the solid–drop contribution in (3.1), but gives rise to added-back integrals, which need to be handled either analytically (possible only for a few shapes), or, in general, numerically with high accuracy. This desingularization can be written as

(4.4)\begin{align} & \int_{\hat{S}_\beta} \boldsymbol{q}(\boldsymbol{x})\boldsymbol{\cdot}\boldsymbol{\tau}(\boldsymbol{r}) \boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x}) \,\mbox{d} S_x \nonumber\\ &\quad =\int_{\hat{S}_\beta}\left[ \boldsymbol{q}(\boldsymbol{x})-\boldsymbol{q}(\boldsymbol{x}^{*})- \sum_{j\in\mathcal{A}(\boldsymbol{x}^{*})} (\boldsymbol{c}_j \boldsymbol{\cdot}(\boldsymbol{x}-\boldsymbol{x}^{*})){\rm \Delta}\boldsymbol{q}(\boldsymbol{x}^{j})\right] \boldsymbol{\cdot} \boldsymbol{\tau}(\boldsymbol{r}) \boldsymbol{\cdot} \boldsymbol{n}(\boldsymbol{x}) \,\mbox{d} S_x \nonumber\\ &\qquad +\sum_{j\in\mathcal{A}(\boldsymbol{x}^{*})} [\boldsymbol{c}_j{\rm \Delta}\boldsymbol{q}(\boldsymbol{x}^{j})]\boldsymbol{\cdot}\frac{3}{4{\rm \pi}} \int_{\hat{S}_\beta} \frac{[\boldsymbol{r}\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})] \boldsymbol{r}\boldsymbol{r}\boldsymbol{r}}{r^{5}}\,\mbox{d} S_x, \end{align}

where $\mathcal {A}(\boldsymbol {x}^{*})$ is the set of nodes directly connected to $\boldsymbol {x}^{*}$; the vector coefficients $\boldsymbol {c}_j$ depend only on the local mesh geometry near $\boldsymbol {x}^{*}$ and are pre-calculated as detailed in § 3 of Zinchenko & Davis (Reference Zinchenko and Davis2006).

So, high-accuracy computation of three integral terms is required. Up to a factor, we denote these terms as

(4.5a,b)\begin{gather} G_1(\boldsymbol{y},\beta) \equiv \int_{\hat{S}_\beta}\frac{1}{r}\,\mbox{d} S_x,\quad \boldsymbol{G}_2(\boldsymbol{y},\beta) \equiv \int_{\hat{S}_\beta}\frac{\boldsymbol{r}\boldsymbol{r}}{r^{3}}\,\mbox{d} S_x, \end{gather}

(4.6)\begin{gather}\boldsymbol{\varPsi}(\boldsymbol{y},\beta) \equiv \int_{\hat{S}_\beta}\frac{(\boldsymbol{r}\boldsymbol{\cdot}\boldsymbol{n}(\boldsymbol{x})) \boldsymbol{r}\boldsymbol{r}\boldsymbol{r}}{r^{5}}\,\mbox{d} S_x. \end{gather}

It is advantageous that the second-rank tensor $\boldsymbol {G}_2$ and the third-rank tensor $\boldsymbol {\varPsi }$ are fully permutable in their indices, thus reducing the amount of work to calculate these tensors.

4.2. Semi-analytical desingularization for axisymmetric particles

Axisymmetric particles permit semi-analytical calculation of the integrals (4.5a,b) and (4.6), by closed-form azimuthal integration around the particle axis of rotation, with external numerical integration along the particle half-contour. This kind of approach has been widely used in axisymmetric boundary-integral solutions (e.g. Rallison & Acrivos Reference Rallison and Acrivos1978; Lee & Leal Reference Lee and Leal1982; Davis Reference Davis1999; Ratcliffe et al. Reference Ratcliffe, Zinchenko and Davis2010) by reducing the azimuthal integration to elliptic integrals (as discussed in detail by Pozrikidis Reference Pozrikidis1992). In the present 3-D work, analytical azimuthal integration for the high-order subtraction tensor (4.6) is a more cumbersome, but still very efficient approach. As shown in figure 2, a representative integration circumference $C$ is the perimeter of a representative disk of radius $\rho$, a sufficient number of which form an approximation to the axisymmetric volume. The surface integral of a function $f(\boldsymbol {x})$ is simply

(4.7)\begin{equation} \int_{\hat{S}_\beta} f(\boldsymbol{x})\,\mbox{d} S_x = \int \rho\,\mbox{d}s \left[\oint_{C}f(\boldsymbol{x})\,\mbox{d} \varPsi \right], \end{equation}

where $\varPsi \in [-{\rm \pi} ,{\rm \pi} ]$ is the azimuthal angle of rotation around the particle axis, and $\mbox {d}s$ is the contour length element along the particle's half-contour.

Figure 2.

An axisymmetric particle oriented as a surface of revolution $\hat {S}_\beta$ around the $z$-axis. The single- and double-layer contributions of circular segments can be calculated analytically using elliptic integrals. For a given representative circle $C$, the origin $\boldsymbol {O}$ may be placed at its centre, with the $x$-axis in plane but oriented away from the observation point $\boldsymbol {x}_0$. Cylindrical coordinates ($\rho ,\psi ,z$) are used for deriving elliptic integral calculations.

Now, we outline derivations for closed-form expressions for the integrals (4.5a,b) and (4.6) over a representative circumference $C$, in terms of complete elliptic integrals of the first and second kind. Let the origin $\boldsymbol {O}$ of the temporary coordinate system $(\rho \cos {\psi },\rho \sin {\psi },z)$ be at the centre of $C$, with its axis of revolution along the $z$-axis (figure 2). Without loss of generality, the $x$-axis can be defined as antiparallel to the shortest vector from the $z$-axis to the observation point $\boldsymbol {x}_0$. The observation point is then $\boldsymbol {x}_0 = (-\rho _0,0,z_0)$ and the vector $\boldsymbol {r}=\boldsymbol {x}-\boldsymbol {x}_0$ between $\boldsymbol {x}_0$ and a integration point $\boldsymbol {x}$ on $C$ is given by

(4.8)\begin{equation} \boldsymbol{r} = (\rho\cos{\psi}+\rho_0,\rho\sin{\psi},z-z_0), \end{equation}

and

(4.9)\begin{align} r^{2} = \left[(\rho+\rho_0)^{2}+(z-z_0)^{2}\right]\left[1-k^{2}\sin^{2}{\frac{\psi}{2}}\right],\quad \mbox{where}\ k^{2} = \frac{4\rho\rho_0}{(\rho+\rho_0)^{2}+(z-z_0)^{2}}. \end{align}

The simplest contour integral over $C$ relevant to computing $G_1$ may be expressed as

(4.10)\begin{equation} \int_{-{\rm \pi}}^{\rm \pi} \frac{1}{r} \mbox{d} \psi = c_0\int_{-{\rm \pi}}^{\rm \pi} \left[1-k^{2}\sin^{2}{\frac{\psi}{2}}\right]^{{-}1/2}\mbox{d} \psi = 4c_0F, \end{equation}

where $c_0 = ((\rho +\rho _0)^{2}+(z-z_0)^{2})^{-1/2}$ and $F$ is the complete elliptic integral of the first kind with modulus $k$. The relevant integrals over $C$ for $\boldsymbol {G}_2$ and $\boldsymbol {\varPsi }$ can be derived in a similar manner and final expressions are provided in appendix A. Once $\boldsymbol {G}_2$ and $\boldsymbol {\varPsi }$ are computed in the special coordinate basis defined above, they are then transformed by usual tensor transformation rules to the global basis used in the squeezing simulation.

The present semi-analytical approach results in very fast and accurate calculation of single- and double-layer desingularization integrals (4.5a,b) and (4.6) for axisymmetric particles. Surface integrals calculated by (4.7) converge quickly with respect to (with second-order Euler external integration). For example, using a capsule ($\mathcal {L}=6$) as a test particle and an observation point $0.01L$ away from the cylindrical section of particle surface, all $\boldsymbol {\varPsi }$ terms agree to within at least the fifth decimal point for values of 1000 and 10 000. The former corresponds to a representative disk width of 0.006. Further convergence tests and comparisons with the multimesh method are provided in § 5.1.

4.3. Multimesh desingularization

Desingularization for arbitrarily shaped particles can be achieved with direct numerical calculation of $G_1$, $\boldsymbol {G}_2$ and $\boldsymbol {\varPsi }$ (4.5a,b) and (4.6), if sufficiently high resolution meshes are used to discretize these surface integrals. Note that this ‘multimesh desingularization method’ still requires particles to have smoothly varying normals (i.e. to be Lyapunov surfaces). Extremely high resolutions are required for single- and double-layer contributions to approach exact values, in order to accurately and stably reproduce analytical or semi-analytical tight-squeezing results. Due to the singular nature of the integrals, accuracy is more sensitive to mesh resolution as the observation point approaches the particle surface; therefore, a hierarchy of embedded meshes is used for multimesh desingularization.

Two high-resolution auxiliary meshes are used to compute the desingularization integrals. The first high-resolution mesh typically has four times more triangles than the basic mesh used for BIM calculations. The second is an ultra-high-resolution mesh, typically 64 times denser than the basic mesh, as shown in figure 3(a,b). For particles with higher surface area, the ultra-high-resolution mesh may have in excess of 2.5 million triangles. Therefore, the two meshes are chained to each other: the indices of the 16 ultra-high-resolution triangles within each high-resolution triangle are pre-calculated and stored. In contrast to boundary-integral calculations for desingularized terms on the basic mesh, with mesh triangle contribution reassigned to vertices (a procedure due to Rallison (Reference Rallison1981) used in many prior simulations), such reassignment is disadvantageous for multimesh calculations of the desingularization tensors. A primitive form is used instead:

(4.11)\begin{equation} \int_{\hat{S}_\beta} f(\boldsymbol{x}) \,\mbox{d}S_x \approx \sum_i f(\boldsymbol{x}_i){\rm \Delta} S_i, \end{equation}

where the summation is over all triangles for a given triangulation of the surface $\hat {S}_\beta$, $\boldsymbol {x}_i$ is the triangle centroid and ${\rm \Delta} S_i$ is the flat triangle area. This form makes it straightforward to use one mesh within a given cutoff distance $r_c$ from the observation point, and another one beyond that radius, as illustrated in figure 3(c). Specifically, the distance between a given observation point $\boldsymbol {x}_0$ and high-resolution triangle is determined using the triangle centroid. If this distance is greater than the prescribed distance cutoff $r_c$, then this triangle is used as an area element ${\rm \Delta} S_i$ of the surface integral in (4.11). Otherwise, the 16 corresponding ultra-high-resolution triangles are substituted to represent this triangle contribution to (4.11). The desingularization integrals for solid self-interaction and solid–solid interactions are calculated only once at the start of each simulation, so a high distance cutoff ($r_c=1.0L$) was used with a negligible penalty on total CPU time. The solid–drop contribution to the desingularization integrals is computed outside of GMRES iterations, but still represents the slowest operation of the multimesh method. A distance cutoff of $r_c=0.4L$ for the solid–drop contribution was found sufficient to stably model tight-squeezing or trapped drops with an accuracy approaching that of analytical methods, as discussed in § 5.1.

Figure 3.

Solid-particle meshes used for the $\mathcal {L}=4$ capsule by the multimesh desingularization method. (a) The high-resolution ($\hat {N}_{\triangle }=52\ \text {K}$) and ultra-high-resolution ($\hat {N}_{\triangle }=840\ \text {K}$) meshes used for numerical desingularization. The additional edges of the ultra-high-resolution mesh are coloured cyan. (b) The basic mesh ($\hat {N}_{\triangle }=13\ \text {K}$) used by the BIM simulation for desingularized terms. (c) A given observation point $\boldsymbol {x}_0$ calculates per-triangle values from the ultra-high-resolution mesh within a cutoff $r_c$.

Several technical issues arise due to the requirements of (i) quality meshing of arbitrary shapes, (ii) obtaining machine precision for node positions, (iii) memory requirements of ultra-high-resolution meshes and (iv) mesh chaining. An approximate mesh for each particle resolution was created using the powerful open-source software Blender, which has no limitation on mesh size and can easily handle meshes with $N_{\triangle }>2.5\ \mbox {M}$, but uses single-precision floating-point format for coordinates (which may create unwanted numerical errors in tight-squeezing simulations). After node positions and connectivity are exported from Blender, coordinates are read into the BIM program and adjusted to double precision according to piecewise analytical functions that define a given particle shape. Similarly, the surface normals at each node are calculated using analytical formulae for a particular surface. A more general approach would be to use just node positions on the solid surfaces and calculate the normals by the same tools as we employ for the triangulated drop surface (appendix C), but such generality was not pursued in the present calculations. Due to memory considerations, only one instance of each desingularization mesh is saved in memory, and is translated and/or mirrored (see § 4.4) during calculations over a given particle. Finally, an efficient algorithm was written to chain two meshes to each other, i.e. to determine which ultra-high-resolution triangles lie within each high-resolution triangle. The algorithm requires knowledge only of the node positions and connectivity of the two meshes, assuming every high-resolution node coincides with an ultra-high-resolution node. Briefly, coinciding nodes between the meshes are determined by position, and a nearby ultra-high-resolution triangle is associated with a given high-resolution triangle if its centroid (projected along the surface normal) lies within the perimeter of the high-resolution triangle.

4.4. Additional numerical considerations

Two configurations, drop flow around a single capsule and between two flat plates, presented particular challenges regarding numerical stability and drop mesh quality. A useful extension of the commonly used triangle subdivision scheme introduced by Cristini, Bławzdziewicz & Loewenberg (Reference Cristini, Bławzdziewicz and Loewenberg2001) was developed, which is particularly effective when remeshing contiguous regions on a drop. The simple idea is, given a list of all triangles to be subdivided within a given region, a uniform mesh four times denser than the original can immediately result, by adding nodes at all edge midpoints and appropriately connecting them (see appendix B). Connectivity around the perimeter of this region can also be determined simply, e.g. in the manner of Cristini et al. (Reference Cristini, Bławzdziewicz and Loewenberg2001) (§ A.1.2. therein). Occasional node addition on regions of contiguous triangles had two benefits. First, a much cleaner mesh immediately resulted, such that further relaxation tangential to the drop surface was unneeded. Second, updating the mesh topology was required much less often, resulting in a more efficient algorithm. For example, in the case of a single drop approaching a single fibre, passive mesh stabilization alone (in the form of Zinchenko & Davis Reference Zinchenko and Davis2006) was found to be inadequate due to the large drop deformation. This simulation was greatly extended by defining an upper limit on triangle area, and whenever this limit was reached, executing node addition on all triangles within a cutoff distance of the triangle with maximal area. In this case, similar simulation times were reached by actively redistributing nodes using an advanced potential energy function (Zinchenko & Davis Reference Zinchenko and Davis2013). A fourth-order variation of the best-paraboloid method (Zinchenko & Davis Reference Zinchenko and Davis2006) was used to used to choose coordinates for the inserted nodes. Including the second- and third-level neighbours in the fit improved the local curvatures at inserted nodes, in areas of both low and high curvature, and allowed the simulation to continue stably without the intermediate step of relaxing inserted nodes. Note that a four-times-denser triangle mesh corresponds to mesh edges one half their original length; the tighter Courant stability limitation on the time step was not prohibitive for any systems presented herein. To maintain mesh quality on the whole drop, passive mesh stabilization (Zinchenko & Davis Reference Zinchenko and Davis2006) is also enforced for all nodes at every time step.

The above ‘coordinated adaptive remeshing’ scheme was extended to allow added nodes to be subsequently removed. In the case of a drop squeezing between thin plates, especially at low capillary numbers, it was necessary to include a spatially dependent remeshing scheme for the drop surface, with criteria defined to increase mesh resolution within near-contact regions of the drop (figure 4a,b). In this case, the drop does not experience particularly high deformation or curvatures, yet mesh stability greatly benefits from increased resolution in specific regions. This criteria-based local remeshing is achieved by remembering the original mesh connectivity, restoring it and re-evaluating the given criteria for refinement (every so many time steps). Therefore, the region of higher mesh resolution shown in figure 4 remains constant as the drop moves through the constriction. This technique was also necessary for certain half-capsule configurations at low capillary number.

Figure 4.

Adaptive mesh schemes for droplet squeezing between plates. (a) Coordinated adaptive meshing of drop surface based on spatial criteria. By storing the original mesh connectivity and re-evaluating the criteria before every remeshing, the region of denser triangulation remains constant as the drop moves through the constriction. (b) Orthogonal view of the same drop, with solid particles hidden.

To prevent symmetry breaking during drop flow, it was necessary that individual solid-particle meshes were made symmetric with respect to all three orthogonal planes that pass through the constriction centre. Also, for the flat-plate system, a fixed adaptive mesh was used for solid particles to increase node density in near-contact regions. Mesh density within two unit lengths of the constriction centre was increased between four and sixteen times, depending on proximity to the constriction centre. This base mesh was refined to create the higher-resolution desingularization meshes. In order to maintain the aforementioned symmetry, while also storing only one copy of the higher-resolution meshes in memory, a value is associated with each particle to indicate whether or not to use the mirror image of a given mesh. Inverting a coordinate is not computationally intensive, allowing for symmetric systems with non-uniform meshes to be efficiently treated within the framework of the multimesh desingularization method.

5.1. Convergence with mesh resolution and desingularization methods

Convergence of results for the default capsule length ($\mathcal {L} = 6$) with respect to solid-particle mesh resolution is shown in figure 5. The $\hat {N}_{\triangle } = 35.6\ \text {K}$ mesh utilized adaptive resolution such that the mesh was four times denser within a distance of ${\approx }1.4$ from the constriction centre, compared to the $\hat {N}_{\triangle } = 20.2\ \text {K}$ mesh. Although was set throughout this work for the number of characteristic disks used in the semi-analytical method (increased for longer capsules to provide an equivalent linear density), figure 5 also shows that half that number produces reasonable results (resulting in an eight-times speedup for desingularization as compared to the multimesh method). Indeed, halving the number of characteristic disks () results in agreement near the solution tolerance ($1\times 10^{-5}$) for the majority of the squeezing process.

The multimesh implementation was validated against previous results that used analytical desingularization (Zinchenko & Davis Reference Zinchenko and Davis2006), and the multimesh and semi-analytical methods were compared against each other. As shown in figure 6, good agreement is seen for drops squeezing between spheres and spheroidal particles. These two-particle constrictions are prone to symmetry breaking and sensitive to mesh resolutions, so reproducing these tight-squeezing results efficiently and without analytical expressions was unexpected. In particular, the spheroidal system in figure 6(b) is metastable and undergoes symmetry breaking shortly after the minimum drop velocity, with one side of the drop advancing slightly farther downstream than the other, even for the analytical system. For these tests, each higher-resolution mesh was 16 times denser than the previous mesh in the hierarchy. Further increasing these desingularization-mesh densities negatively impacted performance while negligibly affecting droplet dynamics. As shown in figure 9, close agreement was also observed between the multimesh and semi-analytical desingularization techniques for drop squeezing between axisymmetric particles, even for tight squeezing.

Figure 6.

Comparison between two types of desingularization methods ($\tilde {N}_{\triangle }=8.6\ \text {K}$, $\hat {N}_{\triangle }=5\ \text {K}\times 16\times 16$). Insets shows snapshot of simulation at $U_{min}$. (a) Multimesh and analytical desingularization methods agree well for drop squeezing between two spheres ($\epsilon =0.5,\tilde {a}=0.9,Ca=0.63,\lambda =4.0$). (b) Drop squeezing between two parallel oblate spheroids with major and minor axes of 1 and 0.4, respectively ($\epsilon =0.5,\tilde {a}=0.75,Ca=0.4,\lambda =1.0$). (c) Absolute error for various methods of calculating the desingularization tensor $\boldsymbol {\varPsi }$, as compared to the analytical values for the oblate spheroid shown in (b). The semi-analytical method uses 3000 characteristic disks with adaptive resolution for the disk spacing.

The absolute error of the semi-analytical and multimesh methods for $\boldsymbol {\varPsi }(\boldsymbol {y})$-calculation is shown in figure 6(c) for an oblate spheroid (one of those in the insert of figure 6b), with the error of the basic mesh (without desingularization) shown for comparison. This error is quantified through root-mean-square deviation of the components $\varPsi _{ijk}$ from the analytical values available for this shape (Zinchenko & Davis Reference Zinchenko and Davis2006). Unlike with the cylindrical particles discussed below, $\boldsymbol {\varPsi }$ is calculated for observation points $\boldsymbol {y}$ near characteristic disks with radii approaching zero, where the characteristic disk radius is indicated by the vertical height of the spheroid in figure 6(c). It was found that a higher density of characteristic disks is required when the radius is near zero. For example, when using 10 000 equally spaced disks for the oblate spheroid (disk width of $8\times 10^{-5}$), the semi-analytical method was conspicuously less accurate than the multimesh method, especially near the horizontal axis of figure 6(c). It was necessary to use an adaptive density of disks in this case, whose positions were defined by starting with 1500 points equally spaced on the vertical axis, and determining their (horizontal) location on the spheroid surface using the equation for the spheroid itself. This discretization resulted in 3000 total disks with a minimum width of $1.8\times 10^{-7}$ and a maximum width of $7.3\times 10^{-3}$ (near the top of the spheroid). The values for the arc lengths $\mbox {d}s$ in (4.7) were tabulated beforehand using incomplete elliptic integrals of the second kind. In short, the semi-analytical implementation herein is most efficient, and straightforward, when desingularization accuracy is not critical near the high-curvature solid-particle tips.

5.2. Single capsule

Capsules are presented as a prototype for fibrous material. First, consider a drop approaching a single capsule ($\mathcal {L}=6$). Patel et al. (Reference Patel, Shaqfeh, Butler, Cristini, Bławzdziewicz and Loewenberg2003) introduced this configuration experimentally, defining the capillary number the same way as herein, and provided boundary-integral results for a drop approaching an infinitely long cylinder. Patel et al. avoided the Stokes paradox (due to an infinite cylinder extent in the third dimension) by using the flow field around the cylinder, which would exist in the absence of the drop, based on the Brinkman equation for an isotropic porous medium, to approximately model the far-field effect of the surrounding fibres in the dilute limit. The boundary-integral problem for the drop was then solved assuming that the drop is immersed into this field, without actual interaction with the cylinder. In the present case, we consider a capsule of finite length only, but solve the problem with full hydrodynamical interaction between the drop and the capsule. The effect of particle length is discussed in more detail in the context of the two-capsule constriction.

Patel et al. (Reference Patel, Shaqfeh, Butler, Cristini, Bławzdziewicz and Loewenberg2003) identified two modes of drop breakup: a grazing mechanism in which drop breakup occurs downstream from the fibre, and hairpin formation when the drop passes around either side, resulting in two bulbous ends, or lobes, connected by a thin filament. The mode of breakup primarily depends on the initial horizontal offset ${\rm \Delta} y_o$ of the drop centre from the centreline of the capsule. In all cases for the present single-capsule study, system parameters are $Ca=3.0, \lambda =1.0$ and $\tilde {a}=1.5$, which are typical values used by Patel et al. (Reference Patel, Shaqfeh, Butler, Cristini, Bławzdziewicz and Loewenberg2003), and the drop is initially placed $10L$ upstream from the capsule. Hairpin formation was observed for ${\rm \Delta} y_o\leq 0.6$, the grazing breakup mechanism was observed for ${\rm \Delta} y_o = 0.7$ and the drop remains stable (separation without breakup) for ${\rm \Delta} y_o\geq 0.75$. The drop velocity $U$ corresponding to various mechanisms is shown vs. time in figure 7(a). Despite the apparent similarity between ${\rm \Delta} y_o = 0$ and ${\rm \Delta} y_o = 0.1$, overall drop behaviour is highly sensitive to the initial offset. For example, for an offset of ${\rm \Delta} y_o = 0.1$, 59 % of the drop (by volume) passes around one side of the capsule. Typical configurations during hairpin formation are shown in figure 7(b). As discussed further below, a large gap (${\approx }0.2$) persists between the filament and the solid particle, even after the majority of drop has passed downstream. It is most noticeable for ${\rm \Delta} y_o = 0$ in figure 7(c), that the gap does not go to zero while the drop is able to stay wrapped around the capsule for a long time. The drop diameter, defined as the maximum distance between any points on the drop and normalized by drop radius, is shown in figure 7(d), highlighting the system's sensitivity to initial drop offset.

Figure 7.

Droplet flow around a single capsule ($Ca=3.0, \lambda =1.0, \tilde {a}=1.5, \mathcal {L}=6.0$). (a) Volume-averaged drop velocity $vs.$ time as a function of drop offset. The drop remains stable for ${\rm \Delta} y_o\geq 0.75$. (b) Typical configurations during the hairpin breakup mechanism. (c) Drop–particle gap as a function of time and drop offset. (d) Drop diameter is sensitive to offset. Length is the maximum distance between any two points on the drop normalized by the undeformed drop radius (curves for ${\rm \Delta} y_o = 0.6$ and $0.7$ are very similar to ${\rm \Delta} y_o = 0.5$).

The lubrication layer between the filament and the fibre influences the drop–particle gap, and this layer is affected by the shape of the filament. This effect is demonstrated using a drop offset of ${\rm \Delta} y_o = 0.3$ as shown in figure 8(a), which results in a drop shape resembling those observed experimentally by Patel et al. (Reference Patel, Shaqfeh, Butler, Cristini, Bławzdziewicz and Loewenberg2003) (see figure 4 therein). As the filament initially nears the particle face, it spreads axially and an elongated dimple forms along the transverse direction. During this phase, the gap plateaus with respect to time, before decreasing abruptly as the filament thins. This decrease in separation corresponds to the flattening of the dimple into a convex filament. Cross sections of the filament directly above the capsule centre are shown as insets in figure 8(b). The dimple formation shown at $t=30$ is dependent on the capsule length; for example, using a capsule of length $\mathcal {L}=10$, the filament remains less concave at it approaches the cylindrical surface.

Figure 8.

Droplet undergoing large deformation around a capsule ($Ca=3.0$, $\lambda =1.0$, $\tilde {a}=1.5$, $\mathcal {L}=6.0$). (a) Droplet approaching breakup via the hairpin mechanism. (b) Rate of lubrication layer thinning correlates with the cross-sectional shape of the filament (cross-sectional plane indicated by grey line in (a)).

Our simulations have also revealed the importance of full hydrodynamical interactions for accurate prediction of the drop trajectories around the capsule. When $\lambda =1$, neglecting the capillary contribution to (3.4) when $\boldsymbol {y}$ is on the solid surface, simply makes (3.2) and (3.1) the boundary-integral formulation for the solitary drop immersed in the flow field, which would exist around the solitary capsule. Such an approximation, though, is no longer accurate when the drop comes near a solid surface. For ${\rm \Delta} y_o=0.3$, large deviations in the drop shape, position, and velocity between the rigorous and approximate solutions were observed even at $t=19$, when the drop is just beginning to wrap around the capsule and the gap is still $O(1)$.

5.3. Two-capsule constriction

Consider a drop approaching two parallel capsules separated by a gap $\epsilon =0.5$, where the non-deformed drop diameter is twice the gap and half the solid-particle diameter. In contrast to drop interactions with a single capsule in unconfined flow, the drop tends to remain compact throughout the squeezing process. An important binary property of this system is whether or not the drop will pass through the constriction; as shown in figure 9(a), the drop is trapped below a critical capillary number. There is excellent agreement between the multimesh ($\hat {N}_{\triangle }=20.2\ \text {K}\times 4\times 16$) and semi-analytical () methods, while the semi-analytical calculations were four to eight times faster. As observed for drops squeezing between spheres, the minimum drop–particle gap is almost insensitive to $Ca$ (figure 9b), although this separation remains slightly larger compared to drop squeezing through a three-sphere constriction (Zinchenko & Davis Reference Zinchenko and Davis2006). The tendency of the trailing end of the drop to approach the solid surface more closely as it is exiting, as compared to the majority of the squeezing process (gap decrease ${\approx }30\,\%$), is more pronounced for capsules than, e.g. drop squeezing between three spheres. Presumably, the differing nature of the lubrication layers is responsible for this effect; as discussed in § 5.5, the effect is even more pronounced for drops squeezing between flat plates. A snapshot of a supercritical drop when it is approximately centred in the constriction is shown in figure 9(c). The drop profile at the centre of the constriction is outlined by the dashed line in the inset. This dimple in the drop surface is less pronounced during the early and late stages of squeezing. A steady-state trapped drop is also shown from axial and transverse perspectives (figure 9d). The transverse drop profile is approximately circular, as is the case for supercritical drops at a similar position.

Figure 9.

Droplet squeezing between two capsules ($\lambda =4.0$, $\tilde {a}=0.5$, $\epsilon =0.5$, $\mathcal {L}=6.0$). (a) Droplet becomes trapped below a critical capillary number. An agreement is observed between semi-analytical and multimesh methods ($\hat {N}_{\triangle }=20.2\ \text {K}\times 4\times 16$). (b) The minimum drop–particle gap is very small for tight-squeezing simulations, and not sensitive to $Ca$. (c) Snapshot as drop passes through a two-capsule constriction above $Ca_{crit}$. Dashed line in inset outlines drop profile at constriction centre. Two orthogonal views shown. (d) Steady-state drop trapped between two capsules.

Drop squeezing is significantly affected by capsule length, even for long capsules with ends far from the drop. For example, repeating the two-sphere simulation above ($\epsilon =0.5$, $\tilde {a}=0.9$, $Ca=0.63$, $\lambda =4.0$), but replacing the spheres with capsules of length $\mathcal {L}=4$, results in vastly different behaviour. As shown in figure 10(a), the minimum drop velocity $U_{min}$ is $5\,\%$ of far-field velocity for the two-sphere constriction, whereas the drop becomes trapped in the two-capsule constriction. A parametric study of the effect of cylinder length on droplet trapping was completed at $Ca=1.5$, $\lambda =4.0$, $\tilde {a}=0.5$ and $\epsilon =0.5$. One might expect drop velocities to converge relatively quickly after the capsules reach a certain length, due the diminishing end effects. However, a notable difference in drop velocities and squeezing times was observed between $\mathcal {L}$ values of 4, 6, 8 and 10, as seen in figure 10(b). This disparity was not significantly alleviated by rescaling the capillary number individually for each system, based on the flow velocity at the centre of the constriction (in an attempt to minimize the effect of the far-field flow $\boldsymbol {u}_{\infty }$). The lack of convergence is attributed to variations in the undisturbed (single-phase) flow and can be directly related to the size-effect problems that plague both experimental and numerical studies of the Stokes paradox for flow past long cylinders. For example, Khalili & Liu (Reference Khalili and Liu2017) demonstrate that, for 3-D cylinders, large aspect ratios are required for convergence of results as Reynolds number goes to zero. Furthermore, velocity fields do not converge with decreasing Reynolds number (although unintuitive, this lack of convergence is less surprising given there is no solution for unbounded Stokes flow around an infinite cylinder). Indeed, the aspect ratios of cylinders tractable for this study are orders of magnitude smaller than Khalili & Liu's recommendations for converging velocity and pressure fields. A related work is that of Smith (Reference Smith1990) on the Jeffery paradox (an unintuitive result that involves two cylinders in two dimensions, at least one of which is rotating). Smith came to a similar conclusion that, when attempting to approach this 2-D result using finite 3-D cylinders, cylinder aspect ratios must be many orders of magnitude before end effects become negligible. However, it is not the objective of this work to approximate drop squeezing between infinite cylinders, but rather to model a local environment within fibrous material. Furthermore, a model of porous media that consists of particles in a periodic simulation box, where the flow is either pressure-driven or constant flow rate, is no longer subject to the Stokes paradox, even for effectively infinite cylinders that connect across periodic boundaries. Such a simulation would allow for a more direct comparison to models of flow through porous media, such as the Brinkman equations. The default $\mathcal {L}=6$ herein was chosen as a compromise between capsule length and simulation speed.

Figure 10.

Effect of capsule aspect ratio on drop squeezing. (a) Cylindrical particles, even at small aspect ratios ($\mathcal {L}=4$), significantly enhance drop trapping ($Ca=0.63$, $\lambda =4.0$, $\tilde {a}=0.9$, $\epsilon =0.5$). (b) Drop trajectories do not converge with respect to capsule length, $\mathcal {L}=4$ to $10$ ($Ca=1.5$, $\lambda =4.0$, $\tilde {a}=0.5$).

Drop behaviour is also dependent on viscosity ratio $\lambda$ and non-deformed droplet radius $\tilde {a}$. A parametric study of $\lambda$ at $Ca=1.5$ is shown in figure 11(a). At low $\lambda$ and high $Ca$, the drop undergoes greater variation in velocity during the squeezing process, and the overall trajectory is more symmetric with respect to time (including the appearance of two global minima for $U$). Although squeezing times increase with $\lambda$, there is no indication of switching from droplet pass through to trapping within reasonable $\lambda$ values. As observed for the three-sphere constriction (Zinchenko & Davis Reference Zinchenko and Davis2006), the drop–particle gap is sensitive to $\lambda$, decreasing from about 0.021 for $\lambda =10$, to 0.015 for $\lambda =0.25$. Drop velocity as a function of $Ca$ at low viscosity ratio ($\lambda = 0.25$) is shown in figure 11(b) (for comparison at $\lambda =4.0$, see figure 9a). Parametric studies for $Ca$ also were completed at $\lambda =1.0$ and 10; in all cases, the viscosity ratio was not observed to alter pass through $vs.$ trapping, as has been noted for drop squeezing between spheres (Zinchenko & Davis Reference Zinchenko and Davis2006). Drop size, in contrast, notably affects $Ca_{crit}$. Increasing the drop radius by $50\,\%$ results in $Ca=1.5$ being subcritical (figure 11c). For $\tilde {a}/\epsilon =2.0$, where the drop diameter is four times the gap diameter, the drop must undergo large deformation to pass through. A steady-state trapped drop is shown in figure 11(d). For this drop size, our simulations indicated pass-through at $Ca=3$. At $Ca=4$, the drop flows around the capsules rather than through the constriction (similar to the head on approach of a drop toward a single capsule, discussed above).

Figure 11.

Drop squeezing behaviour with respect to various system parameters ($\mathcal {L}=6$). (a) Drop velocity is more symmetric with respect to time at lower viscosity ratio $\lambda$ and high $Ca$ ($Ca=1.5$, $\tilde {a}/\epsilon =1.0$). (b) Effect of $Ca$ at low $\lambda =0.25$; $Ca_{crit}$ is not sensitive to $\lambda$, e.g. compared to $\lambda =4.0$ (figure 9a). (c) Tendency for drop trapping increases with droplet radius ($Ca=1.5$, $\lambda =4.0$). (d) Typical trapped state for a droplet radius of twice the gap diameter ($Ca=2.0$, $U=1{\times }10^{-4}$).

5.4. Half-capsule constriction

To demonstrate the ability of the multimesh desingularization method to handle more complex geometries and new related physical phenomena, consider a capsule that is cut in half lengthwise with rounded edges (as described in § 2). Now, the upstream single-phase flow field depends on the particle orientation, as the particles are no longer symmetric with respect to the constriction's basal plane. Therefore, droplet shape and velocity also depend on the particle orientation (or equivalently, the direction from which the drop approaches the constriction). This dependence is shown in figure 12(a), where ‘upper half-capsule’ refers to a capsule whose volume in the upper half-space has been retained. For the supercritical case $Ca=1.5$, a slight but notable difference in $U_{min}$ and squeezing times is observed (${\approx }33\,\%$ decrease in $U_{min}$). For the subcritical $Ca=0.9$, particle orientation manifests as a disparity in droplet deceleration as it approaches the constriction (figure 12a). The drops have very different shapes at analogous stages of the squeezing process. For example, the drops are compared in figure 12(c) for each orientation at their respective $U_{min}$. These configurations correspond to different lubrication behaviour. At their respective $U_{min}$, the drop–particle gap for the lower half-capsule is ${\approx }40\,\%$ lower than the upper half-capsule (similar to the subcritical $Ca=0.9$ case), and the trajectory-minimum gap disparity is ${\approx }20\,\%$ (figure 12b). While there still occurs the characteristic dimple formation near the high-curvature edges of the lower half-capsule, the surface area of this region is smaller than for the upper half-capsule case. This finding affects the interfacial velocity field, as shown in figure 12(d). As compared to a drop in near contact with a sphere, the surface-divergent flow near the centre of the near-contact region is elongated in the axial direction.

Figure 12.

Drop squeezing between two half-capsules ($\lambda =4.0$, $\tilde {a}=0.6$, $\epsilon =0.5$, $\mathcal {L}=6$, $\hat {N}_{\triangle }=32.8\ \text {K}\times 4\times 16$). (a) Drop velocity is sensitive to particle orientation (or equivalently, direction of droplet approach). The disparity between minimum velocity $U_{min}$ and squeezing times ($Ca=1.5$) indicates weak flow rectifying behaviour. (b) Minimum drop–particle gap as a function of time and particle orientation. (c) Snapshots of each drop at their respective global minimum velocities with respect to particle orientation. (d) The interfacial velocity fields of the drops pictured in (c), from a transverse view.

5.5. Two-plate constriction

Multiphase flow past fixed surfaces with high curvature presents particular numerical challenges for boundary-integral methods. In this case, the undeformed drop has radius $\tilde {a}=0.75$ and the rounded edge of the plate has a radius of 0.2, representing a curvature ratio of 3.75 (where the characteristic length $L$ is the gap $\epsilon =1$ between plates). Note that our algorithm requires smooth versus sharp edges, and that a sharp edge would lead to physical contact due to reduced lubrication. We did not explore the effect of the rounding radius. As described is § 2, the major plate dimensions are $4{\times }8$, and the constriction is formed by placing two plates on a plane with the gap $\epsilon$ defined as the minimum distance between their longest sides. This configuration has a lower $Ca_{crit}$ than the systems above, despite the lower drop-radius-to-gap ratio. As shown in figure 13(a), $Ca=0.9$ is supercritical and has a similar $U_{min}$ and squeezing time as $Ca=1.5$. The drop also passed through the constriction with capillary numbers of 0.7 and 0.5, but reached subcritical behaviour with $Ca=0.3$. The features described in § 4.4 are collectively able to demonstrate drop trapping in this constriction to sufficiently low velocities, though the simulation was not stable indefinitely; eventually the drop mesh started to overlap with the solid particle. The ability of the multimesh method to model this challenging system efficiently, without code parallelization, is promising, especially considering the extreme approach of the drop interface to the solid particle surfaces. As shown in figure 13(b), the trailing end of the drop approaches the particle to within ${\approx }0.006$ when exiting the constriction, similar to the steady-state separation of the trapped drop. This contrasts with the two-capsule system above, where the minimum gap is ${\approx }0.02$, i.e. a decrease of ${\approx }70\,\%$. Snapshots during the squeezing process are shown in figure 13(c), including the close approach at points of high drop curvature while exiting the constriction. Close inspection of the drop at $t=150$ also shows a profile with sharp angles in the near-contact region. This profile is not a mesh effect (the coordinated adaptive meshing scheme ensures high drop mesh resolution in this region at all times), but rather is due to the elongated dimple in the drop interface, which approaches the solid particle most closely along its perimeter. The evolution of these drops, suspended in far-field flow and squeezing between rectangular plates, is visually similar from this axial perspective to those of buoyant drops squeezing through circular orifices, see, e.g. Bordoloi & Longmire (Reference Bordoloi and Longmire2014) (figure 13a therein) and Ratcliffe, Zinchenko & Davis (Reference Ratcliffe, Zinchenko and Davis2012).

Figure 13.

Drop squeezing between two rounded plates ($\lambda =4.0$, $\tilde {a}=0.75$, $\epsilon =1.0$, $\hat {N}_{\triangle }=45.7\ \text {K}\times 4\times 16$). (a) Supercritical drop velocity vs. time for squeezing between flat plates is notably different than those for more rounded shapes. (b) Drops closely approach rounded flat plates, especially while the drop is exiting the constriction, or for subcritical $Ca$. (c) Evolution of drop shape during the tight-squeezing process ($Ca=0.9$).

Finally, consider a drop approaching two plates that are angled 45$^{\circ }$ with respect to the far-field flow. The objective of this configuration was to emphasize the slight drop rectifying flow displayed by the half-capsule system, perhaps to the extent of demonstrating drop pass through when approaching from one direction, and trapping or breakup from the other. However, for $Ca = 0.9\text {--}1.5$, stable supercritical behaviour is observed with either plate orientation, and furthermore $U_{min}$ is not sensitive to orientation. At $Ca=2.0$, the global velocity minimum is lower for upward-angled plates (each plate rotated such that its edge near the constriction centre is moved upstream). A similar disparity is observed at $Ca=3.0$, as shown in figure 14(a). However, in this case, while the drop passes stably through the downward-angled plates (figure 14b), it will break when flowing through the upward plates. The mechanism for this breakup is shown in figure 14(c), suggesting the majority of drop volume passes through the constriction, while two significantly smaller drops will flow around the outer faces of each plate.

Figure 14.

Drop squeezing between two plates angled at 45$^{\circ }$ with respect to the far-field flow ($Ca=3.0$, $\lambda =4.0$, $\tilde {a}=0.75$, $\epsilon =1.0$, $\hat {N}_{\triangle }=45.7\ \text {K}\times 4\times 16$). (a) Drop flows stably between downward-angled plates. The $X$ mark indicates simulation exit before imminent drop breakup in the case of plates angled upward. (b) Snapshot of supercritical drop when approximately centred in the downward-angled plate constriction. (c) Mechanism of drop breakup in the case of a drop squeezing between upward-angled plates ($Ca=3.0$).