2.1. Implementation of periodic boundaries

Simple shear. The reproducible lattice of points $\{m^{i}\boldsymbol{e}_{i}\}$ , where $m^{i}$ are arbitrary integers, has long been known for shear flow (2.1a ), and is extensively used in molecular dynamics and suspension flow simulations. It has also received applications in emulsion rheology studies (Loewenberg & Hinch Reference Loewenberg and Hinch1996, paper I). At $t=0$ , the basic cell $V$ is a cube $[0,L)^{3}$ , with the side $L$ , which we use as the length scale. For $t>0$ the basis vectors $\boldsymbol{e}_{1}=(1,0,0)$ and $\boldsymbol{e}_{3}=(0,0,1)$ are unchanged. Once the strain $\dot{{\it\gamma}}t$ crosses a half-integer value, the lattice is restructured by taking $\boldsymbol{e}_{2}-\boldsymbol{e}_{1}$ as the new $\boldsymbol{e}_{2}$ -vector, to make, in general, $\boldsymbol{e}_{2}=({\it\gamma},1,0)$ , where ${\it\gamma}=\dot{{\it\gamma}}t-[\dot{{\it\gamma}}t+1/2]\in [-1/2,1/2)$ and the brackets denote the integer part. Figure 1(a) (two-dimensional slice) shows the simple mapping to translate the drop centres from the old cell $V$ to the new cell to continue the simulation. Such restructuring ensures the minimal possible skewage at all $t>0$ .

Figure 1. Mapping between the old and new (shaded) periodic cells at the restructuring moments for (a) simple shear and (b) planar extension. Region $\text{I}$ is translated to $\text{I}^{\prime }$ , region $\text{II}$ to $\text{II}^{\prime }$ , etc. The third dimension $x_{3}$ is not shown.

Planar extension. Far less obvious is the existence of reproducible lattices in PE found by Kraynik & Reinelt (Reference Kraynik and Reinelt1992). Following their results for a ‘critical lattice’ (with minimal skewage), we choose a square with side $L$ inclined at angle ${\it\theta}_{o}=\tan ^{-1}[(\sqrt{5}-1)/2]\approx 31.7^{\circ }$ to the $x_{1}$ -axis as the two-dimensional slice of the initial cubic periodic cell $V$ (figure 2 a). Again, $L$ is used as the length scale. The initial lattice vectors then become $\boldsymbol{e}_{1}^{o}=(\cos {\it\theta}_{o},\sin {\it\theta}_{o},~0)$ and $\boldsymbol{e}_{2}^{o}=(-\text{sin}\,{\it\theta}_{o},\cos {\it\theta}_{o},~0)$ . Figure 2 shows how the lattice reproducibility is achieved in the Hencky strain interval of $\dot{{\it\Gamma}}t={\it\varepsilon}_{p}$ , with ${\it\varepsilon}_{p}=\ln [(3+\sqrt{5})/2]\approx 0.9624$ . By $\dot{{\it\Gamma}}t={\it\varepsilon}_{p}/2$ , the lattice dynamics $\text{d}\boldsymbol{e}_{i}/\text{d}t=\text{diag}(\dot{{\it\Gamma}},-\dot{{\it\Gamma}},0)\boldsymbol{e}_{i}$ brings the basic periodic cell $V$ to a shape shown as the dashed-line parallelogram in figure 2(c). At this time moment, switching to $\boldsymbol{e}_{1}+\boldsymbol{e}_{2}$ and $\boldsymbol{e}_{1}+2\boldsymbol{e}_{2}$ as the new lattice vectors and integrating the lattice dynamics further to $\dot{{\it\Gamma}}t={\it\varepsilon}_{p}$ reproduces the initial cell shown in figure 2(a), due to the choice of ${\it\theta}_{o}$ and ${\it\varepsilon}_{p}$ . In simulations, such a lattice restructuring must be carried out at all half-integer values of $\dot{{\it\Gamma}}t/{\it\varepsilon}_{p}=n+1/2$ , but there is no need to match $(n+1/2)$ precisely. The vectors $\boldsymbol{e}_{1}$ and $\boldsymbol{e}_{2}$ are replaced with $\boldsymbol{e}_{1}+\boldsymbol{e}_{2}$ and $\boldsymbol{e}_{1}+2\boldsymbol{e}_{2}$ once $\dot{{\it\Gamma}}t$ crosses a half-integer value, to make, in general,

(2.3a,b ) $$\begin{eqnarray}\boldsymbol{e}_{1}=(\cos {\it\theta}_{o}\text{e}^{{\it\Gamma}},\sin {\it\theta}_{o}\text{e}^{-{\it\Gamma}},0),\quad \boldsymbol{e}_{2}=(-\text{sin}\,{\it\theta}_{o}\text{e}^{{\it\Gamma}},\cos {\it\theta}_{o}\text{e}^{-{\it\Gamma}},0),\end{eqnarray}$$

where

(2.4) $$\begin{eqnarray}{\it\Gamma}=\dot{{\it\Gamma}}t-{\it\varepsilon}_{p}\left[\frac{1}{2}+\frac{\dot{{\it\Gamma}}t}{{\it\varepsilon}_{p}}\right]\in \left[-\frac{{\it\varepsilon}_{p}}{2},\frac{{\it\varepsilon}_{p}}{2}\right).\end{eqnarray}$$

Figure 2. Dynamics of the computational cell in PE. At the Hencky strain of $\dot{{\it\Gamma}}t={\it\varepsilon}_{p}/2$ , the cell is restructured (the shaded area is the new cell) to continue the simulation. The third dimension is not shown; (a) $\dot{{\it\Gamma}}t=0$ , (b) $\dot{{\it\Gamma}}t={\it\varepsilon}_{p}/4$ , (c) $\dot{{\it\Gamma}}t={\it\varepsilon}_{p}/2$ , (d) $\dot{{\it\Gamma}}t=3{\it\varepsilon}_{p}/4$ .

Prior to a jump in $\boldsymbol{e}_{1}$ and $\boldsymbol{e}_{2}$ , mapping between the old and new periodic cells (figure 1 b) is used to recalculate the new drop positions after the cell restructuring. Namely, drops with centroids in region I are translated by $\boldsymbol{e}_{1}+\boldsymbol{e}_{2}$ to region $\text{I}^{\prime }$ of the new cell. Region II is translated by $\boldsymbol{e}_{2}$ to $\text{II}^{\prime }$ , region III – by $\boldsymbol{e}_{1}+2\boldsymbol{e}_{2}$ to $\text{III}^{\prime }$ .

Strong mixed flows. Kraynik & Reinelt (Reference Kraynik and Reinelt1992) noted that their reproducible lattice dynamics for PE can be generalized to other macroscopic flows with a constant velocity gradient matrix $\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{\infty }$ , if a diagonalizing transformation

(2.5) $$\begin{eqnarray}\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{\infty }=\unicode[STIX]{x1D64E}~\text{diag}(\dot{{\it\Gamma}},-\dot{{\it\Gamma}},0)\unicode[STIX]{x1D64E}^{-1}\end{eqnarray}$$

exists between $\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{\infty }$ and the velocity gradient matrix for PE, with a real-valued matrix $\unicode[STIX]{x1D64E}$ . Indeed, the required lattice vector dynamics $\text{d}\boldsymbol{e}_{i}/\text{d}t=\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{\infty }\boldsymbol{\cdot }\boldsymbol{x}$ is reduced to the PE flow dynamics $\text{d}\widetilde{\boldsymbol{e}}_{i}/\text{d}t=\text{diag}(\dot{{\it\Gamma}},-\dot{{\it\Gamma}},0)\widetilde{\boldsymbol{e}}_{i}$ for auxiliary vectors $\widetilde{\boldsymbol{e}}_{i}=\unicode[STIX]{x1D64E}^{-1}\boldsymbol{e}_{i}$ . The restructuring transformation $(\,\widetilde{\boldsymbol{e}}_{1},~\widetilde{\boldsymbol{e}}_{2})\rightarrow (\,\widetilde{\boldsymbol{e}}_{1}+\widetilde{\boldsymbol{e}}_{2},~\widetilde{\boldsymbol{e}}_{1}+2\widetilde{\boldsymbol{e}}_{2})$ to the new basis vectors in the PE flow space at half-integer values of $\dot{{\it\Gamma}}t/{\it\varepsilon}_{p}$ begets the same transformation of the lattice vectors $\boldsymbol{e}_{i}=\unicode[STIX]{x1D64E}\widetilde{\boldsymbol{e}}_{i}$ .

For a strong mixed flow, the transformation (2.5) was found by Hunt et al. (Reference Hunt, Bernardi and Todd2010) and used in sample molecular dynamics simulations. Their definition of mixed flows, though, differs from our (2.1c ) in the choice of the coordinate axes (rotated by $45^{\circ }$ ). For our definition (2.1c ) and ${\it\chi}>0$ (characterizing ‘strong’ flows), the relation between $\dot{{\it\Gamma}}$ in (2.5) and $\dot{{\it\gamma}}_{m}$ is $\dot{{\it\Gamma}}=\dot{{\it\gamma}}_{m}\sqrt{{\it\chi}}$ ; the matrix $\unicode[STIX]{x1D64E}$ and its inverse take the forms

(2.6a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D64E}=\left(\begin{array}{@{}ccc@{}}(4{\it\chi})^{-1/4} & -(4{\it\chi})^{-1/4} & 0\\ ({\it\chi}/4)^{1/4} & ({\it\chi}/4)^{1/4} & 0\\ 0 & 0 & 1\end{array}\right),\quad \unicode[STIX]{x1D64E}^{-1}=\left(\begin{array}{@{}ccc@{}}({\it\chi}/4)^{1/4} & (4{\it\chi})^{-1/4} & 0\\ -({\it\chi}/4)^{1/4} & (4{\it\chi})^{-1/4} & 0\\ 0 & 0 & 1\end{array}\right).\end{eqnarray}$$

