3.1 Boundary-integral equations

Like in paper I, highly efficient multipole-accelerated boundary-integral techniques are used to solve the hydrodynamical multidrop problem at each time step, given the surface tension distributions on the drop surfaces. For either of the database flows, the length scale $L$ is defined as in § 2, the time scale is $(2I_{2})^{-1/2}$ and the velocity scale is $(2I_{2})^{1/2}L$ . The non-dimensional periodic cell $V$ has unit volume at all times. The surface tension is scaled with $\unicode[STIX]{x1D70E}_{eq}$ . The non-dimensional system of BI equations for the fluid velocity $\boldsymbol{u}(\boldsymbol{y})$ on drop surfaces $S_{\unicode[STIX]{x1D6FC}}$ to be solved iteratively is written as

(3.1) $$\begin{eqnarray}\displaystyle \boldsymbol{u}(\boldsymbol{y}) & = & \displaystyle \frac{(\unicode[STIX]{x1D706}-1)}{(\unicode[STIX]{x1D706}+1)}\left[2\mathop{\sum }_{\unicode[STIX]{x1D6FD}=1}^{N}\int _{S_{\unicode[STIX]{x1D6FD}}}\boldsymbol{Q}(\boldsymbol{x})\boldsymbol{\cdot }\widetilde{\unicode[STIX]{x1D749}}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}(\boldsymbol{x})\,\text{d}S_{x}+\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D6FC}}\right.\nonumber\\ \displaystyle & & \displaystyle +\,\left.2\mathop{\sum }_{\unicode[STIX]{x1D6FD}=1}^{N}\int _{S_{\unicode[STIX]{x1D6FD}}}(\boldsymbol{Q}\boldsymbol{\cdot }\boldsymbol{n})(\boldsymbol{x}-\boldsymbol{x}_{\unicode[STIX]{x1D6FD}}^{c})\,\text{d}S_{x}\right]+\boldsymbol{F}(\boldsymbol{y}),\quad \boldsymbol{y}\in S_{\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$

(Here and in the rest of § 3, the index $\unicode[STIX]{x1D6FD}$ should not be confused with the elasticity parameter (2.3).) Except for the expression for the inhomogeneous term $\boldsymbol{F}(\boldsymbol{y})$ (see below), (3.1) is identical to the one used in paper I for clean drops. In particular, velocity fluctuations $\boldsymbol{Q}(\boldsymbol{x})=\boldsymbol{u}(\boldsymbol{x})-\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D6FD}}$ from the surface-averaged values $\langle \boldsymbol{u}\rangle _{\unicode[STIX]{x1D6FD}}$ are used on $S_{\unicode[STIX]{x1D6FD}}$ instead of $\boldsymbol{u}(\boldsymbol{x})$ to significantly reduce the integrand variations and thereby accelerate convergence in the most difficult, multipole part of the algorithm. Instead of the fundamental stresslet $\unicode[STIX]{x1D749}(\boldsymbol{r})$ in the double-layer integrals (3.1), we work with the auxiliary third-rank tensor $\widetilde{\unicode[STIX]{x1D749}}(\boldsymbol{r})=\unicode[STIX]{x1D749}(\boldsymbol{r})-\boldsymbol{r}\unicode[STIX]{x1D644}$ which is (unlike $\unicode[STIX]{x1D749}$ ) triply periodic and fully symmetric in all three indices (Zinchenko & Davis Reference Zinchenko and Davis2000, Reference Zinchenko and Davis2002, Reference Zinchenko and Davis2003, paper I).

The non-dimensional form of $\boldsymbol{F}(\boldsymbol{y})$ in the present problem takes into account the variable surface tension and the Marangoni stress in the interfacial stress jump condition:

(3.2) $$\begin{eqnarray}\displaystyle \boldsymbol{F}(\boldsymbol{y})=\frac{2}{\unicode[STIX]{x1D706}+1}\left[\boldsymbol{u}_{\infty }(\boldsymbol{y})+\mathop{\sum }_{\unicode[STIX]{x1D6FD}=1}^{N}\int _{S_{\unicode[STIX]{x1D6FD}}}\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x}\right]. & & \displaystyle\end{eqnarray}$$

Here, the vector field

(3.3) $$\begin{eqnarray}\displaystyle \boldsymbol{f}(\boldsymbol{x})=\frac{a}{Ca}\{2[\unicode[STIX]{x1D70E}(\boldsymbol{x})k(\boldsymbol{x})-\langle \unicode[STIX]{x1D70E}k\rangle _{\unicode[STIX]{x1D6FD}}]\boldsymbol{n}(\boldsymbol{x})-\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D70E}\}\quad (\boldsymbol{x}\in S_{\unicode[STIX]{x1D6FD}}) & & \displaystyle\end{eqnarray}$$

represents the stress jump on $S_{\unicode[STIX]{x1D6FD}}$ (the surface average of $\unicode[STIX]{x1D70E}k$ over $S_{\unicode[STIX]{x1D6FD}}$ is subtracted to reduce the variation of $\boldsymbol{f}$ without affecting the integrals (3.2)), $a=a_{o}/L$ is the non-dimensional non-deformed drop radius and $\unicode[STIX]{x1D642}$ is the triply periodic, symmetric second-rank tensor of Hasimoto (Reference Hasimoto1959). With our normalization, this tensor satisfies

(3.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D642}(\boldsymbol{x})-\unicode[STIX]{x1D735}\boldsymbol{q}(\boldsymbol{x})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D749}(\boldsymbol{x})=\mathop{\sum }_{\boldsymbol{m}}\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{P}(\boldsymbol{m}))\unicode[STIX]{x1D644}, & & \displaystyle\end{eqnarray}$$

where $\boldsymbol{q}(\boldsymbol{x})$ is the vector of pressures for $\unicode[STIX]{x1D642}(\boldsymbol{x})$ , $\unicode[STIX]{x1D6FF}$ is the delta function 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}$ .

Once the interfacial velocity is found iteratively from (3.1), the macroscopic instantaneous stress tensor $\unicode[STIX]{x1D72E}$ for the emulsion can be expressed via surface integrals by Batchelor’s (Reference Batchelor1970) method. With $\unicode[STIX]{x1D72E}$ scaled by $\unicode[STIX]{x1D707}_{e}(2I_{2})^{1/2}$ , the relevant non-dimensional form becomes

(3.5) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D72E}=\text{I.T.}+2\unicode[STIX]{x1D640}_{\infty }+\mathop{\sum }_{\unicode[STIX]{x1D6FC}=1}^{N}\int _{S_{\unicode[STIX]{x1D6FC}}}[(\boldsymbol{x}-\boldsymbol{x}_{\unicode[STIX]{x1D6FC}}^{c})\boldsymbol{f}+(\unicode[STIX]{x1D706}-1)(\boldsymbol{u}\boldsymbol{n}+\boldsymbol{n}\boldsymbol{u})]\,\text{d}S, & & \displaystyle\end{eqnarray}$$

where I.T. is an insignificant isotropic part; tensor $\unicode[STIX]{x1D72E}$ is symmetric upon accurate surface discretization. From (3.5), all the rheological functions (2.5) and (2.6) are calculated.

3.2 Surface discretization, mesh control 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; in the present work, it was sufficient to use $N_{\triangle }=2160$ and $2880$ for systematic simulations, and $N_{\triangle }=3840$ for selected refinement and convergence testing. Initially, the drops are spheres with standard, nearly uniform triangulations. As the drops move and deform, the quality of triangulations is maintained by a combination of three different methods as detailed in paper I, namely (i) ‘passive mesh stabilization’ (in the form of Zinchenko & Davis Reference Zinchenko and Davis2013) through global minimization of the ‘kinetic energy’ function (separately on each surface), (ii) active mesh restructuring through minimization of the ‘potential energy’ function (ibid.) and (iii) node reconnection (edge swapping). The swapping strategy mostly follows Cristini, Blawzdziewicz & Loewenberg (Reference Cristini, Bławzdziewicz and Loewenberg2001). In the runs with modest drop deformation, the technique (i) alone, which maintains fixed mesh topology, has soft stability limitations on the time step and avoids any surface interpolations, was sufficient to maintain mesh quality for an indefinite time. Even for larger deformations, the power of passive mesh stabilization made the use of other tools (ii) and (iii) very infrequent, thus reducing the number of BI iterations – the main contributor to computational cost. More complex topological changes (node addition/subtraction) are avoided altogether in our simulations.

Surface normals and curvatures are calculated by the best paraboloid spline method of Zinchenko & Davis (Reference Zinchenko and Davis2000). Surface integrations (after removing/alleviating singularities and near singularities in the boundary integrals) 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.

Such desingularizations are paramount prior to numerical implementation of (3.1), especially at high drop volume fractions, because of the singular behaviour of the kernels $\unicode[STIX]{x1D642}(\boldsymbol{r})$ and $\widetilde{\unicode[STIX]{x1D749}}(\boldsymbol{r})$ at all lattice points $\boldsymbol{r}=\boldsymbol{P}(\boldsymbol{m})$ . According to (3.4),