An arbitrary factor in $\unicode[STIX]{x1D64E}$ is chosen to make $\det (\unicode[STIX]{x1D64E})=1$ , which preserves the cell volume in transformations between $\widetilde{\boldsymbol{e}}_{i}$ and $\boldsymbol{e}_{i}$ . The side of the initial cubic box span on $\widetilde{\boldsymbol{e}}_{1}^{\,o}$ , $\widetilde{\boldsymbol{e}}_{2}^{\,o}$ and $\widetilde{\boldsymbol{e}}_{3}^{\,o}$ (see figure 2 a) is used as the length scale in our mixed-flow algorithm for ${\it\chi}>0$ . It can be derived from (2.3)–(2.6) that $\boldsymbol{e}_{1}^{2}+\boldsymbol{e}_{2}^{2}=(1/2)({\it\chi}^{1/2}+{\it\chi}^{-1/2})$ $(\,\widetilde{\boldsymbol{e}}_{1}^{2}+\widetilde{\boldsymbol{e}}_{2}^{2})$ , indicating only a modest additional skewage of the periodic cell in mixed-flow simulations compared with PE simulations (except for very small ${\it\chi}$ with a stronger effect). For ${\it\chi}=0.25$ , the dynamics of the periodic cell $V$ in the interval $0\leqslant \dot{{\it\gamma}}_{m}t<2{\it\varepsilon}_{p}$ is shown in figure 3; at $\dot{{\it\gamma}}_{m}t=2{\it\varepsilon}_{p}$ , the cell $V$ returns to its initial shape. It can be noted that the formal mapping (2.6) between PE and mixed flow is only used to implement the periodic boundary conditions for the latter; of course, mixed-flow rheological simulations are irreducible to PE simulations.

Figure 3. The geometry of the computational cell in a strong mixed flow ( ${\it\chi}=0.25$ ) at $\dot{{\it\gamma}}_{m}t/{\it\varepsilon}_{p}=0$ (a), 1∕2 (b), 2∕3 (c), 1 (d), 4∕3 (e) and 5∕3 (f). In (d), the cell is restructured (the shaded area is the new cell) to continue the simulation.

Weak mixed flows. For ${\it\chi}<0$ (when the transformation (2.5) does not exist), implementation of periodic boundaries is much simpler and does not require the computational cell restructuring. As for simple shear, we choose the cube $[0,L)^{3}$ as the initial periodic cell and use $L$ as the length scale. In contrast to simple shear and PE, though, the lattice vector dynamics $\text{d}\boldsymbol{e}_{i}/\text{d}t=\boldsymbol{{\rm\nabla}}\boldsymbol{u}_{\infty }$ . $\boldsymbol{e}_{i}$ does not lead to indefinite stretching/contraction (which would require cell restructuring), but gives the periodic solution instead:

(2.7a,b ) $$\begin{eqnarray}\boldsymbol{e}_{1}=(\cos {\it\omega}t,-|{\it\chi}|^{1/2}\sin {\it\omega}t,0),\quad \boldsymbol{e}_{2}=(|{\it\chi}|^{-1/2}\sin {\it\omega}t,\cos {\it\omega}t,0),\end{eqnarray}$$

with ${\it\omega}=|{\it\chi}|^{1/2}\dot{{\it\gamma}}_{m}$ . Figure 4 illustrates the degree of compactness of the periodic cell for ${\it\chi}=-0.25$ ; the selected shapes at ${\it\omega}t=0$ , ${\rm\pi}/4$ , ${\rm\pi}/2$ and $3{\rm\pi}/4$ also give the periodic cells at ${\it\omega}t={\rm\pi}$ , $5{\rm\pi}/4$ , $3{\rm\pi}/2$ and $7{\rm\pi}/4$ respectively, through central symmetry about the origin $\boldsymbol{O}$ . Difficulties with bad aspect ratios would appear for much smaller $|{\it\chi}|$ (not considered herein); for the shear-flow limit ${\it\chi}=0$ , an alternative standard implementation of periodic boundaries (see above) is used instead.

Figure 4. The geometry of the computational cell in a weak mixed flow ( ${\it\chi}=-0.25$ ). The directions of the $x_{1}$ and $x_{2}$ coordinate axes are the same as in figure 3; (a)  ${\it\omega}t=0$ , (b) ${\rm\pi}/4$ , (c) ${\rm\pi}/2$ , (d)  $3{\rm\pi}4$ .

2.2. Boundary-integral formulation

For all three flows (2.1), the time scale is defined as $(2I_{2})^{-1/2}$ , the velocity scale is $(2I_{2})^{1/2}L$ . The non-dimensional system of BI equations for the interfacial fluid velocity $\boldsymbol{u}(\boldsymbol{y})$ on drop surfaces $S_{{\it\alpha}}~({\it\alpha}=1,\ldots ,N)$ , to be solved iteratively at each time step, is used in the form of paper I:

(2.8) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(\boldsymbol{y}) & = & \displaystyle \frac{({\it\lambda}-1)}{({\it\lambda}+1)}\left[2\mathop{\sum }_{{\it\beta}=1}^{N}\int _{S_{{\it\beta}}}\boldsymbol{Q}(\boldsymbol{x})\boldsymbol{\cdot }\widetilde{{\bf\tau}}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x})\,\text{d}S_{x}+\langle \boldsymbol{u}\rangle _{{\it\alpha}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,2\mathop{\sum }_{{\it\beta}=1}^{N}\int _{S_{{\it\beta}}}(\boldsymbol{Q}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{x}-\boldsymbol{x}_{{\it\beta}}^{c})\,\text{d}S_{x}\right]+\boldsymbol{F}(\boldsymbol{y}),\quad \boldsymbol{y}\in S_{{\it\alpha}}.\end{eqnarray}$$

Here, the inhomogeneous term is

(2.9) $$\begin{eqnarray}\boldsymbol{F}(\boldsymbol{y})=\frac{2}{{\it\lambda}+1}\left[\boldsymbol{u}_{\infty }(\boldsymbol{y})+\frac{2a}{\mathit{Ca}}\mathop{\sum }_{{\it\beta}=1}^{N}\int _{S_{{\it\beta}}}(k(\boldsymbol{x})-\langle k\rangle _{{\it\beta}})\boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x}\right],\end{eqnarray}$$

where $a=a_{o}/L$ is the non-dimensional non-deformed drop radius, $k(\boldsymbol{x})$ is the local mean of the principal surface curvatures, $\langle k_{{\it\beta}}\rangle$ is the average of $k$ over $S_{{\it\beta}}$ , $\boldsymbol{n}(\boldsymbol{x})$ is the outward unit normal at $\boldsymbol{x}\in S_{{\it\beta}}$ and $\unicode[STIX]{x1D642}(\boldsymbol{x})$ is the triply periodic symmetric second-rank Green tensor. This tensor satisfies

(2.10) $$\begin{eqnarray}{\rm\nabla}^{2}\unicode[STIX]{x1D642}(\boldsymbol{x})-\boldsymbol{{\rm\nabla}}\boldsymbol{q}(\boldsymbol{x})=\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }{\bf\tau}(\boldsymbol{x})=\mathop{\sum }_{\boldsymbol{m}}{\it\delta}(\boldsymbol{x}-\boldsymbol{P}(\boldsymbol{m}))\unicode[STIX]{x1D644},\end{eqnarray}$$

where $\boldsymbol{q}(\boldsymbol{x})$ is the vector of pressures for $\unicode[STIX]{x1D642}(\boldsymbol{x})$ and the summation is over all lattice points $\boldsymbol{P}(\boldsymbol{m})=m^{1}\boldsymbol{e}_{1}+m^{2}\boldsymbol{e}_{2}+m^{3}\boldsymbol{e}_{3}$ with integer $m^{1}$ , $m^{2}$ , $m^{3}$ . The Ewald-like form for $\unicode[STIX]{x1D642}(\boldsymbol{x})$ for a lattice with arbitrary periods $\boldsymbol{e}_{1}$ , $\boldsymbol{e}_{2}$ , $\boldsymbol{e}_{3}$ follows from Hasimoto (Reference Hasimoto1959) and is given in appendix B of paper I. It is inconvenient, though, that the associated third-rank stress tensor ${\bf\tau}(\boldsymbol{x})$ is not periodic (and not symmetric), since it contains a linearly growing part from the pressure $\boldsymbol{q}(\boldsymbol{x})$ . Following Zinchenko & Davis (Reference Zinchenko and Davis2000) and paper I, we work instead with the periodic and fully symmetric tensor $\widetilde{{\bf\tau}}(\boldsymbol{r})={\bf\tau}(\boldsymbol{r})-\boldsymbol{r}\unicode[STIX]{x1D644}$ in the double-layer integrals (2.8). Additionally, velocity fluctuations $\boldsymbol{Q}(\boldsymbol{x})=\boldsymbol{u}(\boldsymbol{x})-\langle \boldsymbol{u}\rangle _{{\it\beta}}$ from the surface-averaged values $\langle \boldsymbol{u}\rangle _{{\it\beta}}$ are used on $S_{{\it\beta}}$ instead of $\boldsymbol{u}(\boldsymbol{x})$ to significantly reduce the integrand variations and thereby accelerate convergence in the most difficult multipole part (§ 3) of the algorithm. For the same reason, we work with $k(\boldsymbol{x})-\langle k\rangle _{{\it\beta}}$ instead of $k(\boldsymbol{x})$ in the single-layer integrals (2.9). It is also convenient that $\boldsymbol{Q}(\boldsymbol{x})$ (unlike $\boldsymbol{u}(\boldsymbol{x})$ ) is the same for all periodic images of $S_{{\it\beta}}$ .

Once the interfacial velocity is found, the instantaneous volume-averaged stress tensor ${\it\bf\Sigma}$ in the emulsion can be expressed via surface integrals (e.g. Batchelor Reference Batchelor1970; Kennedy et al. Reference Kennedy, Pozrikidis and Skalak1994). With ${\it\mu}_{e}(2I_{2})^{1/2}$ as the stress scale, our non-dimensional stress for all three flows (2.1) is calculated as

(2.11) $$\begin{eqnarray}\displaystyle {\it\Sigma}_{ij} & = & \displaystyle \text{I.T.}+\boldsymbol{{\rm\nabla}}_{i}(u_{\infty })_{j}+\boldsymbol{{\rm\nabla}}_{j}(u_{\infty })_{i}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{{\it\alpha}=1}^{N}\int _{S_{{\it\alpha}}}\left\{\frac{2a}{\mathit{Ca}}[k(\boldsymbol{x})-\langle k\rangle _{{\it\alpha}}]n_{i}(x-x_{{\it\alpha}}^{c})_{j}+({\it\lambda}-1)(u_{i}n_{j}+u_{j}n_{i})\!\right\}\,\text{d}S,\qquad\end{eqnarray}$$

where I.T. is an insignificant isotropic term; in deriving (2.11), the unit volume of the periodic cell has been taken into account. The tensor ${\it\bf\Sigma}$ is symmetric upon accurate surface discretization. From (2.11), all the rheological functions (2.2) are calculated.

3.1. Surface discretization and BI desingularizations

An unstructured mesh of $N_{\triangle }\gg 1$ flat triangles with $N_{\triangle }/2+2$ vertices (called nodes) is used to represent each drop surface, including the choices $N_{\triangle }=2160$ , 2880, 3840 and 6000 (in one simulation) used in the present calculations. Initial nearly uniform triangulations were prepared in a standard way from sphere tesselations with mesh refinements. As the drops move and deform, the quality of triangulations is maintained as described in § 4 (quite differently from paper I). Surface normals and curvatures are calculated by the best paraboloid-spline method of Zinchenko & Davis (Reference Zinchenko and Davis2000). Surface integrations were performed by a standard trapezoidal rule, with reassignment of triangle contributions to vertices for efficiency; the latter procedure follows Rallison (Reference Rallison1981) and was used in much of our prior work.

The double- and single-layer kernels in (2.8) and (2.9) are singular at all lattice points $\boldsymbol{x}-\boldsymbol{y}=\boldsymbol{P}(\boldsymbol{m})$ due to free-space contributions:

(3.1a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D642}(\boldsymbol{r})\simeq \unicode[STIX]{x1D642}_{o}(\boldsymbol{r}-\boldsymbol{P}(\boldsymbol{m})),\quad \widetilde{{\bf\tau}}(\boldsymbol{r})\simeq {\bf\tau}_{o}(\boldsymbol{r}-\boldsymbol{P}(\boldsymbol{m})),\end{eqnarray}$$

where

(3.2a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D642}_{o}(\boldsymbol{r})=-\frac{1}{8{\rm\pi}}\left(\frac{\unicode[STIX]{x1D644}}{r}+\frac{A\boldsymbol{A}}{r^{3}}\right),\quad {\bf\tau}_{o}(\boldsymbol{r})=\frac{3}{4{\rm\pi}}~\frac{rA\boldsymbol{A}}{r^{5}}.\end{eqnarray}$$

This behaviour necessitates efficient desingularizations of the boundary integrals prior to numerical solution. For a mesh node $\boldsymbol{y}\in S_{{\it\alpha}}$ on a drop surface different from an integration surface $S_{{\it\beta}}$ , the double-layer contribution from $S_{{\it\beta}}$ is written as (paper I)

(3.3) $$\begin{eqnarray}\displaystyle \int _{S_{{\it\beta}}}\boldsymbol{Q}(\boldsymbol{x})\boldsymbol{\cdot }\widetilde{{\bf\tau}}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x})\,\text{d}S_{x} & {\approx} & \displaystyle \mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\beta}}}\boldsymbol{Q}(\boldsymbol{x}_{j})\boldsymbol{\cdot }\widetilde{{\bf\tau}}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}_{j}){\rm\Delta}S_{j}-\mathop{\sum }_{\boldsymbol{m}}{\it\Theta}(\boldsymbol{m},\boldsymbol{y})\nonumber\\ \displaystyle & & \displaystyle \times \,\boldsymbol{Q}^{\ast }\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\beta}}}{\bf\tau}_{o}(\boldsymbol{x}_{j}+\boldsymbol{P}(\boldsymbol{m})-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}_{j}){\rm\Delta}S_{j}.\end{eqnarray}$$

Here, the summation is over all mesh nodes $\boldsymbol{x}_{j}\in S_{{\it\beta}}$ and all integer vectors $\boldsymbol{m}$ , and ${\rm\Delta}S_{j}$ is the surface area associated with node $\boldsymbol{x}_{j}$ . The barrier function ${\it\Theta}$ is designed to act only when $\boldsymbol{y}$ is close to a periodic image $S_{{\it\beta}}+\boldsymbol{P}(\boldsymbol{m})$ of $S_{{\it\beta}}$ . Namely, ${\it\Theta}$ is close to unity when the distance from $\boldsymbol{y}$ to the nearest node $\boldsymbol{x}^{\ast }+\boldsymbol{P}(\boldsymbol{m})$ on $S_{{\it\beta}}+\boldsymbol{P}(\boldsymbol{m})$ is much less than a threshold $h_{o}$ , and ${\it\Theta}$ gradually falls to zero as this distance approaches $h_{o}$ . For the precise definition of ${\it\Theta}$ , see paper I; $h_{o}$ is set to $0.3a$ in the present calculations. The choice of $\boldsymbol{Q}^{\ast }$ is crucial for successful large-strain simulations at high drop volume fractions and significant viscosity contrast ${\it\lambda}$ . The simplest near-singularity subtraction $\boldsymbol{Q}^{\ast }=\boldsymbol{Q}(\boldsymbol{x}^{\ast })$ is insufficient to overcome divergence of BI iterations even for moderate ${\it\lambda}\neq 1$ . The more robust (yet more complex) variational technique of paper I, which we follow herein, requires the subtracted quantity $\boldsymbol{Q}^{\ast }$ to minimize the Euclidean norm of the discrete free-space double-layer contribution from $S_{{\it\beta}}+\boldsymbol{P}(\boldsymbol{m})$ to $\boldsymbol{y}$ after subtraction. When ${\it\Theta}\approx 1$ , the subtracted term in (3.3) effectively alleviates the near-singular behaviour of the double-layer integrand (although it does not eliminate it completely). Of course, this subtracted term, helping calculations at finite triangulations, disappears in the limit $N_{\triangle }\rightarrow \infty$ . The optimized way of calculating the subtracted term (Zinchenko & Davis Reference Zinchenko and Davis2006, Reference Zinchenko and Davis2008) had a relatively small cost in the present simulations.

The single-layer contribution in (2.9) from a surface $S_{{\it\beta}}\not\ni \boldsymbol{y}$ is calculated as

(3.4) $$\begin{eqnarray}\displaystyle \int _{S_{{\it\beta}}}f(\boldsymbol{x})\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x})\,\text{d}S_{x} & {\approx} & \displaystyle \mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\beta}}}f(\boldsymbol{x}_{j})\unicode[STIX]{x1D642}(\boldsymbol{x}_{j}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}_{j}){\rm\Delta}S_{j}-\mathop{\sum }_{\boldsymbol{m}}{\it\Theta}(\boldsymbol{m},\boldsymbol{y})\nonumber\\ \displaystyle & & \displaystyle \times ~f^{\ast }\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\beta}}}\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}_{j}+\boldsymbol{P}(\boldsymbol{m})-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}_{j}){\rm\Delta}S_{j},\end{eqnarray}$$

with the same barrier ${\it\Theta}$ as in (3.3), and $f(\boldsymbol{x})=k(\boldsymbol{x})-\langle k_{{\it\beta}}\rangle$ for brevity. We found the simplest choice of $f^{\ast }=f(\boldsymbol{x}^{\ast })$ to be sufficient in the present calculations, although a more complex scheme was offered in paper I. The subtracted term in (3.4) completely eliminates near-singular behaviour of the integrand (thus helping calculations at finite triangulations), but, again, disappears in the limit $N_{\triangle }\rightarrow \infty$ . Due to the barrier ${\it\Theta}$ , typically, only one periodic image of $S_{{\it\beta}}$ contributes to the subtracted terms in (3.3) and (3.4), and so it was inexpensive to calculate these terms by direct summations.

For self-interactions $(\boldsymbol{y}\in S_{{\it\alpha}}=S_{{\it\beta}})$ , the contributions of $S_{{\it\beta}}$ to the BIs in (2.8) and (2.9) are computed as

(3.5) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{S_{{\it\alpha}}}\boldsymbol{Q}(\boldsymbol{x})\widetilde{{\bf\tau}}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x})\,\text{d}S_{x}\nonumber\\ \displaystyle & & \displaystyle \quad \approx \frac{1}{2}\boldsymbol{Q}(\boldsymbol{y})+\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\alpha}}}\boldsymbol{Q}(\boldsymbol{x}_{j})\boldsymbol{\cdot }[\widetilde{{\bf\tau}}(\boldsymbol{x}_{j}-\boldsymbol{y})-{\bf\tau}_{o}(\boldsymbol{x}_{j}-\boldsymbol{y})]\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}_{j}){\rm\Delta}S_{j}\nonumber\\ \displaystyle & & \displaystyle \qquad +\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\alpha}}}[\boldsymbol{Q}(\boldsymbol{x}_{j})-\boldsymbol{Q}(\boldsymbol{y})]\boldsymbol{\cdot }{\bf\tau}_{o}(\boldsymbol{x}_{j}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x}_{j}){\rm\Delta}S_{j}\end{eqnarray}$$

and

(3.6) $$\begin{eqnarray}\displaystyle \int _{S_{{\it\alpha}}}f(\boldsymbol{x})\boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x} & {\approx} & \displaystyle \mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\alpha}}}f(\boldsymbol{x}_{j})\boldsymbol{n}(\boldsymbol{x}_{j})\boldsymbol{\cdot }[\unicode[STIX]{x1D642}(\boldsymbol{x}_{j}-\boldsymbol{y})-\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}_{j}-\boldsymbol{y})]{\rm\Delta}S_{j}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{{\it\alpha}}}[f(\boldsymbol{x}_{j})-f(\boldsymbol{y})]\boldsymbol{n}(\boldsymbol{x}_{j})\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}_{j}-\boldsymbol{y}){\rm\Delta}S_{j}.\end{eqnarray}$$