(3.6a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D642}(\boldsymbol{r})\simeq \unicode[STIX]{x1D642}_{o}\left(\boldsymbol{r}-\boldsymbol{P}(\boldsymbol{m})\right),\quad \widetilde{\unicode[STIX]{x1D749}}(\boldsymbol{r})\simeq \unicode[STIX]{x1D749}_{o}\left(\boldsymbol{r}-\boldsymbol{P}(\boldsymbol{m})\right), & & \displaystyle\end{eqnarray}$$

where

(3.7a,b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D642}_{o}(\boldsymbol{r})=-\frac{1}{8\unicode[STIX]{x03C0}}\left(\frac{\unicode[STIX]{x1D644}}{r}+\frac{\boldsymbol{r}\boldsymbol{r}}{r^{3}}\right),\quad \unicode[STIX]{x1D749}_{o}(\boldsymbol{r})=\frac{3}{4\unicode[STIX]{x03C0}}\frac{\boldsymbol{r}\boldsymbol{r}\boldsymbol{r}}{r^{5}} & & \displaystyle\end{eqnarray}$$

are the free-space Stokeslet and the corresponding stresslet, respectively.

To this end, for the double-layer integrals (3.1), the approach of paper I is followed without changes. For self-interactions (when $\boldsymbol{y}\in S_{\unicode[STIX]{x1D6FC}}=S_{\unicode[STIX]{x1D6FD}}$ in (3.1)), complete desingularization is easily achieved. When the observation point $\boldsymbol{y}$ is away from, but close to, an integration surface $S_{\unicode[STIX]{x1D6FD}}$ (or its periodic image $S_{\unicode[STIX]{x1D6FD}}+\boldsymbol{P}(\boldsymbol{m})$ ), near-singular behaviour of the integrand in (3.1) is only alleviated, not eliminated completely. Nevertheless, our variational approach to the choice of the subtracted quantity $\boldsymbol{Q}^{\ast }$ greatly improves the near-singularity subtraction (compared to the standard choice) and allows one to overcome the divergence of BI iterations in a wide range of viscosity ratios $\unicode[STIX]{x1D706}$ , with successful simulations to large times. The tool of ‘barrier functions’ (paper I) facilitates efficient implementation of near-singularity subtractions, activating them only when necessary.

For single-layer integrals (3.2), the singularities and near singularities are weaker than for the double-layer counterparts. However, from our experience, (3.2) must be desingularized completely for successful long-time simulations at high emulsion concentrations. The approach of paper I is not applicable here, because of the tangential component in $\boldsymbol{f}$ due to the Marangoni stress. A more complex technique herein is based on the following general result.

Let $\boldsymbol{f}(\boldsymbol{x})=\boldsymbol{f}_{\Vert }(\boldsymbol{x})+f_{n}(\boldsymbol{x})\boldsymbol{n}(\boldsymbol{x})$ be an arbitrary vector field on a smooth closed surface $S$ , decomposed into the tangential ( $\boldsymbol{f}_{\Vert }$ ) and normal components ( $f_{n}=\boldsymbol{f}\boldsymbol{\cdot }\boldsymbol{n}$ ). Then, for an observation point $\boldsymbol{y}\in S$ or outside $S$ ,

(3.8) $$\begin{eqnarray}\displaystyle \int _{S}\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x} & = & \displaystyle \int _{S}[\,\boldsymbol{f}-(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{n}^{\ast })\boldsymbol{f}_{\Vert }^{\ast }-(\boldsymbol{f}_{\Vert }^{\ast }\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}^{\ast }-f_{n}^{\ast }\boldsymbol{n}]\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{S}[(\boldsymbol{x}-\boldsymbol{x}^{\ast })\boldsymbol{\cdot }\boldsymbol{n}^{\ast }]\,\boldsymbol{f}_{\Vert }^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{o}(\boldsymbol{x}-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}\,\text{d}S_{x}.\end{eqnarray}$$

Here, $\boldsymbol{x}^{\ast }=\boldsymbol{y}$ when $\boldsymbol{y}\in S$ , and $\boldsymbol{x}^{\ast }\in S$ is arbitrary when $\boldsymbol{y}$ is outside $S$ . On the right-hand side of (3.8), for brevity, $\boldsymbol{f}=\boldsymbol{f}(\boldsymbol{x})$ and $\boldsymbol{n}=\boldsymbol{n}(\boldsymbol{x})$ . The asterisk denotes the values of $\boldsymbol{n}$ , $\boldsymbol{f}_{\Vert }$ and $f_{n}$ at $\boldsymbol{x}^{\ast }$ . Our derivation of (3.8) in the appendix A borrows some elements from the regularization techniques of Klaseboer, Sun & Chan (Reference Klaseboer, Sun and Chan2012) for boundary-integral equations, but, in other respects, our simple logic is different and gives the new general formula (3.8) unavailable from their paper. In particular, (3.8) is unrelated to any BI equation.

Owing to $(\boldsymbol{x}-\boldsymbol{x}^{\ast })\boldsymbol{\cdot }\boldsymbol{n}^{\ast }=O(|\boldsymbol{x}-\boldsymbol{x}^{\ast }|^{2})$ at $\boldsymbol{x}\rightarrow \boldsymbol{x}^{\ast }$ , both integrals in the right-hand side of (3.8) are obviously non-singular when $\boldsymbol{x}^{\ast }=\boldsymbol{y}\in S$ . To make them non-singular when $\boldsymbol{y}\rightarrow S$ from outside, $\boldsymbol{x}^{\ast }$ is chosen as the mesh node on $S$ nearest to $\boldsymbol{y}$ .

Farutin, Biben & Misbah (Reference Farutin, Biben and Misbah2014) offered a quite different regularization of single-layer BI with a general force density $\boldsymbol{f}$ . Unlike our (3.8), however, their transformation of the left-hand side results in non-singular integrals only for $\boldsymbol{y}\in S$ ; when $\boldsymbol{y}$ is away from (but close to) $S$ , desingularization is not achieved. For this reason, (3.8) is the clear choice in the present simulations for highly concentrated emulsions, with strong hydrodynamical interactions between nearly touching drops.

Based on (3.6) and (3.8), the single-layer contribution in (3.2) from a surface $S_{\unicode[STIX]{x1D6FD}}\not \ni \boldsymbol{y}$ is calculated as

(3.9) $$\begin{eqnarray}\displaystyle & & \displaystyle \int _{S_{\unicode[STIX]{x1D6FD}}}\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x}\approx \mathop{\sum }_{\boldsymbol{x}_{j}\in S_{\unicode[STIX]{x1D6FD}}}\boldsymbol{f}(\boldsymbol{x}_{j})\boldsymbol{\cdot }\unicode[STIX]{x1D642}(\boldsymbol{x}_{j}-\boldsymbol{y})\triangle S_{j}\nonumber\\ \displaystyle & & \displaystyle \quad -\,\mathop{\sum }_{\boldsymbol{m}}\unicode[STIX]{x1D6E9}(\boldsymbol{m},\boldsymbol{y})\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{\unicode[STIX]{x1D6FD}}}[\boldsymbol{f}_{\Vert }^{\ast }(\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{n}^{\ast })+(\boldsymbol{f}_{\Vert }^{\ast }\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}^{\ast }+f_{n}^{\ast }\boldsymbol{n}]\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}_{j}+\boldsymbol{P}(\boldsymbol{m})-\boldsymbol{y})\triangle S_{j}\nonumber\\ \displaystyle & & \displaystyle \quad +\,\mathop{\sum }_{\boldsymbol{m}}\unicode[STIX]{x1D6E9}(\boldsymbol{m},\boldsymbol{y})\mathop{\sum }_{\boldsymbol{x}_{j}\in S_{\unicode[STIX]{x1D6FD}}}[(\boldsymbol{x}_{j}-\boldsymbol{x}^{\ast })\boldsymbol{\cdot }\boldsymbol{n}^{\ast }]\boldsymbol{f}_{\Vert }^{\ast }\boldsymbol{\cdot }\unicode[STIX]{x1D749}_{o}(\boldsymbol{x}_{j}+\boldsymbol{P}(\boldsymbol{m})-\boldsymbol{y})\boldsymbol{\cdot }\boldsymbol{n}\triangle S_{j}.\end{eqnarray}$$

Here, the summation is over all mesh nodes $\boldsymbol{x}_{j}\in S_{\unicode[STIX]{x1D6FD}}$ and all integer vectors $\boldsymbol{m}$ , and $\triangle S_{j}$ is the surface area associated with node $\boldsymbol{x}_{j}$ in the trapezoidal integration rule. The barrier function $\unicode[STIX]{x1D6E9}$ acts only when $\boldsymbol{y}$ is close to a periodic image $S_{\unicode[STIX]{x1D6FD}}+\boldsymbol{P}(\boldsymbol{m})$ of $S_{\unicode[STIX]{x1D6FD}}$ . Namely, $\unicode[STIX]{x1D6E9}\approx 1$ when the distance from $\boldsymbol{y}$ to the nearest node $\boldsymbol{x}^{\ast }+\boldsymbol{P}(\boldsymbol{m})$ on $S_{\unicode[STIX]{x1D6FD}}+\boldsymbol{P}(\boldsymbol{m})$ is much less than the threshold $h_{o}$ (set to $0.3a$ ), and $\unicode[STIX]{x1D6E9}$ gradually falls to zero as this distance approaches $h_{o}$ . The precise definition of $\unicode[STIX]{x1D6E9}$ (as in paper I) follows Zinchenko & Davis (Reference Zinchenko and Davis2002). For brevity, $\boldsymbol{n}=\boldsymbol{n}(\boldsymbol{x}_{j})$ in (3.9). The asterisk marks the values of $\boldsymbol{n}$ , $\boldsymbol{f}_{\Vert }$ and $f_{n}$ at $\boldsymbol{x}=\boldsymbol{x}^{\ast }\in S_{\unicode[STIX]{x1D6FD}}$ . This node $\boldsymbol{x}^{\ast }$ depends on $\boldsymbol{y}$ , $S_{\unicode[STIX]{x1D6FD}}$ and $\boldsymbol{m}$ . In the limit of fine triangulations $N_{\triangle }\rightarrow \infty$ , the additional terms (represented by the last two lines of (3.9)) total to zero. At finite $N_{\triangle }$ , though, the second line fully cancels the near-singular behaviour of the single-layer integrand in the first line, thus greatly accelerating the convergence to $N_{\triangle }=\infty$ ; without near-singularity subtractions in (3.9), our simulations could not succeed beyond very small times. For large enough systems of drops, not more than one integer vector $\boldsymbol{m}$ contributes to (3.9).

For self-interactions ( $\boldsymbol{y}\in S_{\unicode[STIX]{x1D6FC}}=S_{\unicode[STIX]{x1D6FD}}$ ), the necessary desingularization proceeds as

(3.10) $$\begin{eqnarray}\displaystyle \int _{S_{\unicode[STIX]{x1D6FC}}}\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x} & = & \displaystyle \int _{S_{\unicode[STIX]{x1D6FC}}}\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot }[\unicode[STIX]{x1D642}(\boldsymbol{x}-\boldsymbol{y})-\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}-\boldsymbol{y})]\,\text{d}S_{x}\nonumber\\ \displaystyle & & \displaystyle +\,\int _{S_{\unicode[STIX]{x1D6FC}}}\boldsymbol{f}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D642}_{o}(\boldsymbol{x}-\boldsymbol{y})\,\text{d}S_{x}.\end{eqnarray}$$

The first integral in the right-hand side of (3.10) is non-singular; the second one is desingularized using (3.8) with $\boldsymbol{x}^{\ast }=\boldsymbol{y}$ . After that, the trapezoidal rule is applied to all the integrals.

3.3 Multipole acceleration

Using direct point-to-point summations to compute all the boundary integrals in (3.1) and (3.2) (after desingularizations) would be prohibitively expensive, of course, in large dynamical simulations. This operation, required at each time step, has an $O(N^{2}N_{\triangle }^{2})$ cost, with a large coefficient due to complex nature of the kernels $\widetilde{\unicode[STIX]{x1D749}}$ and $\unicode[STIX]{x1D642}$ and BI iterations. We use instead, without changes, the multipole-accelerated scheme of paper I, which is a development from our earlier studies (Zinchenko & Davis Reference Zinchenko and Davis2000, Reference Zinchenko and Davis2002, Reference Zinchenko and Davis2003, Reference Zinchenko and Davis2008, Reference Zinchenko and Davis2013) on multidrop and multidrop/multiparticle interactions. This highly efficient, but involved technique detailed (and verified) in paper I and our earlier work, does not have to be explained herein at any length. Only a brief overview of the main features follows.

A minor portion of the boundary integrals containing the free-space kernels $\unicode[STIX]{x1D642}_{0}$ and $\unicode[STIX]{x1D749}_{o}$ (e.g. those appearing in the second and third lines of (3.9), and in the last line of (3.10) due to desingularizations) are inexpensively calculated by direct point-to-point summations. For the other boundary integrals, a hybrid multipole-accelerated scheme is used. Each drop, when strongly elongated, is divided into a few compact blocks; a compact drop coincides with its only block. A minimal shell is constructed around each block with sufficient accuracy. The periodic kernels $\unicode[STIX]{x1D642}(\boldsymbol{r})$ and $\widetilde{\unicode[STIX]{x1D749}}(\boldsymbol{r})$ are partitioned into the near-field and far-field parts (termed $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ and $\unicode[STIX]{x1D749}_{1}(\boldsymbol{r})$ for the far field) by subtracting the free-space contributions $\unicode[STIX]{x1D642}_{o}\left(\boldsymbol{r}+\boldsymbol{P}(\boldsymbol{m})\right)$ and $\unicode[STIX]{x1D749}_{o}\left(\boldsymbol{r}+\boldsymbol{P}(\boldsymbol{m})\right)$ from several (approximately ten) lattice vectors $\boldsymbol{P}(\boldsymbol{m})$ with small $|\boldsymbol{P}(\boldsymbol{m})|<1.4$ . The free-space contribution of each block is expanded as Lamb’s singular series about the block centre. Now, the near-field contribution from a block ${\mathcal{B}}_{\unicode[STIX]{x1D6FE}}$ can be calculated on another block ${\mathcal{B}}_{\unicode[STIX]{x1D6FF}}$ as a sum of free-space contributions from several periodic images of ${\mathcal{B}}_{\unicode[STIX]{x1D6FE}}$ . For each such image, this operation is either done through singular-to-regular re-expansions between Lambs’ series with subsequent pointwise calculation of the regular series on ${\mathcal{B}}_{\unicode[STIX]{x1D6FF}}$ , or by direct pointwise computation of Lamb’s singular series, or (in rare cases) by direct point-to-point summations, depending on which mode of calculation is appropriate/most efficient. Singular-to-regular re-expansions are based, for maximum efficiency, on rotational transformations of spherical harmonics by Wigner functions (Zinchenko Reference Zinchenko1994, Reference Zinchenko1998; Zinchenko & Davis Reference Zinchenko and Davis2000). These near-field operations essentially employ high-order multipoles when necessary. The far-field contribution from block ${\mathcal{B}}_{\unicode[STIX]{x1D6FE}}$ to block ${\mathcal{B}}_{\unicode[STIX]{x1D6FF}}$ is computed through the special form of Taylor double series near the block centres. Since the singularities of $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ and $\unicode[STIX]{x1D749}_{1}(\boldsymbol{r})$ are far away from the origin, these far-field expansions are always rapidly convergent and require just a few terms for each pair of blocks. The far-field operations are facilitated through a dynamically formed table of $\unicode[STIX]{x1D642}_{1}(\boldsymbol{r})$ , related pressures $\boldsymbol{q}_{1}(\boldsymbol{r})$ and their (relatively low-order) derivatives. This table, in turn, is efficiently computed using the Ewald-like form of the Green tensor $\unicode[STIX]{x1D642}$ which follows from Hasimoto (Reference Hasimoto1959).

Also, the vital feature of the multipole-accelerated scheme is the ‘economical truncation’ of multipole expansions/re-expansions to greatly speed up the performance. It is taken into account, in a nearly optimal way, that the rate of convergence of Lamb’s free-space singular series for a block strongly depends on the clearance from the surrounding minimal shell. Likewise, the convergence rate for free-space singular-to-regular re-expansion between two blocks strongly depends on the clearance between the surrounding minimal shells. For the far-field expansion between two blocks, the rate of (much faster) convergence approximately depends on the ratio of the minimal shell size around the blocks to the minimum distance between their periodic images not included in the near field. All these considerations lead to the construction of a broad spectrum of truncation bounds for near-field and far-field operations strongly dependent on an intuitive ‘precision parameter’ $\unicode[STIX]{x1D700}$ . As $\unicode[STIX]{x1D700}\rightarrow 0$ , the standard (much simpler) non-multipole solution of the multidrop BI problem (3.1)–(3.2) is approached, but the total computational cost is weakly sensitive to $\unicode[STIX]{x1D700}$ . So, the two solutions can be made practically indistinguishable, with the multipole-accelerated solution being two orders of magnitude faster. In the test dynamical simulations of paper I for clean strongly deformed drops with $c=0.55$ , $\unicode[STIX]{x1D706}=3$ , $N=100$ and $N_{\triangle }=2160$ interrupted at an advanced stage for comparisons with the direct summations method, the multipole-accelerated scheme was seen to perform one BI iteration of (3.1) at least 140 times faster. For less deformed drops (table 1 therein), this gain is even bigger (at least 230-fold). Similar dramatic acceleration was observed in the calculation of the inhomogeneous term $\boldsymbol{F}(\boldsymbol{y})$ ; using direct summations would not be an option in large-scale dynamical simulations even at $\unicode[STIX]{x1D706}=1$ when BI iterations are not required. The far-field part of our multipole acceleration scheme has an $O(N^{2})$ scaling, but with a relatively very small coefficient which decreases with the increase in the system size $N$ . There are also small parts of the coding that scale as $O(NN_{\triangle }^{2})$ . As long as the number $N$ of drops is within several hundred and the resolution $N_{\triangle }$ does not exceed several thousand elements per drop, the overall scaling of the algorithm remains approximately linear in both $N$ and $N_{\triangle }$ . Our approach to multipole acceleration has proved to be efficient in many dynamical multidrop and multidrop/multiparticle simulations; as discussed elsewhere (Zinchenko & Davis Reference Zinchenko and Davis2000, Reference Zinchenko and Davis2008), it is vastly different from the general fast multipole method, in both the logic and technical tools.