For the second lines in (3.5) and (3.6) (with $\boldsymbol{x}_{j}=\boldsymbol{y}$ excluded), direct summations are always used. Unlike in (3.3), complete desingularization is achieved in (3.5).

3.2. Multipole-accelerated scheme

For the first lines in (3.3) and (3.4), it would be, of course, prohibitively expensive to use direct summations in dynamical simulations, with a cost of $O(N^{2}N_{\triangle }^{2})$ per time step and a large coefficient, due to iterations and the complex nature of $\widetilde{{\bf\tau}}$ and $\unicode[STIX]{x1D642}$ . A highly efficient multipole-accelerated scheme is used instead. Each drop, when elongated, is divided into a few (not more than three) compact blocks; see paper I for the procedure and figure 2 therein. A compact drop coincides with its only block. As a result, all mesh nodes on all drops are partitioned into compact blocks $\mathscr{B}_{1},\ldots ,\mathscr{B}_{N_{B}}~(N_{B}\geqslant N)$ . A minimal spherical shell $\mathscr{D}_{{\it\gamma}}$ with a centre $\boldsymbol{x}_{{\it\gamma}}^{o}$ and radius $d_{{\it\gamma}}^{o}$ is constructed around each block $\mathscr{B}_{{\it\gamma}}$ with sufficient accuracy.

The next step is to partition the periodic kernels $\unicode[STIX]{x1D642}(\boldsymbol{r})$ and $\widetilde{{\bf\tau}}(\boldsymbol{r})$ into the ‘near-field’ and ‘far-field’ parts. In paper I, only the $\unicode[STIX]{x1D642}_{o}(\boldsymbol{r})$ and ${\bf\tau}_{o}(\boldsymbol{r})$ parts of $\unicode[STIX]{x1D642}(\boldsymbol{r})$ and $\widetilde{{\bf\tau}}(\boldsymbol{r})$ were qualified as near-field. Such partition causes difficulties and/or performance degradation when $N$ is not too large ( ${\sim}100$ ), periodic cells are significantly more stretched than in shear flow (see figures 2–4) and some drops are strongly deformed. A more flexible and robust scheme herein is to assign free-space contributions from lattice points to near-field, if their distance from the origin is within a given threshold $R_{sub}$ :

(3.7a,b ) $$\begin{eqnarray}\unicode[STIX]{x1D642}(\boldsymbol{r})=\mathop{\sum }_{|\boldsymbol{P}(\boldsymbol{m})|\leqslant R_{sub}}\!\!\!\!\unicode[STIX]{x1D642}_{o}(\boldsymbol{r}+\boldsymbol{P}(\boldsymbol{m}))+\unicode[STIX]{x1D642}_{1}(\boldsymbol{r}),\quad \widetilde{{\bf\tau}}(\boldsymbol{r})=\mathop{\sum }_{|\boldsymbol{P}(\boldsymbol{m})|\leqslant R_{sub}}\!\!\!\!{\bf\tau}_{o}(\boldsymbol{r}+\boldsymbol{P}(\boldsymbol{m}))+{\bf\tau}_{1}(\boldsymbol{r}).\end{eqnarray}$$

The number of free-space contributions to be subtracted from $(\unicode[STIX]{x1D642},\widetilde{{\bf\tau}})$ to form the far-field $(\unicode[STIX]{x1D642}_{1},{\bf\tau}_{1})$ depends, therefore, on both $R_{sub}$ and the instantaneous lattice geometry. For $R_{sub}=1.4$ , which we found to be in the optimal range in the present simulations with $N=100{-}400$ drops, there are 7–9, 7–13 and 9–13 subtractions for simple shear, PE and mixed-flow ( ${\it\chi}=0.25$ ) lattices respectively. The partition (3.7) with $R_{sub}=1.4$ moves the singularities of $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ and ${\bf\tau}_{1}(\boldsymbol{r})$ farther away from the origin, compared with the $R_{sub}=0$ choice in paper I.

The free-space contributions of each block $\mathscr{B}_{{\it\gamma}}~({\it\gamma}=1,\ldots ,N_{B})$ to the summations in the first lines of (3.3) and (3.4) (i.e. with $\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}_{j}-\boldsymbol{y})$ and ${\bf\tau}_{o}(\boldsymbol{x}_{j}-\boldsymbol{y})$ in place of $\unicode[STIX]{x1D642}(\boldsymbol{x}_{j}-\boldsymbol{y})$ and $\widetilde{{\bf\tau}}(\boldsymbol{x}_{j}-\boldsymbol{y})$ ) are expanded in Lamb’s singular series about the block centre $\boldsymbol{x}_{{\it\gamma}}^{o}$ . These expansions are generated to a sufficient order by the alternative algorithm of Zinchenko & Davis (Reference Zinchenko and Davis2008, § 4.2 therein), which is ${\approx}2.5$ times faster for double-layer expansions than the original algorithm of Zinchenko & Davis (Reference Zinchenko and Davis2000) used in paper I for this purpose.

Now, to calculate the sums in the first lines of (3.3) and (3.4) for $\boldsymbol{y}\in \mathscr{B}_{{\it\delta}}\subset S_{{\it\alpha}}$ , we consider, for each block $\mathscr{B}_{{\it\gamma}}\not\subset S_{{\it\alpha}}$ , a periodically shifted block $\mathscr{B}_{{\it\gamma}}^{min}=\mathscr{B}_{{\it\gamma}}+\boldsymbol{P}(\boldsymbol{m}_{{\it\delta}{\it\gamma}})$ , which minimizes the centre-to-centre distance $|\boldsymbol{x}_{{\it\gamma}}^{o,min}-\boldsymbol{x}_{{\it\delta}}^{o}|$ (with $\boldsymbol{x}_{{\it\gamma}}^{o,min}=\boldsymbol{x}_{{\it\gamma}}^{o}+\boldsymbol{P}(\boldsymbol{m}_{{\it\delta}{\it\gamma}})$ ). In what follows, the index min will supplement quantities related to $\mathscr{B}_{{\it\gamma}}^{min}$ . According to (3.7) and periodicity of the kernels $\unicode[STIX]{x1D642}$ and $\widetilde{{\bf\tau}}$ , the contribution of $\mathscr{B}_{{\it\gamma}}$ to either (3.3) or (3.4) can be written as the sum of free-space contributions from $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ (with $|\boldsymbol{P}(\boldsymbol{m})|\leqslant R_{sub}$ ) plus the far-field contribution from $\mathscr{B}_{{\it\gamma}}^{min}$ using ${\bf\tau}_{1}$ and $\unicode[STIX]{x1D642}_{1}$ . The latter is calculated by a special form of Taylor double series for ${\bf\tau}_{1}(\boldsymbol{x}-\boldsymbol{y})$ and $\unicode[STIX]{x1D642}_{1}(\boldsymbol{x}-\boldsymbol{y})$ in powers of $\boldsymbol{x}-\boldsymbol{x}_{{\it\gamma}}^{o,min}$ and $\boldsymbol{y}-\boldsymbol{x}_{{\it\delta}}^{o}$ (Zinchenko & Davis Reference Zinchenko and Davis2000), which only requires knowledge of $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ , the associated vector of pressures $\boldsymbol{q}_{1}(\boldsymbol{r})$ and their high-order derivatives at $\boldsymbol{r}=\boldsymbol{R}_{{\it\gamma}{\it\delta}}=\boldsymbol{x}_{{\it\delta}}^{o}-\boldsymbol{x}_{{\it\gamma}}^{o,min}$ (appendix B of paper I). The optimized form of far-field expansions (appendix A of Zinchenko & Davis Reference Zinchenko and Davis2003) valid for any lattice geometry (with a unit cell volume) is used herein in favour of the algorithm from paper I. The necessary functions $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r}),\boldsymbol{q}_{1}(\boldsymbol{r})$ and their derivatives are pretabulated (prior to iterations of (2.8)) to a sufficient order $n_{tab}(n_{1},n_{2},n_{3})$ on a skewed uniform mesh $\boldsymbol{x}=hn_{1}\boldsymbol{e}_{1}+hn_{2}\boldsymbol{e}_{2}+hn_{3}\boldsymbol{e}_{3}$ , with $h=0.5/N_{T}$ and the integers $n_{i}$ in the range $|n_{1}|\leqslant N_{T}$ , $0\leqslant n_{2}$ , $n_{3}\leqslant N_{T}$ . Extension to negative $n_{2}$ or $n_{3}$ simply follows from the symmetry properties of $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ and $\boldsymbol{q}_{1}(\boldsymbol{r})$ . Since the singularities of $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ and $\boldsymbol{q}_{1}(\boldsymbol{r})$ are far away from the origin, these functions are very smooth in the tabulation range $V_{tab}=\{{\it\xi}^{i}\boldsymbol{e}_{i},|{\it\xi}^{i}|\leqslant 1/2\}$ , so only a small value of $N_{T}$ (typically 7) and a small order $n_{tab}$ of derivatives suffice. Still, due to the lattice time dependence, the tabulation is required at each time step. The most efficient Ewald-like forms for $\unicode[STIX]{x1D642}_{1}$ and $\boldsymbol{q}_{1}$ follow from Hasimoto (Reference Hasimoto1959), like in paper I. These forms, together with fast high-order differentiation techniques (appendix C of paper I), allow dynamical tabulation of $\unicode[STIX]{x1D642}_{1},\boldsymbol{q}_{1}$ and their high-order derivatives at a negligible cost. For a skewed lattice, the minimal vector $\boldsymbol{R}_{{\it\gamma}{\it\delta}}$ is not always in the tabulation range $V_{tab}$ . In general, the periodic shift $\boldsymbol{P}(\boldsymbol{m}^{\ast })$ is found to map $\boldsymbol{R}_{{\it\gamma}{\it\delta}}$ on $V_{tab}$ . The necessary values of $\unicode[STIX]{x1D642}(\boldsymbol{r})$ , $\boldsymbol{q}_{1}(\boldsymbol{r})$ and their high-order derivatives at $\boldsymbol{r}=\boldsymbol{R}_{{\it\gamma}{\it\delta}}$ are easily linked to those at $\boldsymbol{r}=\boldsymbol{R}_{{\it\gamma}{\it\delta}}^{\ast }=\boldsymbol{R}_{{\it\gamma}{\it\delta}}+\boldsymbol{P}(\boldsymbol{m})\in V_{tab}$ , by adding/subtracting a few free-space contributions in an efficient manner (appendix D of paper I). In turn, the values of $\unicode[STIX]{x1D642}_{1},\boldsymbol{q}_{1}$ and their derivatives at $\boldsymbol{r}=\boldsymbol{R}_{{\it\gamma}{\it\delta}}^{\ast }$ are obtained from the nearest table node by Taylor expansions (typically of second or third order). Again, due to the smoothness of ${\bf\tau}_{1}(\boldsymbol{r})$ and $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ near $\boldsymbol{r}=\boldsymbol{R}_{{\it\gamma}{\it\delta}}$ , the far-field double series in powers of $\boldsymbol{x}-\boldsymbol{x}_{{\it\gamma}}^{min}$ and $\boldsymbol{y}-\boldsymbol{x}_{{\it\delta}}^{o}$ for the ( ${\bf\tau}_{1},\unicode[STIX]{x1D642}_{1}$ ) contribution to (3.3) and (3.4) are rapidly convergent, and just a few terms suffice.

If the shell $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ around $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ and the shell $\mathscr{D}_{{\it\delta}}$ have ‘enough clearance’ (see below), the free-space contributions from $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ to (3.3) and (3.4) are evaluated at $\boldsymbol{y}\in \mathscr{B}_{{\it\delta}}$ by first re-expanding Lamb’s singular series for $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ as Lamb’s regular series near $\boldsymbol{x}_{{\it\delta}}^{o}$ . All such re-expansions are greatly accelerated (Zinchenko 1994, Zinchenko & Davis Reference Zinchenko and Davis2000) by rotational transformations of spherical harmonics with Wigner functions, and the cumulative Lamb’s regular series from all blocks ‘sufficiently separated’ from $\mathscr{B}_{{\it\delta}}$ are then efficiently evaluated at all mesh nodes $\boldsymbol{y}\in \mathscr{B}_{{\it\delta}}$ .

If the block $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ is not ‘sufficiently separated’ from $\mathscr{B}_{{\it\delta}}$ , but node $\boldsymbol{y}\in \mathscr{B}_{{\it\delta}}$ is ‘well outside’ the shell $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ (as defined below), Lamb’s singular series for $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ are used directly to compute the contributions of $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ to (3.3) and (3.4). Only in rare cases, when $\boldsymbol{y}$ is too close to $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ from outside, or is inside $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ , should we use direct summations for these contributions.

A separate scheme (also different from the one in paper I) is used for self-interactions, i.e. to calculate the sums in the first lines of (3.5) and (3.6). It was observed in our PE and strong mixed-flow simulations in near-critical conditions that one (or a few) drop(s) can acquire extreme elongation (even comparable with the size of the periodic cell), still without necking, while all the other drops remain compact. To handle such cases in a simple and efficient manner, we replace, for self-interactions only, the uniform partition (3.7) by a variable partition with the parameter $R_{sub}^{{\it\alpha}}$ (instead of $R_{sub}$ ) now dependent on the drop elongation. Drop partitioning into blocks is not used in the self-interactions (3.5) and (3.6); we work with entire drops instead. Let ${\it\zeta}_{o}$ be the minimum distance from the origin to the lattice points $\boldsymbol{P}(\boldsymbol{m}),\boldsymbol{m}\neq \mathbf{0}$ , and $d_{{\it\alpha}}$ be the radius of the minimal spherical shell around $S_{{\it\alpha}}$ centred at $\boldsymbol{x}_{{\it\alpha}}^{c}$ . The convergence rate of the far-field expansion (with ${\bf\tau}_{1}$ and $\unicode[STIX]{x1D642}_{1}$ in place of $\widetilde{{\bf\tau}}$ and $\unicode[STIX]{x1D642}$ in (3.5) and (3.6)) is conservatively measured by the progression exponent $2d_{{\it\alpha}}/\max (R_{sub}^{{\it\alpha}},{\it\zeta}_{o})$ . We require this quantity to stay within $(1+2a/{\it\zeta}_{o})/2$ , which determines $R_{sub}^{{\it\alpha}}$ . Upon the partition of $\unicode[STIX]{x1D642}$ and $\widetilde{{\bf\tau}}$ , the remaining free-space contributions from surfaces $S_{{\it\alpha}}+\boldsymbol{P}(\boldsymbol{m})$ (with $0<\boldsymbol{P}(\boldsymbol{m})\leqslant R_{sub}^{{\it\alpha}}$ ) to (3.5) and (3.6) are handled by direct summations. For an extremely deformed drop, a fairly large number of free-space contributions need to be handled in this manner; still, the procedure remains fast compared with the rest of the code, since it is used for self-interactions only.

In BI simulations of emulsion flow through a packed bed (Zinchenko & Davis Reference Zinchenko and Davis2008, Reference Zinchenko and Davis2013) with ultrahigh surface resolution $(N_{\triangle }>10\,000)$ on solid boundaries, it was found to be useful to additionally divide each surface into a large number of ‘patches’ for multipole acceleration. This second level of mesh-node decomposition would not be beneficial in the present simulations with typically $N_{\triangle }\leqslant 3840$ .

Economical truncation bounds. As in paper I, the vital feature of the present multipole-accelerated scheme is the ‘economical truncation’ of multipole expansions/re-expansions to greatly boost the performance. It is taken into account, in nearly optimal fashion, that the rate of convergence of Lamb’s free-space singular series for block $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ is a strong function of the clearance from the surrounding minimal spherical shell $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ . Likewise, the rate of convergence of singular-to-regular free-space re-expansion between $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ and $\mathscr{B}_{{\it\delta}}$ strongly depends on the clearance between the surrounding shells $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ and $\mathscr{D}_{{\it\delta}}$ . In constructing the economical truncation bounds, we follow the technique of paper I, based, in part, on plausible estimations for the rate of decay of multipole coefficients. The only difference is that, according to the partition (3.7), there are now more free-space interactions qualified as ‘near-field’. Accordingly, for blocks $\mathscr{B}_{{\it\gamma}}^{min}$ and $\mathscr{B}_{{\it\delta}}$ separated (centre-to-centre) by $\boldsymbol{R}_{{\it\gamma}{\it\delta}}$ , the appropriate distance to affect the rate of convergence of far-field series contributions $\mathscr{B}_{{\it\gamma}}^{min}\rightarrow \mathscr{B}_{{\it\delta}}$ is now ${\it\zeta}_{{\it\alpha}{\it\beta}}=\min |\boldsymbol{R}_{{\it\gamma}{\it\delta}}+\boldsymbol{P}(\boldsymbol{m})|$ (over $|\boldsymbol{P}(\boldsymbol{m})|>R_{sub}$ ). Our truncation scheme, inevitably semi-intuitive, contains a single heuristic ‘precision parameter’ ${\it\varepsilon}$ . Most importantly, as ${\it\varepsilon}\rightarrow 0$ , all multipoles in near-field expansions/re-expansions and all terms in the far-field Taylor series are eventually included. However, for optimal performance, a threshold order $k_{o}$ is set to limit the use of multipoles in near-field operations. For a given ${\it\varepsilon}$ , blocks $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ and $\mathscr{B}_{{\it\delta}}$ are qualified as ‘sufficiently separated’ if the multipole truncation order required for the block-to-block re-expansion is within $k_{o}$ . Likewise, the node $\boldsymbol{y}$ is considered to be ‘well outside’ the shell $\mathscr{D}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})$ if the necessary truncation order for Lamb’s singular series for $\mathscr{B}_{{\it\gamma}}^{min}+\boldsymbol{P}(\boldsymbol{m})\rightarrow \boldsymbol{y}$ free-space contribution does not exceed $k_{o}$ . The value of $k_{o}=25$ was about optimal in the present calculations. In the limit ${\it\varepsilon}\rightarrow 0$ , though, the same non-multipole solution is reached regardless of $k_{o}$ and other artificial parameters $R_{sub},R_{sub}^{{\it\alpha}}$ . The correctness and high efficiency of our multipole-accelerated code are demonstrated in § 5.

7.1. Results for $c=0.45$

Figure 13(a–d) summarizes our calculations for the non-dimensional steady state viscometric functions – the excess shear viscosity $\langle {\it\mu}_{sh}-1\rangle$ and normal stress differences $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ (defined in (2.2a )) at 45 % drop volume fractions and different viscosity ratios ${\it\lambda}$ . Resolutions $N_{\triangle }=2160{-}3840$ and system sizes $N=100{-}200$ were used, as specified in the caption. For ${\it\lambda}=10$ (figure 13 a) and $\mathit{Ca}=0.06{-}0.02$ , the runs were stopped at non-dimensional times (i.e. strains) $t=60{-}66$ , to average ${\it\mu}_{sh}-1$ , $N_{1}^{sh}$ and $N_{2}^{sh}$ with estimated absolute errors of 0.01–0.02. As an example of the difficulties of such ${\it\lambda}=10$ simulations, the run $N=100$ , $N_{\triangle }=3840$ for the smallest $\mathit{Ca}=0.02$ in figure 13(a) took 160 000 time steps with ${\approx}10$ BI iterations per time step. For the same ${\it\lambda}$ and higher $\mathit{Ca}=0.08$ , 0.10 and 0.16, the runs could not be continued beyond $t=46$ , 42 and 27 respectively, due to divergence of the BI iterations. Obviously, our variational technique (§ 3.1) of the double-layer near-singularity subtraction, although highly beneficial, may be not powerful enough for concentrated emulsions at extreme viscosity ratios and large deformations. These shorter runs still gave the averaging errors likely within 0.01–0.02, since the increase of $\mathit{Ca}$ reduces the statistical fluctuations. For all $\mathit{Ca}$ in figure 13(a), relaxation to the statistical steady state is fast, and cut-out times $t_{cut}=2{-}12$ sufficed. The results in figure 13(a) show significant shear thinning at the smallest $\mathit{Ca}$ . The first normal stress difference $\langle N_{1}^{sh}\rangle$ , smaller but comparable with $\langle {\it\mu}_{sh}-1\rangle$ , is also a strong function of $\mathit{Ca}$ when $\mathit{Ca}\ll 1$ , but it levels off at larger drop deformations. The negative second normal stress difference is small compared with $\langle {\it\mu}_{sh}-1\rangle$ .