As in paper I, the iterative scheme to solve the BI problem (3.1) is the generalized minimal residual method, with the initial approximation constructed from two preceding time steps; for viscosity ratios $0.25\leqslant \unicode[STIX]{x1D706}\leqslant 3$ considered herein, typically a few iterations per time step suffice.

6.1 Planar mixed flow $\unicode[STIX]{x1D712}=0.16$

Figure 3 presents the final snapshots of our dynamical simulations at 60 % drop volume fraction and matching viscosities, with (a) $\unicode[STIX]{x1D6FD}=0.2$ , $Ca=0.0116$ , $N=100$ drops with centres in a periodic cell, $N_{\triangle }=3840$ triangular elements per drop and (b) $\unicode[STIX]{x1D6FD}=0.05$ , $Ca=0.0331$ , $N=200$ , $N_{\triangle }=2160$ . In both snapshots, the cube shown is not a periodic cell, but it has the same volume and the same number of drops with centres in it as does a periodic cell. For the large non-dimensional times shown (see § 3.1 for the definition), the systems are well at statistical steady state; there are no signs of drop ordering. Simulations with high concentrations $c$ and small drop deformations, like those in figure 3, provide the most stringent testing for the effects of surface resolution and drop jamming, owing to strongest hydrodynamical interactions and stress localization in small lubrication areas. In the run for figure 3(a), where the drops are nearly spherical with the average length $L_{d}\approx (2.20{-}2.22)a$ , only about 0.25 % of the total number of mesh nodes required the geometrical correction (§ 5) at each time step to prevent surface overlapping. In the run for figure 3(b), with somewhat larger drop deformation ( $\langle L_{d}\rangle \approx 2.30a$ ), the portion of the nodes requiring such correction was even less (0.1 %). Thus, surface contamination by surfactant appears to increase the thickness of lubrication layers due to Marangoni stresses. Also, the ability of our algorithm to resolve such layers quickly improves with further increase in the drop deformation.

Figure 3.

Final snapshots of the dynamical simulations for steady PM flow ( $\unicode[STIX]{x1D712}=0.16$ ) at $c=0.6$ and $\unicode[STIX]{x1D706}=1$ , with (a) $\unicode[STIX]{x1D6FD}=0.2$ , $Ca=0.0116$ , $t=352$ , $N=100$ , $N_{\triangle }=3840$ and (b) $\unicode[STIX]{x1D6FD}=0.05$ , $Ca=0.0331$ , $t=479$ , $N=200$ , $N_{\triangle }=2160$ .

Figure 4.

Trajectories of the rheological functions $\unicode[STIX]{x1D708}_{12}-1$ , $\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x1D708}_{22}$ (top to bottom) for the PM flow simulations corresponding to figures 3(a) and 3(b), respectively. Horizontal lines represent long-time averages.

Figures 4(a) and 4(b) show the trajectories of the rheological functions $\unicode[STIX]{x1D708}_{12}-1$ , $\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x1D708}_{22}$ (top to bottom), defined by (2.5), for the simulations corresponding to figures 3(a) and 3(b), respectively. Unity is subtracted from $\unicode[STIX]{x1D708}_{12}$ to exclude the continuous-phase contribution to the non-dimensional effective viscosity in the PM flow. Short relaxation times to the statistical steady state confirm the absence of phase transition to ordered flow in these simulations. Excluding the transient parts $t\leqslant t_{cut}$ (with $t_{cut}\approx 5$ for a and $t_{cut}\approx 10$ for b) of the trajectories, the long-time averages are accurately obtained:

(6.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\langle \unicode[STIX]{x1D708}_{12}-1\rangle =7.347\pm 0.008,\quad \langle \unicode[STIX]{x1D708}_{11}\rangle =1.330\pm 0.010,\\ \langle \unicode[STIX]{x1D708}_{22}\rangle =-1.688\pm 0.010\quad \text{for }(a)\end{array}\right\} & & \displaystyle\end{eqnarray}$$

and

(6.2) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\langle \unicode[STIX]{x1D708}_{12}-1\rangle =3.601\pm 0.005,\quad \langle \unicode[STIX]{x1D708}_{11}\rangle =1.373\pm 0.007,\\ \langle \unicode[STIX]{x1D708}_{22}\rangle =-1.356\pm 0.007\quad \text{for }(b).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

These averages are shown in figure 4 by the horizontal lines. The $\pm$ numbers in (6.1) and (6.2) represent estimated statistical errors. To test the effects of surface resolution, the run from figures 3(a) and 4(a) was repeated to $t=637$ using cruder triangulations $N_{\triangle }=2160$ , with the results

(6.3a-c ) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D708}_{12}-1\rangle =7.313\pm 0.006,\quad \langle \unicode[STIX]{x1D708}_{11}\rangle =1.251\pm 0.007,\quad \langle \unicode[STIX]{x1D708}_{22}\rangle =-1.625\pm 0.007. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

The discrepancies between (6.1) and (6.3) must be compared to $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ – the main drop-phase contribution to the deviatoric stress – and they are only 0.5 % for the excess viscosity $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ , and ${\approx}$ 1 % for $\langle \unicode[STIX]{x1D708}_{11}\rangle$ and $\langle \unicode[STIX]{x1D708}_{22}\rangle$ .

To assess the effects of the system size, the run from figures 3(b) and 4(b) was repeated with a smaller number of drops $N=100$ to $t\approx 1000$ , yielding

(6.4a-c ) $$\begin{eqnarray}\displaystyle \langle \unicode[STIX]{x1D708}_{12}-1\rangle =3.596\pm 0.005,\quad \langle \unicode[STIX]{x1D708}_{11}\rangle =1.386\pm 0.007,\quad \langle \unicode[STIX]{x1D708}_{22}\rangle =-1.374\pm 0.007. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Negligible differences between (6.2) and (6.4) confirm our general expectation that excessively large systems ( $N\sim 1000$ or more) are not needed in simulations of bulk rheology, as long as (i) periodic boundaries are implemented and (ii) there is no phase transition to ordered flow. Finally, the $N=100$ run was repeated with higher resolution $N_{\triangle }=3840$ to $t=409$ , to replace the numbers in (6.4) with $3.571\pm 0.008$ , $1.364\pm 0.010$ and $-1.374\pm 0.004$ . Here, the change in the resolution affects the results by only 0.6–0.7 % (relative to $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ ).

For $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.2$ , $N=100$ , $N_{\triangle }=2160$ and larger (but still subcritical) $Ca=0.058$ , the effects of the diffusivity parameter $\unicode[STIX]{x1D702}\ll 1$ in (4.7) on the steady-state rheological functions were also explored. To this end, additional simulations were performed for $\unicode[STIX]{x1D702}=5.7\times 10^{-5}$ and $5.7\times 10^{-4}$ to large times $t=260$ and $315$ , respectively, to exclude any appreciable averaging errors. However, such tenfold variation in $\unicode[STIX]{x1D702}$ (around the default value of $1.9\times 10^{-4}$ ) had only 0.5 % effect on $\langle \unicode[STIX]{x1D708}_{11}\rangle$ (relative to $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ ) and much smaller, undetectable effects on $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ and $\langle \unicode[STIX]{x1D708}_{22}\rangle$ . For comparison, the simulation with zero diffusion (as detailed in § 4) could only reach $t\approx 29.5$ , due to the surface tension becoming negative on one of the drops; repeating such simulation with much higher resolution $N_{\triangle }=5120$ made it far more expensive, but extended the time range only modestly, to $t\approx 40.5$ . It is hard to argue whether the difficulty is due to still insufficient resolution, or to the equation of state for $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D6E4})$ . Quite possibly, in real systems, abnormally small local surface tension will lead to subsequent cusping and tip streaming, with a miniscule amount of surfactant ejected from the tip, sufficient to restore normal values of the surface tension and stop surfactant leakage. Adequate simulations of tip streaming in concentrated multidrop systems would be very problematic. Instead, using a small but non-zero diffusion parameter $\unicode[STIX]{x1D702}$ is offered here as an acceptable and feasible way to bypass the local difficulties, because the hydrodynamical fluctuations around individual drops favouring abnormal surfactant accumulation (and thereby negative surface tension) in certain areas are very infrequent and very short-lived (for subcritical $Ca$ ), and so they do not affect the overall multidrop dynamics and global quantities of interest. Independence of our long-time results for the average stress of the diffusivity parameter $\unicode[STIX]{x1D702}\ll 1$ confirms the validity of the proposed approach. Thus, negative surface tensions are avoided, but catastrophically high surface curvatures still occasionally develop on individual drops (except for small $Ca$ ) requiring the pre-emptive smoothing procedure outlined in § 5. Figure 5 gives an example of such smoothing in our PM flow simulations. (It would not make sense, instead of smoothing, to try numerically resolving cusp shape irregularities in global three-dimensional simulations like ours, since the physical cusp length scale is known to be extremely small, even approaching molecular dimensions, e.g. Jeong & Moffatt Reference Jeong and Moffatt1992; Eggers & du Pont Reference Eggers and du Pont2009.) Of course, when using this artificial smoothing, only simulations with rare cusping can be accepted. For $Ca=0.058$ , there were only 11 cusping events in the whole simulation of 100 drops to $t=260$ when $\unicode[STIX]{x1D702}=5.7\times 10^{-5}$ ; cusping was not observed when $\unicode[STIX]{x1D702}=5.7\times 10^{-4}$ . In the two additional simulations to $t\approx 100$ for larger $Ca=0.0812$ (viewed as critical), there were 19 cusping events when $\unicode[STIX]{x1D702}=1.9\times 10^{-4}$ and $N_{\triangle }=2160$ , and more frequent cusping when $\unicode[STIX]{x1D702}=5.7\times 10^{-5}$ and $N_{\triangle }=5120$ . However, the global stress trajectories were still at statistical steady state allowing for time averaging, with less than 1 % difference in $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ , $\langle \unicode[STIX]{x1D708}_{11}\rangle$ and $\langle \unicode[STIX]{x1D708}_{22}\rangle$ (relative to $\langle \unicode[STIX]{x1D708}_{12}-1\rangle$ ) between the two simulations. The zero-diffusion run failed very early for this $Ca=0.0812$ . In all the cases mentioned, there was very close agreement between the global stress trajectories for the zero-diffusion solution and those for $\unicode[STIX]{x1D702}=5.7\times 10^{-5}$ to times $t=O(1)$ , until the trajectories randomize. This agreement between the two, vastly different ways to handle the surfactant transport confirmed technical correctness of both implementations.