Figure 13. Steady-state viscometric functions at $c=0.45$ for ${\it\lambda}=10$ (a), 3 (b), 1 (c) and 0.25 (d). Squares, $N_{\triangle }=2160$ ; diamonds, $N_{\triangle }=2880$ ; circles, $N_{\triangle }=3840$ . Solid symbols are for $N=100$ , open symbols are for $N=200$ . Small-size symbols are used for $N_{1}^{sh}$ to distinguish them from ${\it\mu}_{sh}-1$ values. Lines connect the data points selected for constitutive modelling in §§ 8 and 9.

When ${\it\lambda}=3$ (figure 13 b), relaxation to the statistical steady state remains fast for all capillary numbers considered. We used $t_{cut}=2.5{-}5$ to exclude the transient regimes and integrated to $t\sim 100{-}300$ in each run to reduce likely statistical errors in $\langle {\it\mu}_{sh}-1\rangle$ , $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ to $\pm 0.005{-}0.01$ . Divergence of BI iterations was never encountered at this ${\it\lambda}$ . At the smallest $\mathit{Ca}=0.02$ (with the strongest surface resolution effect), the results were obtained for both $N_{\triangle }=2160$ and 3840:

(7.1a−c ) $$\begin{eqnarray}\langle {\it\mu}_{sh}-1\rangle =2.539,\quad \langle N_{1}^{sh}\rangle =0.540,\quad \langle N_{2}^{sh}\rangle =-0.446\quad (N_{\triangle }=2160),\end{eqnarray}$$

(7.2a−c ) $$\begin{eqnarray}\langle {\it\mu}_{sh}-1\rangle =2.594,\quad \langle N_{1}^{sh}\rangle =0.534,\quad \langle N_{2}^{sh}\rangle =-0.449\quad (N_{\triangle }=3840).\end{eqnarray}$$

The deviations between $\langle {\it\mu}_{sh}-1\rangle$ are only 2 % and are negligible for $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ . The triangulation effects should be much smaller for the other $\mathit{Ca}$ values in figure 13(b). The qualitative trends of the viscometric functions are roughly the same as for ${\it\lambda}=10$ .

For ${\it\lambda}=1$ (figure 13 c), relaxation to the statistical steady state is fast ( $t_{cut}=3{-}4$ ) in all cases considered, except for $\mathit{Ca}=0.025$ , $N=200$ when $t_{cut}=45$ had to be excluded from the averaging. Even with fast relaxation in other runs, the statistical fluctuations grow dramatically as $\mathit{Ca}=0.025$ is approached. These are signs of configurational ‘phase transition’ to a partially ordered state. The final simulation times were $T\sim 100$ for all $\mathit{Ca}\geqslant 0.11$ , but 400–900 for all $\mathit{Ca}\leqslant 0.065$ . For the most difficult case $\mathit{Ca}=0.025$ , the resolution and system-size effects on the steady-state viscometric functions were explored:

(7.3) $$\begin{eqnarray}\displaystyle & \langle {\it\mu}_{sh}-1,N_{1}^{sh},N_{2}^{sh}\rangle =1.591,0.462,-0.269\quad (N=100,N_{\triangle }=2160), & \displaystyle\end{eqnarray}$$

(7.4) $$\begin{eqnarray}\displaystyle & \langle {\it\mu}_{sh}-1,N_{1}^{sh},N_{2}^{sh}\rangle =1.549,0.477,-0.226\quad (N=100,N_{\triangle }=3840), & \displaystyle\end{eqnarray}$$

(7.5) $$\begin{eqnarray}\displaystyle & \langle {\it\mu}_{sh}-1,N_{1}^{sh},N_{2}^{sh}\rangle =1.502,0.488,-0.233\quad (N=200,N_{\triangle }=3840). & \displaystyle\end{eqnarray}$$

$\langle {\it\mu}_{sh}-1\rangle$

$\mathit{Ca}=0.025$

$N=200$

$x_{1}$

$x_{2}$

$x_{3}$

$\mathit{Ca}\geqslant 0.0375$

The case of low viscosity ratio ${\it\lambda}=0.25$ (figure 13 d) stands out, by far, as the most difficult, since anomalies are observed in a wide range of capillary numbers. The simulation times were $T\sim 60{-}80$ for $\mathit{Ca}\geqslant 0.24$ , but much larger ( $T=200{-}800$ ) for all smaller $\mathit{Ca}$ . Rather fast relaxation ( $t_{cut}\sim 10$ ) to the statistical steady state for $\mathit{Ca}\geqslant 0.24$ dramatically slows down with $\mathit{Ca}$ reduced to 0.16, and it remains slow ( $t_{cut}\sim 100{-}200$ ) for all smaller $\mathit{Ca}$ . To illustrate this point, figure 14 presents the trajectories of ${\it\mu}_{sh}-1$ , $N_{1}^{sh}$ and $N_{2}^{sh}$ for (a)  $\mathit{Ca}=0.16$ , (b)  $\mathit{Ca}=0.12$ and (c,d)  $\mathit{Ca}=0.1$ , all at $N=100$ and $N_{\triangle }=2160$ . For (c,d), two different well-mixed initial configurations of spherical drops were used; to avoid overlapping of data, the trajectories of ${\it\mu}_{sh}-1$ are shown as insets. By excluding the transient parts ( $t_{cut}=125$ , 150, 230 and 100 for (a–d) respectively), excellent averaging quality is achieved, with the absolute errors in $\langle {\it\mu}_{sh}-1\rangle$ , $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ within 0.004 for (a), and within 0.002 for (b–d). However, the two $\mathit{Ca}=0.1$ runs clearly belong to two different ergodic classes. While the effect of the initial configuration on $\langle {\it\mu}_{sh}-1\rangle$ is small (0.675 in c versus 0.694 in d), it is more noticeable for $\langle N_{2}^{sh}\rangle$ ( $-0.202$ versus  $-0.253$ ) and particularly large for $\langle N_{1}^{sh}\rangle$ (0.594 versus 0.747). In contrast, we did not find any appreciable effects of the initial drop configuration on the time averaging for $\mathit{Ca}=0.12$ and 0.16 at $N=100$ . Ergodic difficulties akin to those at $\mathit{Ca}=0.1$ are familiar in phase transition simulations for molecular systems by molecular dynamics or Monte Carlo methods. In principle, increasing the system size should merge separate ergodic classes inaccessible to each other, and make average results independent of the initial conditions. In practice, though, such system-size increase can make the averaging task much harder, if not insurmountable, since the system must now transition with enough frequency between metastable quasiergodic classes for good statistics. For $N=200$ , we successfully obtained steady-state viscometric functions at $\mathit{Ca}=0.04$ , 0.08, 0.1, 0.12, 0.14, 0.16 and 0.2 (although with larger statistical uncertainties of $\pm 0.01{-}0.02$ in $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ ) to compare with those for $N=100$ . The only large differences were found for $\langle N_{1}^{sh}\rangle$ at $\mathit{Ca}=0.1$ and 0.12. From this analysis, $\mathit{Ca}\approx 0.1{-}0.12$ is believed to be the left end of the phase transition interval (while the right end is estimated as $\mathit{Ca}\approx 0.18$ ). The excess viscosity $\langle {\it\mu}_{sh}-1\rangle$ in figure 13(d) is the least sensitive to $N$ and ergodic difficulties, and is reasonably accurate in the whole range of $0.04\leqslant \mathit{Ca}\leqslant 0.34$ considered. In particular, the kink of the viscosity curve near $\mathit{Ca}=0.18$ is not a numerical artefact.

Figure 14. Trajectories of the viscometric functions for $c=0.45$ , ${\it\lambda}=0.25$ , $N=100$ , $N_{\triangle }=2160$ at (a) $\mathit{Ca}=0.16$ , (b) $\mathit{Ca}=0.12$ and (c,d) $\mathit{Ca}=0.1$ .

Our results in figure 13 confirm the hypothesis (suggested in paper I, although with little support therein) that increasing ${\it\lambda}$ shifts the phase transition in the emulsion shear flow to smaller capillary numbers. Indeed, the related ergodic difficulties, so pronounced at ${\it\lambda}=0.25$ , are not felt at ${\it\lambda}=1$ until $\mathit{Ca}$ is reduced to 0.025, and these difficulties are absent altogether for ${\it\lambda}=3$ and 10 in the whole range of $\mathit{Ca}$ considered. One reason for this shift is the increase in the drop deformation in the flow with higher ${\it\lambda}$ (in the dilute limit, the trend would be the opposite for ${\it\lambda}>1$ ). This larger deformation helps to eliminate geometrical blockages and merge different ergodic classes. Another reason is, presumably, more efficient mixing in the flow with higher ${\it\lambda}$ acting the same way.

Figure 15. Steady-state viscometric functions at $c=0.55$ for ${\it\lambda}=10$ (a), 3 (b), 1 (c) and 0.25 (d). Squares, $N_{\triangle }=2160$ ; diamonds, $N_{\triangle }=2880$ ; circles, $N_{\triangle }=3840$ . Solid symbols are for $N=100$ , open symbols are for $N=200$ . Plus symbols are for $N=400$ with $N_{\triangle }=2160$ . In (c), ${\it\mu}_{sh}-2$ is presented (instead of ${\it\mu}_{sh}-1$ ) to avoid data overlapping. Small-size symbols are used for $N_{1}^{sh}$ to distinguish them from the viscosity values. Lines connect the data points selected for constitutive modelling in §§ 8 and 9.

The classical model of a liquid of molecules interacting as hard spheres is known to exhibit phase transition to the solid state in the interval $0.494\leqslant c\leqslant 0.545$ of volume fractions $c$ (e.g. Hoover & Ree Reference Hoover and Ree1968; Hansen & McDonald Reference Hansen and McDonald1976), with discontinuities of the first derivatives of pressure versus  $c$ at the end points. The present case of dynamical phase transition in emulsion shear flow is certainly different. Still, it is unexpected that we observe phase-transitional behaviour at significantly smaller volume fraction $c=0.45$ without non-hydrodynamical attractive forces between the drops.

7.2. Results for $c=0.55$

Our calculations of the steady-state non-dimensional viscometric functions $\langle {\it\mu}_{sh}\!-\!1\rangle$ , $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ at $c=0.55$ are summarized in figure 15(a–d). Again, at ${\it\lambda}=10$ (figure 15 a), proceeding to large simulation times is most prohibitive due to smaller time steps and a larger number of BI iterations. For this ${\it\lambda}$ , the runs $\mathit{Ca}=0.02$ and 0.04 were stopped at $t=63$ and 74 respectively. At larger $\mathit{Ca}=0.06$ and 0.08, though, eventual divergence of the BI iterations did not allow us to proceed beyond $t=35$ and 45 respectively, even with the resolution $N_{\triangle }=2880$ . Still, the absolute statistical errors for all the data in figure 15(a) are expected to be within $\pm 0.02$ . Relaxation to the statistical steady state was fast for all capillary numbers considered, and $t_{cut}=2{-}5$ was sufficient to exclude the transient stages from averaging. As well as for $c=0.45$ (figure 13 a), the ${\it\lambda}=10$ emulsion shows most pronounced shear thinning, but is also significantly more viscous with the drop volume fraction increased to 55 %.

Figure 16. Two trajectories of the excess shear viscosity ${\it\mu}_{sh}-1$ for $c=0.55$ , ${\it\lambda}=3$ , $\mathit{Ca}=0.02$ and $N=100$ , starting from the same random configuration of spherical drops. Thin line, $N_{\triangle }=2160$ ; thick line, higher resolution $N_{\triangle }=3840$ .

For ${\it\lambda}=3$ (figure 15 b), the steady-state results could not be extended to $\mathit{Ca}=0.02$ for an unexpected reason. In figure 16, the trajectories of the excess viscosity ${\it\mu}_{sh}-1$ are shown for two simulations with $\mathit{Ca}=0.02$ and $N=100$ . Both simulations start from the same random configuration of spherical drops, the thin (blue) line being for $N_{\triangle }=2160$ , the thick (red) line for $N_{\triangle }=3840$ . The two runs appear to quickly reach statistical steady states and closely agree, on average, up to $t\sim 65$ . After $t\sim 70$ , though, the lower-resolution run $N_{\triangle }=2160$ abruptly settles to a different state with much lower viscosity, while the run $N_{\triangle }=3840$ remains at the higher viscosity. We believe that the only explanation for this erratic behaviour and the unusually strong effect of the resolution is that the conditions are on the verge of phase transition, and a small numerical perturbation (change in $N_{\triangle }$ ) brings the system to a different quasiergodic class after a long time. Graphical analysis of drop configurations confirms that the $N_{\triangle }=2160$ run settles to a partially ordered state, unlike the other run. To make sure that such difficulties (compromising time averaging) do not appear for larger $\mathit{Ca}$ in figure 15(b), we performed detailed testing of the effects of the initial conditions and programme parameters on the steady-state functions $\langle {\it\mu}_{sh}-1\rangle$ , $\langle N_{1}^{sh}\rangle$ , $\langle N_{2}^{sh}\rangle$ at $\mathit{Ca}=0.04$ . Two-level behaviour (akin to that in figure 16 for $N_{\triangle }=2160$ ) was never observed in these test runs; the simulated times were $T=230{-}470$ , with always fast relaxation ( $t_{cut}=4$ ) to the statistical steady state. For one initial drop configuration, the $N_{\triangle }=2160$ results (solid squares in figure 15 b at $\mathit{Ca}=0.04$ ) are (3.254, 1.745, $-0.817$ ). Repeating the run with $N_{\triangle }=3840$ gave (3.307, 1.793, $-0.879$ ) (solid circles in figure 15 b). Tightening the ‘precision’ ${\it\varepsilon}$ to $10^{-4}a$ in the $N_{\triangle }=2160$ run yielded (3.267, 1.733, $-0.852$ ). Coarsening the barrier ${\it\delta}_{min}$ (which controls drop overlap, see § 4) to $0.003a$ resulted in (3.271, 1.743, $-0.846$ ). Finally, replacing the initial drop configuration in the $N_{\triangle }=2160$ run gave (3.269, 1.770, $-0.832$ ). The deviations from the base results (shown by solid squares in figure 15 b) are within 1.8 % of $\langle {\it\mu}_{sh}-1\rangle$ , the main drop-phase contribution to the stress, and often much less. Therefore (unlike for $\mathit{Ca}=0.02$ ), the steady-state results for $\mathit{Ca}\geqslant 0.04$ are stable to the programme parameters and do not suffer from ergodic difficulties. For $\mathit{Ca}\geqslant 0.06$ , shorter simulation times $T\sim 110{-}160$ were sufficient (with $t_{cut}=4{-}5$ ). Although the results in figure 15(b) (and 15 a) are for $N=100$ only, any significant $N$ dependence is not suspect; the latter could only result from phase transition and related ergodic difficulties.

Figure 17. The trajectories of the viscometric functions ${\it\mu}_{sh}$ , $N_{1}^{sh}$ and $-N_{2}^{sh}$ (top to bottom on each graph) for $c=0.55$ , ${\it\lambda}=1$ , $N=200$ and $N_{\triangle }=2160$ ; (a) $\mathit{Ca}=0.075$ , (b) $\mathit{Ca}=0.07$ , (c) $\mathit{Ca}=0.06$ , (d) $\mathit{Ca}=0.05$ .

When ${\it\lambda}$ is decreased to 1 (figure 15 c), the phase transition is shifted to much higher $\mathit{Ca}$ ; the ergodic difficulties in the transition range are severe and appear to be insurmountable, even when using $N=400$ drops. Figure 17(a–d) shows the trajectories of ${\it\mu}_{sh}(t)$ , $N_{1}^{sh}(t)$ and $-N_{2}^{sh}(t)$ at $N=200$ and $N_{\triangle }=2160$ , when $\mathit{Ca}$ is progressively decreased from (a) 0.075 to (d) 0.05. Here, we plot the full viscosity (instead of ${\it\mu}_{sh}-1$ ) to avoid data overlapping. Fast approach to the statistical steady state at $\mathit{Ca}=0.075~(t_{cut}=4)$ changes to very slow relaxation ( $t_{cut}\sim 1000{-}1200$ ) at $\mathit{Ca}=0.06$ and 0.07. Relaxation is still slow ( $t_{cut}=160$ ), although not as slow, at $\mathit{Ca}=0.05$ . The interval $\mathit{Ca}\approx 0.06{-}0.07$ is therefore believed to be phase-transitional. For $\mathit{Ca}\geqslant 0.075$ , very good averaging quality is achieved with simulation times $T\sim 200{-}600$ . The results are well convergent with the number of drops $N=100{-}200$ when $\mathit{Ca}\geqslant 0.08$ (as tested for $\mathit{Ca}=0.08{-}0.11$ ), but the system size has to be increased to $N=200{-}400$ at $\mathit{Ca}=0.075$ for convergence. For $\mathit{Ca}\leqslant 0.07$ , all the runs were extended to much larger times for adequate averaging, but the results still show a strong dependence on the initial drop configuration (as tested for $N=100$ , cf. solid squares in figure 15 c) and on the system size $N=100{-}400$ . Therefore, only the $\mathit{Ca}\geqslant 0.075$ portion of the results in figure 15(c) is trusted as unambiguous and is used for constitutive modelling in §§ 8 and 9.

It is very instructive to compare the transition from fast to slow relaxation in figure 17 with figure 12 for PE at the same $c=0.55$ and ${\it\lambda}=1$ . For shear flow (figure 17), we observe this transition to occur when the average drop length $\langle L_{d}\rangle$ in the steady state (correlating with the drop deformation and the degree of configurational arrest) decreases from $2.42a$ at $\mathit{Ca}=0.075$ to $2.36a$ at $\mathit{Ca}=0.07$ . In contrast, in the PE simulations (figure 12), $\langle L_{d}\rangle$ progressively decreases from $2.42a$ at $\mathit{Ca}=0.06a$ to much smaller, $2.15a$ , at $\mathit{Ca}=0.02$ , but the relaxation remains fast. This comparison reveals a profound difference between the two types of flow: in PE, phase transition does not occur regardless of how small $\mathit{Ca}$ is, at least for drop volume fractions $c\leqslant 0.55$ . The same is true for other viscosity ratios ${\it\lambda}$ .

The small viscosity ratio ${\it\lambda}=0.25$ at $c=0.55$ (figure 15 d) is another very challenging case in the shear-flow simulations. The transition from fast ( $t_{cut}\sim 10$ ) to slow ( $t_{cut}\sim 100$ ) relaxation is observed when $\mathit{Ca}\approx 0.19{-}0.22$ , which estimates the right end of the phase transition interval. The simulations for figure 15(d) were performed to $T=72{-}75$ at $\mathit{Ca}\geqslant 0.25$ , and to $T\sim 200{-}400$ at each smaller $\mathit{Ca}\leqslant 0.22$ . There is no doubt in the accuracy of all steady-state viscometric functions for $\mathit{Ca}>0.22$ . In the transition range, though, strong effects of the initial drop configuration and/or the system size $N=100{-}200$ on the average normal stress differences (especially $\langle N_{1}^{sh}\rangle$ ) are observed. In contrast, the excess viscosity $\langle {\it\mu}_{sh}-1\rangle$ is almost insensitive to these difficulties and is well convergent even in the transition range, at least for $\mathit{Ca}\geqslant 0.07$ . Again, the viscosity curve has a kink at $\mathit{Ca}\approx 0.19$ , near the upper end of the transition range.