Figure 5.

An example of cusp smoothing in the steady PM flow simulation with $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.2$ , $Ca=0.058$ , $N=100$ , $N_{\triangle }=2160$ and $\unicode[STIX]{x1D702}=5.7\times 10^{-5}$ . The fully remeshed drop shape in (b), used to continue the simulation, fits perfectly the original shape in (a), except for the very small area near the cusp of no global significance and a negligible effect on the total drop volume fraction. The view and scale are the same for (a) and (b).

Figure 6.

Trajectories of the rheological functions for steady PM flow ( $\unicode[STIX]{x1D712}=0.16$ ) at $c=0.55$ , $\unicode[STIX]{x1D706}=0.25$ , $\unicode[STIX]{x1D6FD}=0.1$ , $N=100$ and $N_{\triangle }=2880$ , with (a) $Ca=0.0696$ and (b) $Ca=0.0348$ .

Figure 6 presents the trajectories of our database rheological functions for a small viscosity ratio $\unicode[STIX]{x1D706}=0.25$ , moderate contamination $\unicode[STIX]{x1D6FD}=0.1$ and 55 % drop volume fraction, in the simulations with $N=100$ and $N_{\triangle }=2880$ . For both capillary numbers (a) $Ca=0.0696$ (with the average drop length $\langle L_{d}\rangle \approx 2.47a$ ) and (b) $Ca=0.0348$ (with $\langle L_{d}\rangle \approx 2.29a$ ), fast relaxation to statistical steady state is observed, confirming again the absence of phase transition to ordered flow. These behaviours are in stark contrast to those for clean-drop, shear-flow rheological functions in paper I at the same viscosity ratio $\unicode[STIX]{x1D706}=0.25$ and $N=100$ drops. Even for a smaller emulsion concentration $c=0.45$ and $Ca=0.16$ (figure 14a therein), very slow relaxation to statistical steady state was clearly observed (with the average drop deformation close to that in our figure 6 a), requiring at least $t_{cut}\approx 125$ to be excluded from averaging. Reducing $Ca$ to 0.1 (figure 14c of paper I) results in even slower relaxation ( $t_{cut}\approx 230$ ) to a state with drop deformation close to that in our figure 6(b). These pathologies in shear-flow, clean-drop simulations are characteristic of phase transition to partial order and were observed (paper I) in a broad range of flow intensities, so that the whole case $\unicode[STIX]{x1D706}=0.25$ had to be excluded therein from constitutive modelling for concentrated emulsions. In contrast, the present database PM flows withstand much more jamming without transition to order and are therefore well behaved and suited for general constitutive modelling, including the case $\unicode[STIX]{x1D706}=0.25$ (§§ 7 and 8).

Figures 7 and 8 summarize our systematic calculations for the steady-state (i.e. time-averaged) rheological functions $\unicode[STIX]{x1D708}_{12}-1$ , $\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x1D708}_{22}$ versus $Ca$ in the database PM flow with $\unicode[STIX]{x1D712}=0.16$ . (In what follows and where obvious, the angular brackets for time averaging are omitted.) Figure 7 covers two emulsion concentrations $c=0.45$ (a,c,e) and $c=0.55$ (b,d,f), and three viscosity ratios $\unicode[STIX]{x1D706}=3$ (a,b), $\unicode[STIX]{x1D706}=1$ (c,d) and $\unicode[STIX]{x1D706}=0.25$ (e,f), all at moderate surface contamination $\unicode[STIX]{x1D6FD}=0.1$ . Figure 8 presents the steady-state results at 60 % drop volume fraction and matching viscosities, with two surface contaminations (a) $\unicode[STIX]{x1D6FD}=0.2$ and (b) $\unicode[STIX]{x1D6FD}=0.05$ . The simulations were typically performed to non-dimensional times from one hundred to several hundred, always providing reliable time averaging, due to the absence of phase transition difficulties. Large simulation times are increasingly important at small $Ca$ , large emulsion concentrations and large viscosity ratios – all factors that increase statistical fluctuations. Since the results are well convergent with respect to surface triangulations in the most unfavourable case $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.2$ and $Ca=0.0116$ (see above) with the largest excess effective viscosity $\unicode[STIX]{x1D708}_{12}-1\approx 7.3$ , the surface resolutions used in all other runs for figures 7 and 8 are also deemed sufficient.

At high concentrations $c$ , the effective viscosity $\unicode[STIX]{x1D708}_{12}$ exhibits sharp tension thinning for small $Ca$ , while $\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x1D708}_{22}$ are mild functions of $Ca$ (and they approximately cancel each other, with small $\unicode[STIX]{x1D708}_{33}=-\unicode[STIX]{x1D708}_{11}-\unicode[STIX]{x1D708}_{22}$ ). The results in figures 7 and 8 support the views from paper I and Zinchenko & Davis (Reference Zinchenko and Davis2002) that a concentrated emulsion of deformable drops, with or without surfactant, is not a phenomenological Rivlin–Ericksen fluid. For the latter, described asymptotically by a hierarchy of Rivlin–Ericksen models (e.g. Newtonian fluid, second or fourth-order fluids, etc.) in the limit of slow flows (e.g. Astarita & Marrucci Reference Astarita and Marrucci1974), the PM flow viscosity $\unicode[STIX]{x1D708}_{12}$ would approach its limiting value at $Ca=0$ with an error of $O(Ca^{2})$ , and the diagonal elements $\unicode[STIX]{x1D708}_{ii}$ would all vanish in this limit. Such phenomenological predictions would contradict our results in figures 7 and 8. The reason (paper I) why a concentrated emulsion is not a Rivlin–Ericksen fluid is due to drop interactions, making small-deformation analysis a problem of singular perturbations, so that a regular drop shape expansion in $Ca\ll 1$ is not possible even in principle. In contrast, dilute emulsions of single drops behave, of course, like a Rivlin–Ericksen fluid in PM flow (figure 4 of Martin et al. Reference Martin, Zinchenko and Davis2014).

Our emulsion viscosity for the highest drop volume fractions $c$ , with strong thinning at $Ca\ll 1$ (most noticeable in figure 8 b), still does not point to the existence of the yield stress (which would require an even stronger behaviour $\unicode[STIX]{x1D708}_{12}\sim Ca^{-1}$ ); such a yield stress is assumed by some phenomenological models (e.g. Casson model). It is hard to speculate at this point whether the absence of the yield stress in our emulsions is due to still insufficiently high concentrations $c$ , or to the absence of strong colloidal interactions between the drops. For practical emulsions with $c>0.74$ , the existence of the yield stress was reported (Princen & Kiss Reference Princen and Kiss1989).

Figure 7.

Steady-state database rheological functions in the PM flow ( $\unicode[STIX]{x1D712}=0.16$ ) with $\unicode[STIX]{x1D6FD}=0.1$ and $N=100$ . Emulsion concentrations $c=0.45$ (a,c,e) and 0.55 (b,d,f). Viscosity ratios $\unicode[STIX]{x1D706}=3$ (a,b), 1 (c,d) and 0.25 (e,f). Squares, $N_{\triangle }=2160$ ; diamonds, $N_{\triangle }=2880$ .