Since the flow-induced partially ordered drop arrangement has an approximately hexagonal structure in the ( $x_{2}$ , $x_{3}$ ) plane, it would seem advantageous to adapt the shape of the periodic cell to this structure. Such algorithm modification, though, was not attempted in the present work for two reasons. The present version of the code heavily exploits the symmetry properties of the periodic Green function $\unicode[STIX]{x1D642}(\boldsymbol{r})$ about the plane $r_{3}=0$ , so an extension to an arbitrary geometry, with $\boldsymbol{e}_{2}(t)$ non-orthogonal to $\boldsymbol{e}_{3}$ , would not be straightforward. More importantly, while such an algorithm adaptation can alleviate ergodic and convergence difficulties in the simulations well below the transition point to order, it would not help near the transition with the coexistence of ordered and disordered states. Nor would such an adaptation help to identify the onset of transition. As an example, figure 18(a) shows a well-ordered flow achieved at a large strain $t=837$ in the simulation with $c=0.55$ , ${\it\lambda}=1$ , $N=400$ and $\mathit{Ca}=0.05$ . For $\mathit{Ca}=0.06$ , still below the right end of the transition interval, the emulsion is not nearly as ordered (figure 18 b).

Figure 18. Snapshots from the shear-flow simulations at $c=0.55$ , ${\it\lambda}=1$ , $N=400$ and $N_{\triangle }=2160$ with (a) $\mathit{Ca}=0.05$ , $t=837$ and (b) $\mathit{Ca}=0.06$ , $t=1155$ . The centres of $N$ independent drops have been mapped into the cube $(0,1)^{3}$ ; the same view is shown in (a) and (b).

7.3. The intermediate case $c=0.5$

At the intermediate volume fraction of $c=0.5$ , we have studied the steady-state viscometric functions for ${\it\lambda}=1$ only (figure 19). For moderate $\mathit{Ca}\geqslant 0.1$ , simulation times $T=80{-}190$ were sufficient to provide high averaging accuracy. As $\mathit{Ca}$ is decreased to 0.06–0.07, the statistical fluctuations sharply grow, necessitating much larger simulation times (up to $T\sim 3000$ ), but there is still no clear transient regime, so we used $t_{cut}=3{-}5$ for all $0.06\leqslant \mathit{Ca}\leqslant 0.18$ . For all $\mathit{Ca}\geqslant 0.07$ , the results are proven to be practically system-size-independent (compare the solid squares at $N=100$ with the open squares at $N=200$ at $\mathit{Ca}=0.07$ , 0.075 and 0.1; indistinguishable for $\langle N_{2}^{sh}\rangle$ ). At $\mathit{Ca}=0.06$ , the agreement between the $N=100$ and $N=200$ results is less favourable for $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ . At this $\mathit{Ca}$ , there is also a comparable and modest resolution effect when $N_{\triangle }$ is increased from 2160 (solid squares) to 3840 (solid circles). Pronounced slow relaxation to the statistical steady state is first observed when $\mathit{Ca}$ is decreased to 0.05. In the four runs for this $\mathit{Ca}$ , the simulation times were $T=800{-}1900$ with cutoffs $t_{cut}=160{-}400$ . Starting the $N=100$ , $N_{\triangle }=2160$ runs from two different initial configurations of spherical drops gave practically indistinguishable results. Increasing the resolution to $N_{\triangle }=3840$ (solid circles) increases the excess viscosity (from 2.55 to 2.58) and the magnitude of $\langle N_{2}^{sh}\rangle$ (from 0.32 to 0.36), but reduces $\langle N_{1}^{sh}\rangle$ (from 1.10 to 1.03). Increasing the system size to $N=200$ (open squares) has the opposite effect of slightly reducing the viscosity (which we associate with the phase-transitional behaviour) and further reducing $\langle N_{1}^{sh}\rangle$ ; the second normal stress difference is not affected.

Figure 19. Steady-state viscometric functions at $c=0.5$ and ${\it\lambda}=1$ . Squares, $N_{\triangle }=2160$ ; circles, $N_{\triangle }=3840$ . Solid symbols are for $N=100$ , open symbols are for $N=200$ . Plus symbols are for $N=100$ with $N_{\triangle }=6000$ . Small-size symbols are used for $N_{1}^{sh}$ . Lines connect the data selected for constitutive modelling in §§ 8 and 9.

In the three most difficult runs at $\mathit{Ca}=0.025$ , (i)  $N=100$ , $N_{\triangle }=3840$ , (ii)  $N=200$ , $N_{\triangle }=3840$ and (iii)  $N=100$ , $N_{\triangle }=6000$ , our simulation times were limited to $T=900$ , 500 and 400 respectively, with still acceptable quality of averaging. There are only modest variations, if any, in the excess viscosity $\langle {\it\mu}_{sh}-1\rangle$ between the runs, but serious misconvergence in the normal stress differences.

Overall, the accuracy of our viscosity calculations in figure 19 is sufficient to state the existence of a shallow local minimum of $\langle {\it\mu}_{sh}-1\rangle$ versus  $\mathit{Ca}$ near $\mathit{Ca}\approx 0.05{-}0.06$ due to phase transition, followed by increase in the viscosity as the capillary number is further reduced. It is more difficult to argue about the behaviour of $\langle N_{1}^{sh}\rangle$ in the transition region due to larger uncertainties. Only the $\mathit{Ca}\geqslant 0.07$ part of the results from figure 19 is used in the constitutive modelling of §§ 8 and 9.

7.4. Comparison with paper I and other work

Paper I was the first to present large-scale shear-flow simulations and steady-state viscometric functions for concentrated emulsions. In addition to the results for $c\leqslant 0.4$ (not relevant to the present work), $\langle {\it\mu}_{sh}\rangle$ , $\langle N_{1}^{sh}\rangle$ and $\langle N_{2}^{sh}\rangle$ were reported versus  $\mathit{Ca}$ for (i) $c=0.45$ and 0.5 at ${\it\lambda}=1$ and (ii) $c=0.45$ at ${\it\lambda}=3$ . The viscosity values were also reported versus  $\mathit{Ca}$ for $c=0.55$ at ${\it\lambda}=1$ and 3, but it was not possible to even estimate the steady-state normal stress differences in either case due to the short simulation times and large statistical fluctuations.

In addition to substantially new calculations for the difficult viscosity ratios of ${\it\lambda}=0.25$ and 10 at $c=0.45$ and 0.55, the present results for ${\it\lambda}=1$ and 3 stem from a far more systematic analysis, with much larger simulated strains for averaging, finer surface triangulations and (often) larger system size. In paper I with $N\sim 100$ , the strains were limited to 25–45 for ${\it\lambda}=1$ and only ${\sim}15{-}25$ for ${\it\lambda}=3$ , while they are at least of the order of hundreds for all ${\it\lambda}$ (and even thousands for ${\it\lambda}=1$ ) in the present work. For ${\it\lambda}=3$ , the surface resolutions in paper I were limited to $N_{\triangle }\leqslant 1500$ , compared with 2160–3840 used herein for all ${\it\lambda}\neq 1$ . Finally, in the present work the maximum system size was extended to $N=400$ for matching, and to $N=200$ for contrast viscosities, twice as large as in paper I.

For moderate $\mathit{Ca}$ (roughly when the average drop length $\langle L_{d}\rangle$ exceeds $2.6a$ ), the present results for ${\it\lambda}=1$ and 3 are fully consistent with those of paper I, within the statistical errors therein. In the most difficult range of small $\mathit{Ca}$ (and $c\geqslant 0.45$ ), the results of paper I somewhat overpredict the first normal stress difference and (when ${\it\lambda}=1$ ) the shear viscosity. The reason is twofold. First, Zinchenko & Davis (Reference Zinchenko, Davis and Petsev2004) used the same code as in paper I and in the same limited strain range, but with ultrahigh resolutions $N_{\triangle }=6000{-}8640$ to recalculate the ${\it\lambda}=1$ shear-flow results for $c=0.45{-}0.5$ at $\mathit{Ca}=0.025$ and 0.05. This change results, indeed, in noticeably smaller $\langle {\it\mu}_{sh}\rangle$ and $\langle N_{1}^{sh}\rangle$ . Presumably, tiny surface overlaps (allowed by the code of paper I) favour some drop clustering (especially at small $\mathit{Ca}$ and ${\it\lambda}=1$ versus ${\it\lambda}=3$ ), leading to higher emulsion viscosity and $\langle N_{1}^{sh}\rangle$ , unless this artefact is eliminated by ultrahigh surface resolution. The present code, with explicit surface overlap control, achieves much faster convergence for the viscometric functions with respect to surface triangulations at small $\mathit{Ca}$ . Second, the corrected viscosity values from Zinchenko & Davis (Reference Zinchenko, Davis and Petsev2004) are still slightly higher than the more accurate results herein (most noticeably, at $c=0.5$ , ${\it\lambda}=1$ and $\mathit{Ca}=0.025$ , their work gives $\langle {\it\mu}_{sh}\rangle =3.07$ versus the narrow range of 2.69–2.77 from our figure 19). This remaining discrepancy is obviously due to the relatively short simulation times in paper I and in Zinchenko & Davis (Reference Zinchenko, Davis and Petsev2004), so that averaging was, in fact, performed over metastable levels, well before the emulsion viscosity settles down to the true statistical steady state. Accordingly, the phase-transitional behaviour at $c=0.5$ and ${\it\lambda}=1$ could not be detected in those studies.

At $c=0.55$ and ${\it\lambda}=1$ , the viscosity calculations of paper I are indicative of phase transition, but the results for $\mathit{Ca}\leqslant 0.05$ suffer again from limited triangulations and statistical errors due to short simulation times. For $\mathit{Ca}\geqslant 0.1$ , good agreement with the present work is observed in the viscosity values, but paper I does not give the normal stress differences to compare with.

All of our shear-flow simulations in this section are for essentially subcritical conditions, away from drop breakup. When combining the extensiometric and viscometric functions for general constitutive modelling (as suggested in the next section), early breakup in planar-extensional flow is the limiting factor (except for small ${\it\lambda}$ ), which was the main reason not to try extending the present shear-flow simulations to breakup. Incidentally, the results of paper I for ${\it\lambda}=1$ and 3 come closer to breakup conditions, and they do not suffer in that range from limited triangulations or short simulation times.