Figure 8.

Steady-state database rheological functions in the PM flow ( $\unicode[STIX]{x1D712}=0.16$ ) for $c=0.6$ , $\unicode[STIX]{x1D706}=1$ and $N=100$ , with (a) $\unicode[STIX]{x1D6FD}=0.2$ and (b) $\unicode[STIX]{x1D6FD}=0.05$ . Squares, $N_{\triangle }=2160$ ; diamonds, $N_{\triangle }=2880$ ; circles, $N_{\triangle }=3840$ .

Figure 9.

(a) The final snapshot of the steady PE simulation with $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.05$ , $Ca=0.096$ , $N=200$ and $N_{\triangle }=2160$ at $t=346$ (the directions of the coordinate axes are the same as in figure 3); (b) the corresponding trajectories of the rheological functions.

6.2 Planar extension

Figure 9(a) presents the final snapshot of a PE simulation at 60 % drop volume fraction and $\unicode[STIX]{x1D706}=1$ , with $\unicode[STIX]{x1D6FD}=0.05$ , $Ca=0.096$ , $N=200$ strongly deformable drops in a periodic cell and resolution $N_{\triangle }=2160$ per drop. Again, the box shown (with its sides parallel to the coordinate axes) is not a periodic cell, but it has the same number of drops in it and the same volume as does a periodic cell. The capillary number $Ca=0.096$ is close to critical, beyond which massive drop break-up would occur for a non-diffusive surfactant. Figure 9(b) presents the corresponding trajectories of the non-dimensional excess PE viscosity $\unicode[STIX]{x1D707}_{pe}-1$ and cross-difference $N_{pe}$ (defined through (2.6)) versus the non-dimensional time (scaled with $(2\dot{\unicode[STIX]{x1D6E4}})^{-1}$ ). The convergence of the long-time averages $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle =3.207\pm 0.005$ and $\langle N_{pe}\rangle =0.710\pm 0.005$ was confirmed by repeating this run with a smaller system size $N=100$ and different resolution $N_{\triangle }=2880$ to $t=270$ , yielding $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle =3.211\pm 0.005$ and $\langle N_{pe}\rangle =0.709\pm 0.005$ , in excellent agreement with the results for $N=200$ and $N_{\triangle }=2160$ .

Figure 10.

Steady-state database rheological functions in the PE flow with $\unicode[STIX]{x1D6FD}=0.1$ and $N=100$ (solid symbols). Emulsion concentrations $c=0.45$ (a,c,e) and 0.55 (b,d,f). Viscosity ratios $\unicode[STIX]{x1D706}=3$ (a,b), 1 (c,d) and 0.25 (e,f). Squares, $N_{\triangle }=2160$ ; diamonds, $N_{\triangle }=2880$ ; circles, $N_{\triangle }=3840$ . Open symbols in (f) are the clean-drop results from paper I.

Figure 11.

Steady-state database rheological functions in the PE flow for $c=0.6$ , $\unicode[STIX]{x1D706}=1$ and $N=100$ , with (a) $\unicode[STIX]{x1D6FD}=0.2$ and (b) $\unicode[STIX]{x1D6FD}=0.05$ . Squares, $N_{\triangle }=2160$ ; diamonds, $N_{\triangle }=2880$ ; circles, $N_{\triangle }=3840$ .

Figure 10(a–f) summarizes the time-averaged results for $\unicode[STIX]{x1D707}_{pe}-1$ and $N_{pe}$ versus $Ca$ at two emulsion concentrations $c=0.45$ (a,c,e) and $c=0.55$ (b,d,f), and three viscosity ratios $\unicode[STIX]{x1D706}=3$ (a,b), $\unicode[STIX]{x1D706}=1$ (c,d) and $\unicode[STIX]{x1D706}=0.25$ (e,f), all at moderate surface contamination $\unicode[STIX]{x1D6FD}=0.1$ . These calculations are complemented by figure 11 at $c=0.6$ and $\unicode[STIX]{x1D706}=1$ , with two surface contaminations (a) $\unicode[STIX]{x1D6FD}=0.2$ and (b) $\unicode[STIX]{x1D6FD}=0.05$ . For comparison, open symbols and connecting dash lines in figure 10(f) are for clean drops $\unicode[STIX]{x1D6FD}=0$ from paper I. All the runs in figures 10 and 11 were typically performed to $t$ from ${\sim}150$ to several hundred, with always fast relaxation to statistical steady state and reliable time averaging. (This averaging is generally easier for the PE than for the PM database flow.) For the data in figures 10(a–d) and 10(f) at the smallest $Ca=0.02$ (when the surface discretization effect is by far the strongest), two resolutions were used, $N_{\triangle }=2160$ and $N_{\triangle }=3840$ ; the results are compared in table 1. The discrepancies are within 0.6 % (relative to $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle$ ), and often smaller. For conditions closer to jamming ( $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.2$ and $Ca=0.01$ ), such discrepancies are expectedly larger (cf. squares for $N_{\triangle }=2160$ and circles for $N_{\triangle }=3840$ in figure 11 a), but they remain small (1.4 % for $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle$ and 2.1 % for $\langle N_{pe}\rangle$ ). Thus, the present simulations are still seen to be reasonably accurate for emulsions with quite strong drop interactions, when the effective emulsion viscosity is nine times that for the carrier fluid.

Table 1.

The effects of surface resolution $N_{\triangle }$ on the steady-state rheological functions $\unicode[STIX]{x1D707}_{pe}-1$ and $N_{pe}$ in the PE flow with 100 drops, $\unicode[STIX]{x1D6FD}=0.1$ , $Ca=0.02$ and various emulsion concentrations $c$ and viscosity ratios $\unicode[STIX]{x1D706}$ .

All the results in figures 3, 4 and 6–11 and table 1 were obtained with the default value of $\unicode[STIX]{x1D702}=1.9\times 10^{-4}$ in the diffusivity definition (4.7). In the PE simulations for figures 10(b), 10(c) and 10(f) with the intermediate $Ca=0.04$ , two other values $\unicode[STIX]{x1D702}=5.7\times 10^{-4}$ and $\unicode[STIX]{x1D702}=5.7\times 10^{-5}$ around the default value were also tested. Such tenfold variation in $\unicode[STIX]{x1D702}$ had only (0.5–0.8) % effect on the steady state $\unicode[STIX]{x1D707}_{pe}-1$ , and (0.5–1.3) % effect on $\langle N_{pe}\rangle$ (relative to $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle$ ). For near critical $Ca=0.1$ in figure 10(d), the effect of the tenfold diffusivity parameter variation on $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle$ was larger, but still small (1.9 %). Thus, the intermediate default value of $\unicode[STIX]{x1D702}=1.9\times 10^{-4}$ was deemed sufficient for all the simulations to globally represent the case of non-diffusing surfactant.

The following examples show again how beneficial it is to use a small, but finite diffusion parameter $\unicode[STIX]{x1D702}$ in multidrop simulations instead of neglecting diffusion altogether. At (still subcritical) $Ca=0.08$ in figure 10(d), the Lagrange-based scheme for purely non-diffusive surfactant transport could only allow simulations to reach $t\approx 46$ , in order to keep $\unicode[STIX]{x1D70E}$ positive everywhere. Unexpectedly, increasing the surface resolution $N_{\triangle }$ from 2160 to 5120 has limited the zero-diffusion solution even more, to $t<32$ (due to difficulties with a different drop; a likely explanation is a quick randomization of individual drop trajectories in multidrop simulations). In contrast, for this $Ca$ , our simulations with $\unicode[STIX]{x1D702}$ from 5 $.7\times 10^{-5}$ to 5 $.7\times 10^{-4}$ were simply interrupted at $t\approx 400{-}600$ , providing very reliable time averaging and steady-state rheological functions independent of $\unicode[STIX]{x1D702}$ . For $Ca=0.1$ in figure 10(d), the zero-diffusion solution was limited by $t<15$ for $N_{\triangle }=2160$ , and by $t<22$ for $N_{\triangle }=5120$ , again, due to difficulties with the surface tension reaching zero. Such short times were not enough for accurate time averaging.

For large drop volume fractions and small viscosity ratios $\unicode[STIX]{x1D706}$ , adding even a modest amount of surfactant has the strongest effect on the emulsion rheology. In particular, in figure 10(f) for $c=0.55$ , $\unicode[STIX]{x1D706}=0.25$ and a relatively small elasticity parameter $\unicode[STIX]{x1D6FD}=0.1$ , the excess viscosity $\unicode[STIX]{x1D707}_{pe}-1$ for a contaminated emulsion is about double that for an emulsion of clean drops ( $\unicode[STIX]{x1D6FD}=0$ ), and is a steeper function of the capillary number. A similar strong effect of contamination on $N_{pe}$ is observed. The effects of surfactant are not nearly as prominent for highly viscous drops, which was the main reason to exclude high viscosity ratios $\unicode[STIX]{x1D706}\geqslant O(10)$ from the present study. All the figures 10(a–f) and 11(a,b) point to non-zero (albeit relatively small) values of $N_{pe}$ in the limit $Ca\rightarrow 0$ , and also exclude the possibility of the viscosity function $\unicode[STIX]{x1D707}_{pe}$ approaching its $Ca=0$ limiting value with an error of $O(Ca^{2})$ . These are additional indications that a concentrated emulsion of non-Brownian deformable drops, with or without surfactant, is not a Rivlin–Ericksen fluid in the limit of slow flows $Ca\rightarrow 0$ . Also, the presence of non-diffusing surfactant significantly reduces the critical capillary number for massive cusping or drop break-up to occur (figure 10 f). It was not possible to extend the present simulations beyond this threshold, and it remains unclear to how much extent the other equations of state $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D6E4})$ can relax this limitation.

8.1 Steady planar mixed flow $\unicode[STIX]{x1D712}=0.5$

The first test shows how accurately the Oldroyd model (7.1), with our variable coefficients obtained by fitting (7.1) to the two base flows – PM ( $\unicode[STIX]{x1D712}=0.16$ ) and PE ( $\unicode[STIX]{x1D712}=1$ ) – can predict the PM flow results at intermediate $\unicode[STIX]{x1D712}$ . The constitutive model (7.1) for the mixed-flow kinematics with $\unicode[STIX]{x1D712}=0.5$ was integrated at various capillary numbers (2.7) until the steady state. The initial condition was either $\unicode[STIX]{x1D749}^{d}=0$ or the convergent solution at the nearest available $Ca$ , without an effect on the steady state results. At the transient stage, the argument $\unicode[STIX]{x1D701}$ of the Oldroyd coefficients was allowed to be slightly outside the tabulated range (see figures 12 and 13) with extrapolations used, but it was required to fall within the tabulated range at the steady state. This requirement limited the range of $Ca$ accordingly in the constitutive modelling. Precise results for this flow were obtained by long-time multipole-accelerated BI simulations with $N=100$ drops and sufficient resolution $N_{\triangle }$ (from 2160 to 3840) for convergence with respect to triangulations.

Comparison between the two methods in figure 14(a–c) shows a nearly perfect agreement at various viscosity ratios, emulsion concentrations, surface contaminations and capillary numbers. Unlike in the PM flow test of paper I (figure 23 therein), the present constitutive modelling does not have a substantial lower limitation on the flow intensity, owing to the exclusion of simple shear (with its phase transition difficulties) from the database.

Figure 14.

Steady-state rheological functions for PM flow $\unicode[STIX]{x1D712}=0.5$ with (a) $c=0.55$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.1$ (b) $c=0.55$ , $\unicode[STIX]{x1D706}=0.25$ , $\unicode[STIX]{x1D6FD}=0.1$ and (c) $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.2$ . Lines are from constitutive modelling, and symbols are precise BI results.

8.2 Transient PE flow

Consider a macroscopic PE flow $\boldsymbol{v}(t,\boldsymbol{x})=(\dot{\unicode[STIX]{x1D6E4}}(t)x_{1},-\dot{\unicode[STIX]{x1D6E4}}(t)x_{2},0)$ in dimensional variables with a time-dependent deformation rate $\dot{\unicode[STIX]{x1D6E4}}(t)=\dot{\unicode[STIX]{x1D6E4}}_{max}\,\unicode[STIX]{x1D711}(\dot{\unicode[STIX]{x1D6E4}}_{max}t)$ . The non-dimensional positive function $\unicode[STIX]{x1D711}(s)$ , with $\unicode[STIX]{x1D711}(s)\sim s^{-2}$ at $s\rightarrow -\infty$ , is shown in figure 15. This function was generated numerically to mimic the history of the deformation-rate intensity $(\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640})^{1/2}$ along an open trajectory of a material element approaching a solid sphere from infinity in a steady flow at moderate Reynolds numbers $Re\sim 50$ . For the chosen model flow, the deformation rate $\dot{\unicode[STIX]{x1D6E4}}(t)$ starts from near zero far in the past and reaches a maximum $\dot{\unicode[STIX]{x1D6E4}}_{max}$ at $t=0$ . As pointed out in § 2, implementation of periodic boundaries in long-time multidrop simulations for this type of flow is a simple modification of the Kraynik–Reinelt scheme. A flow was considered with $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , large contamination $\unicode[STIX]{x1D6FD}=0.2$ and $Ca_{max}=2\unicode[STIX]{x1D707}_{e}\dot{\unicode[STIX]{x1D6E4}}_{max}\,a_{o}/\unicode[STIX]{x1D70E}_{eq}=0.2$ ; the latter is the maximum capillary number (2.7) along the trajectory. In the BI simulation used as a precise benchmark, $N=200$ drops with resolution $N_{\triangle }=2160$ started from spherical at $\dot{\unicode[STIX]{x1D6E4}}_{max}t\approx -13$ . Due to mounting difficulties with cusping (see § 4), it was not possible to proceed beyond $\dot{\unicode[STIX]{x1D6E4}}_{max}t\approx 0.014$ ; in the final state (figure 16 a), the drops are strongly deformed. The major rheological characteristic of this type of flow is the viscosity function $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}_{pe}(t)=(\unicode[STIX]{x1D70F}_{11}^{d}-\unicode[STIX]{x1D70F}_{22}^{d})/(4\unicode[STIX]{x1D707}_{e}\dot{\unicode[STIX]{x1D6E4}}_{max})$ . The BI results for this function are shown in figure 16(b) by the solid line. The dotted line in figure 16(b) represents the simplest, quasi-steady approximation $\unicode[STIX]{x0394}\unicode[STIX]{x1D707}_{pe}(t)=\langle \unicode[STIX]{x1D707}_{pe}-1\rangle \unicode[STIX]{x1D711}(\dot{\unicode[STIX]{x1D6E4}}_{max}t)$ , where the steady-state excess non-dimensional viscosity $\langle \unicode[STIX]{x1D707}_{pe}-1\rangle$ in the PE flow, as a function of the instantaneous capillary number $Ca=Ca_{max}\unicode[STIX]{x1D711}(\dot{\unicode[STIX]{x1D6E4}}_{max}t)$ , is interpolated from figure 11(a). This solution had to be terminated at $\dot{\unicode[STIX]{x1D6E4}}_{max}t=-1.53$ to keep $Ca$ in the range of figure 11(a). Even in this range, the quasi-steady approximation is unsuccessful. The dashed line in figure 16(b) is the result of integration of the constitutive model (7.1) with our variable coefficients taken from figure 13(a) and $\unicode[STIX]{x1D749}^{d}=0$ as the initial condition at $\dot{\unicode[STIX]{x1D6E4}}_{max}t\approx -13$ . Here, the argument $\unicode[STIX]{x1D701}$ of the Oldroyd parameters stays within the range of figure 13(a) up to $\dot{\unicode[STIX]{x1D6E4}}_{max}t=-1.25$ , but we continued the solution slightly longer, to $\dot{\unicode[STIX]{x1D6E4}}_{max}t=-0.6$ (with extrapolation from figure 13 a used) to confirm the success of the present constitutive modelling scheme. While the PE quasi-steady approximation is poor, combining the steady-state results for PE with those for a very different steady flow – PM with a small $\unicode[STIX]{x1D712}=0.16$ – in our modelling scheme works to greatly improve the prediction for the time-dependent PE. This improvement appears to rationalize the generalized Oldroyd model (7.1) and the proposed scheme of finding the variable Oldroyd parameters by fitting the model to two distinct steady database flows at arbitrary flow intensities.

Figure 15.

The deformation-rate function $\unicode[STIX]{x1D711}(s)$ used in the transient PE and PM flow tests.

Figure 16.

Transient PE flow test for $c=0.6$ , $\unicode[STIX]{x1D706}=1$ and $\unicode[STIX]{x1D6FD}=0.2$ ; (a) the final snapshot of the BI simulation with 200 drops and $N_{\triangle }=2160$ (the directions of the coordinate axes are the same as in figure 3); (b) the viscosity function from the BI simulation (solid line), constitutive modelling (dashed line) and the quasi-steady approximation (dotted line).

8.3 Transient PM flow ( $\unicode[STIX]{x1D712}=0.16$ )

The next test is for the macroscopic PM flow, written in dimensional variables as $\boldsymbol{v}(t,\boldsymbol{x})=\dot{\unicode[STIX]{x1D6FE}}_{m}(x_{2},\unicode[STIX]{x1D712}x_{1},0)$ . This flow has the same parameter $\unicode[STIX]{x1D712}=0.16$ as one of our steady database flows, but it is now time dependent with the deformation rate $\dot{\unicode[STIX]{x1D6FE}}_{m}(t)=\dot{\unicode[STIX]{x1D6FE}}_{max}\,\unicode[STIX]{x1D711}(\dot{\unicode[STIX]{x1D6FE}}_{max}t)$ ; the function $\unicode[STIX]{x1D711}(s)$ is the same as in figure 15. Implementation of reproducible periodic boundaries in multidrop simulations for this type of flow is a direct generalization of the technique for steady PM flow (§ 2), using an artificial time $\widehat{t}$ defined by $\text{d}\widehat{t}=\dot{\unicode[STIX]{x1D6FE}}_{m}(t)\text{d}t$ in the lattice dynamics part of the algorithm. A flow was chosen with $c=0.55$ , $\unicode[STIX]{x1D706}=0.25$ , moderate drop surface contamination $\unicode[STIX]{x1D6FD}=0.1$ and the maximum capillary number $Ca_{max}=(1+\unicode[STIX]{x1D712})\unicode[STIX]{x1D707}_{e}\dot{\unicode[STIX]{x1D6FE}}_{max}a_{o}/\unicode[STIX]{x1D70E}_{eq}=0.29$ along the trajectory. In the benchmark BI simulation with 200 drops starting from spherical at $\dot{\unicode[STIX]{x1D6FE}}_{max}t\approx -13$ , it was not possible to proceed beyond $\dot{\unicode[STIX]{x1D6FE}}_{max}t=1.34$ due to, again, mounting difficulties with multiple drop cusping. The final state with strong drop deformation is shown in figure 17(a).

Figure 17.

Transient PM flow test for $\unicode[STIX]{x1D712}=0.16$ , $c=0.55$ , $\unicode[STIX]{x1D706}=0.25$ and $\unicode[STIX]{x1D6FD}=0.1$ ; (a) the final snapshot of the BI simulation with 200 drops and $N_{\triangle }=2880$ (the directions of the coordinate axes are the same as in figure 3); (b–d) the rheological functions from the BI simulation (solid lines), constitutive modelling (dashed lines) and the quasi-steady approximation (dotted lines).

The non-dimensional rheological functions for this transient flow are defined as $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{ij}=\unicode[STIX]{x1D70F}_{ij}^{d}/[\unicode[STIX]{x1D707}_{e}\dot{\unicode[STIX]{x1D6FE}}_{max}(1+\unicode[STIX]{x1D712})]$ . The precise BI results for three independent components $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{12}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{22}$ are presented in figure 17(b–d), respectively, by solid lines. The dotted lines in figure 17(b–d) are from the quasi-steady approximation, assuming the stationary flow at the instantaneous capillary number $Ca=Ca_{max}\unicode[STIX]{x1D711}(\dot{\unicode[STIX]{x1D6FE}}_{max}t)$ . The quasi-steady results are obtained by interpolating the PM database rheological functions $\unicode[STIX]{x1D708}_{12}-1$ , $\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x1D708}_{22}$ from figure 7(f) and then multiplying them by $\unicode[STIX]{x1D711}(\dot{\unicode[STIX]{x1D6FE}}_{max}t)$ . The quasi-steady solution was terminated at $\dot{\unicode[STIX]{x1D6FE}}_{max}t=-1.27$ to keep $Ca$ within the range of figure 7(f). The quasi-steady model is again very unsuccessful; in particular, it fails to predict the near-zero values of $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{11}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{22}$ for most of the simulation. The dashed lines in figure 17(b–d) are from integrating the constitutive model (7.1) with our variable Oldroyd coefficients taken from figure 12(c) and the initial condition $\unicode[STIX]{x1D749}^{d}=0$ at $\dot{\unicode[STIX]{x1D6FE}}_{max}t\approx -13$ . The argument $\unicode[STIX]{x1D701}$ of the Oldroyd coefficients was within the range $\unicode[STIX]{x1D701}<0.17$ of figure 12(c) for $\dot{\unicode[STIX]{x1D6FE}}_{max}t<-0.95$ only, but the constitutive modelling (7.1) was extended to $\dot{\unicode[STIX]{x1D6FE}}_{max}t=-0.5$ (with extrapolation from figure 12(c) when $0.17<\unicode[STIX]{x1D701}<0.24$ ) for comparison with the BI solution. Again, while the PM quasi-steady approximation is very poor, a combination of the PM steady-state results with those for steady-state PE in our constitutive modelling dramatically improves the prediction for a time-dependent PM flow. On these grounds, the generalized Oldroyd equation (7.1) with variable coefficients depending on a single instantaneous flow invariant appears to be a meaningful constitutive model for concentrated emulsions.

8.4 Oscillatory PM flow ( $\unicode[STIX]{x1D712}=0.16$ )

The final test is for the oscillating PM flow, written in dimensional variables as $\boldsymbol{v}(t,\boldsymbol{x})=\dot{\unicode[STIX]{x1D6FE}}_{m}(x_{2},\unicode[STIX]{x1D712}x_{1},0)$ . The deformation rate is now defined as $\dot{\unicode[STIX]{x1D6FE}}_{m}(t)=\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle \unicode[STIX]{x1D713}(\unicode[STIX]{x1D714}\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle t)$ with constant $\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle$ and the non-dimensional frequency parameter $\unicode[STIX]{x1D714}$ . The positive, periodic non-dimensional function $\unicode[STIX]{x1D713}(s)$ has been numerically generated to mimic the history of the deformation-rate intensity $(\unicode[STIX]{x1D640}\boldsymbol{ : }\unicode[STIX]{x1D640})^{1/2}$ along a closed trajectory of a material element in the eddy zone behind a sphere in finite Reynolds number flow. This function has unit period and unit average $\langle \unicode[STIX]{x1D713}\rangle$ ; one cycle of $\unicode[STIX]{x1D713}(s)$ is shown in figure 18(a). Thus, the chosen kinematics provides strong, non-harmonic oscillations in the deformation rate $\dot{\unicode[STIX]{x1D6FE}}_{m}(t)$ with the average of $\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle$ and the period of $1/(\unicode[STIX]{x1D714}\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle )$ . A flow was considered with $\unicode[STIX]{x1D714}=1/8$ , $c=0.6$ , $\unicode[STIX]{x1D706}=1$ , $\unicode[STIX]{x1D6FD}=0.2$ and $\langle Ca\rangle =(1+\unicode[STIX]{x1D712})\unicode[STIX]{x1D707}_{e}\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle a_{o}/\unicode[STIX]{x1D70E}_{eq}=0.0355$ . The latter is the average capillary number (2.7) along the trajectory. The flow parameter $\unicode[STIX]{x1D712}=0.16$ is again the same as in our steady database PM flows, to assist in constructing the quasi-steady approximation for comparison with the proposed constitutive modelling and benchmark BI simulation.

The non-dimensional rheological functions for the oscillating PM flow are now defined as $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{ij}^{o}=\unicode[STIX]{x1D70F}_{ij}^{d}/[\unicode[STIX]{x1D707}_{e}\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle (1+\unicode[STIX]{x1D712})]$ . The quasi-steady results for $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{ij}^{o}$ , based on interpolation from figure 8(a) at the instantaneous capillary number $Ca=\langle Ca\rangle \unicode[STIX]{x1D713}(\unicode[STIX]{x1D714}\langle \dot{\unicode[STIX]{x1D6FE}}_{m}\rangle t)$ , are strictly periodic functions of time. The constitutive modelling, based on integrating the Oldroyd equation (7.1) with the variable coefficients from figure 13(a), quickly settles to a periodic regime after just one cycle. In contrast, the multidrop BI simulations do not precisely reach a periodic behaviour. Modest stochastic, statistically steady fluctuations around a periodic regime do not seem to abate even at large times. These fluctuations result from finite system size and can be made small by increasing the number $N$ of drops in a periodic cell; we used $N=400$ for this purpose. The long-time behaviour of $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{ij}^{o}$ over one cycle in the $N=400$ and $N_{\triangle }=2160$ BI simulation is shown in figure 18(b–d) by the solid lines. For the chosen, relatively small frequency parameter $\unicode[STIX]{x1D714}=1/8$ , the quasi-steady approximation (dotted lines in figure 18 b–d) could be expected to work; however, it is found to be reasonable only for $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{12}^{o}$ . The two other components $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{11}^{o}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{22}^{o}$ (not small compared to $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{12}^{o}$ ) are predicted very poorly by the quasi-steady approximation. Again, it is remarkable how the improvements are achieved by the present constitutive modelling scheme (dashed lines in figure 18 b–d) for this time-dependent kinematics, based on the rheological functions for two distinct steady flows. Namely, for $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{12}^{o}$ , the maximum and average deviations from the benchmark BI results are approximately three times smaller. Most significantly, for $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{11}^{o}$ , the proposed constitutive modelling reduced the maximum and average errors seven and 3.5 times, respectively, compared to the quasi-steady approximation. For $\unicode[STIX]{x0394}\unicode[STIX]{x1D708}_{22}^{o}$ , the improvement is less pronounced; still, the maximum and average errors are reduced three and 1.6 times, respectively. We arrived at similar conclusions for a different frequency parameter $\unicode[STIX]{x1D714}=1/6$ .

Figure 18.

Oscillatory PM flow test for $\unicode[STIX]{x1D712}=0.16$ , $c=0.6$ , $\unicode[STIX]{x1D706}=1$ and $\unicode[STIX]{x1D6FD}=0.2$ ; (a) the deformation-rate function $\unicode[STIX]{x1D713}(s)$ ; (b–d) the rheological functions from the $N=400$ BI simulation (solid lines), constitutive modelling (dashed lines) and the quasi-steady approximation (dotted lines).