2. Mathematical formulation and numerical simulation

The dispersed continuous system is governed by the incompressible momentum conservation equations for velocity $\boldsymbol{u}$ in the entire domain $\varOmega$

(2.1) \begin{align}& \frac{\partial (\rho {\textbf u})}{\partial t} + \nabla \cdot (\rho {\textbf u}{\textbf u}) = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} - \int _{\partial \boldsymbol{B}} \textrm{d}{x_{\boldsymbol{B}}}\kappa {\boldsymbol{n}}\Gamma \delta (\boldsymbol{x} - {\boldsymbol{x}}_{\boldsymbol{B}}) ,\end{align}

(2.2) \begin{align}&\qquad\qquad\qquad\qquad\,\, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0 . \end{align}

The total stress $\boldsymbol{\tau }$ has pressure, viscoelastic and viscous parts

(2.3) \begin{equation} \boldsymbol{\tau }=-p\boldsymbol{I}+\boldsymbol{T }^{p}+\boldsymbol{T }^{v},\quad \boldsymbol{T }^{v}={\mu}_{s}\boldsymbol{D} ,\end{equation}

where p is the pressure, ${\mu} _{s}$ is the solvent viscosity and $\boldsymbol{D}=(\boldsymbol{\nabla }\boldsymbol{u})+(\boldsymbol{\nabla }\boldsymbol{u})^{T}$ is twice the deformation rate tensor. The superscript T represents the transpose and $\boldsymbol{T }^{p}$ is the viscoelastic stress due to the presence of polymer. In (2.1), $\rho$ is the density, $\Gamma$ is the surface tension (constant) along the drop surface $\partial B$ consisting of $\boldsymbol{x}_{B}$ points, $\kappa$ is the local curvature, $\boldsymbol{n}$ is the outward normal and $\delta (\boldsymbol{x}-\boldsymbol{x}_{B})$ is the three-dimensional Dirac delta function. We use a modified Chilcott–Rallison-type (Chilcott & Rallison Reference Chilcott and Rallison1988; Matos et al. Reference Matos, Alves and Oliveira2009) constitutive equation (also called the finite extensible nonlinear elastic modified Chilcott–Rallison or FENE-MCR) to model the viscoelasticity of the continuous phase. In our investigations of the effects of viscoelasticity, we have always chosen simple constitutive equations to explain the underlying physics. Unlike the Oldroyd-B equation used in our earlier studies (Aggarwal & Sarkar, Reference Aggarwal and Sarkar2007, Reference Aggarwal and Sarkar2008a , Reference Aggarwal and Sarkarb ; Mukherjee & Sarkar Reference Mukherjee and Sarkar2009, Reference Mukherjee and Sarkar2010), the FENE-MCR model in our recent studies (Mukherjee & Sarkar Reference Mukherjee and Sarkar2011, Reference Mukherjee and Sarkar2014) has a finite extensible viscosity; like Oldroyd-B, it has a constant shear viscosity. It has been extensively used in modelling different viscoelastic flows (Szabo et al. Reference Szabo, Rallison and Hinch1997; Ramaswamy & Leal Reference Ramaswamy and Leal1999; Dou & Phan-Thien Reference Dou and Phan-Thien2003; Kim et al. Reference Kim, Kim, Chung, Ahn and Lee2005). Furthermore, this constant viscosity viscoelastic model is applicable to a Boger fluid, which is also routinely used in experiments to investigate viscoelastic behaviours. e.g. for pair interactions between spheres in Boger fluids (Snijkers et al. Reference Snijkers, Pasquino and Vermant2013). Details of the implementation can be found in our previous work (Mukherjee et al. Reference Mukherjee, Tarafder, Malipeddi and Sarkar2022; Tarafder et al. Reference Tarafder, Malipeddi and Sarkar2022). The viscoelastic stress is determined by a simple rate type of equation,

(2.4) \begin{equation} \frac{\partial \boldsymbol{T}^{p}}{\partial t}+\left\{\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol{T}}^{p}-\boldsymbol{\nabla }\boldsymbol{u}\boldsymbol{\cdot }{\boldsymbol{T}}^{p}-{\boldsymbol{T}}^{p}\boldsymbol{\cdot }\boldsymbol{\nabla }{\boldsymbol{u}}^{T}\right\}+\frac{f}{\lambda }{\boldsymbol{T}}^{p}=\frac{f}{\lambda }{\mu} _{p}\boldsymbol{D}, f=\frac{L^{2}+\dfrac{\lambda }{{ \mu} _{p}}\left(\sum T_{ii}^{p}\right)}{L^{2}-3} .\end{equation}

Here, ${ \mu} _{p}$ is the polymeric viscosity, $\lambda$ the relaxation time and $L$ the finite extensibility, introduced by the FENE-CR model, which limits the maximum length of the end-to-end vector for the polymer molecule. In the limit of $L\rightarrow \infty$ we obtain the Oldroyd-B equation with $f\rightarrow 1$ in (2.4). We use $L=20$ , increasing which was shown to make little difference to the results (Mukherjee & Sarkar Reference Mukherjee and Sarkar2013). We have provided additional details about the modified FENE-MCR model, noting that the nonlinearity does not affect the linear viscoelastic response, and results in a finite extensional viscosity and a constant shear viscosity (Oliveira Reference Oliveira2003; Tarafder et al. Reference Tarafder, Malipeddi and Sarkar2022). In spite of its simplicity, we feel that the model is appropriate for the present purpose of explaining a phenomenon. However, note that shear thinning may qualitatively influence the problem of the particle dynamics in a viscoelastic medium. Studies of shear-induced string-like structure formation in a viscoelastic particulate flow have found no such structure in non-shear-thinning Boger fluids (Scirocco et al. Reference Scirocco, Vermant and Mewis2004, Won & Kim Reference Won and Kim2004). Equations (2.1) and (2.2) along with the viscoelastic constitutive equation (2.4) are solved by a semi-implicit finite difference projection method in a Cartesian domain. An alternating direction implicit scheme is applied to ease the restriction on the time step. A multigrid method is used to solve the pressure position equation. A semi-analytic time integration method is used, which automatically achieves an elastic viscous stress splitting, alleviating some of the numerical stability issues common to viscoelastic numerical schemes (Sarkar & Schowalter Reference Sarkar and Schowalter2000; Izbassarov & Muradoglu Reference Izbassarov and Muradoglu2015). Details of the implementation of the above algorithm can be found in Li & Sarkar (Reference Li and Sarkar2005), Aggarwal & Sarkar (Reference Aggarwal and Sarkar2007) and Mukherjee & Sarkar (Reference Mukherjee and Sarkar2013).

The simulation details for computing the shear-induced collective diffusivity are identical to our previous publications on viscous systems (Malipeddi & Sarkar Reference Malipeddi and Sarkar2019a , Reference Malipeddi and Sarkarb , Reference Malipeddi and Sarkar2021). A total of N = 70 drops with radius a are placed in a random configuration in an initially compact layer at the middle (width ∼ 0.2 $L_{y}$ ) of the computational domain (figure 1) of size $L_{y}=28a$ in the velocity gradient direction (y) and $L_{x}=L_{z}=14a$ along the flow (x) and the vorticity (z) directions, leading to an initial volume fraction of ∼25 % in the compact layer. The independence of the result on domain dimensions, drop numbers and initial compact layer volume fractions (in the range of ∼25 %–43 %) has been investigated in a previous study of a viscous system (Malipeddi & Sarkar Reference Malipeddi and Sarkar2019b ), indicating that the current configuration is sufficient for the simulation. Specifically, the domain length of 28a along the velocity gradient direction has been shown to be sufficient to neglect any wall effects (results differed negligibly from those obtained with a domain length of 42a). Note that this approach differs from previous studies of hydrodynamic diffusivities (King & Leighton Reference King and Leighton2001; Hudson Reference Hudson2003), where drops uniformly distributed between walls experience competing effects of shear-induced diffusion and wall-induced migration. The computational domain is discretised with $96\times 192\times 96$ grid points (∼14 grid points per drop diameter), which was shown to be sufficient in our earlier studies (Srivastava et al. Reference Srivastava, Malipeddi and Sarkar2016). The top and the bottom walls move with equal and opposite velocity U to obtain a shear rate $\dot{\gamma }=2U/L_{y}$ . Initially, the drops are placed in a fully developed shear flow (appropriate for Stokes flow) with no polymeric stresses. Periodic boundary conditions have been applied in the x and z-directions. The drop layer has a homogeneous distribution along the x and z directions. The drop radius a and the inverse shear rate $\dot{\gamma }^{-1}$ are used as the length and the time scales to define the Reynolds number ${Re}=\rho _{m}\dot{\gamma }a^{2}/{ \mu} _{m}$ , capillary number $Ca={ \mu} _{m}\dot{\gamma }a/{\unicode[Arial]{x0393}}$ and Weissenberg number $Wi=\lambda \dot{\gamma }$ . The viscosity ratio $\lambda _{{ \mu} }={ \mu} _{d}/{ \mu} _{m}$ and the density ratio $\lambda _{\rho }=\rho _{d}/\rho _{m}$ have been kept at unity. The polymeric to total matrix viscosity ratio $\beta ={ \mu} _{pm}/{ \mu} _{m}$ is fixed at 0.5 (in Appendix, we briefly investigated the effects of $\beta$ variation). Subscripts m and d denote fluid and drop phases, respectively. The total viscosity of the surrounding fluid ${ \mu} _{m}$ comprises the polymeric and the solvent viscosities ${ \mu} _{m}={ \mu} _{sm}+{ \mu} _{pm}$ . Because of the explicit nature of the code and thereby the diffusion limitation on time stepping, Re has been kept to 0.1 as a proxy for Stokes flow. We have previously shown it to be sufficient for matching with boundary element simulation of Stokes flow of viscous emulsions (Srivastava et al. Reference Srivastava, Malipeddi and Sarkar2016). Note that the front-tracking implementation is a smoothed-interface method, with the drop’s surface moving with the local velocity interpolated from the Eulerian grid to the Lagrangian drop front, naturally avoiding interpenetration; no physical drop coalescence is considered. The above code was run with the help of George Washington University’s High Performance Computing cluster PEGASUS.

Figure 1. Schematic of the computational set-up.

The shear-induced collective diffusion of the drops in an unbounded shear (ensured by a large enough computational domain; see above) leads to a spreading of the initial layer of drops concentrated in the y-direction (homogeneous in the other two directions). It is governed by a coarse-grained one-dimensional diffusion equation for the drop volume fraction $\phi (y,t)$ (Hudson Reference Hudson2003)

(2.5) \begin{equation} \frac{\partial \phi }{\partial t'}=\frac{\partial }{\partial y'}\left(\frac{D_{yy}^{c}\left(\phi \right)}{\dot{\gamma }a^{2}}\frac{\partial \phi }{\partial y'}\right) ,\end{equation}

where $D_{yy}^{c}=\dot{\gamma }\phi a^{2}f_{2}$ is the collective diffusivity coefficient, which is a function of local volume fraction $\phi (y,t)$ and $f_{2}$ is the non-dimensional collective diffusivity. The imposed unbounded shear here eliminates the advection term due to drop migration found in Hudson (Reference Hudson2003). As has been discussed in detail in our previous article (Malipeddi & Sarkar, Reference Malipeddi and Sarkar2019b ), the linear dependence of $D_{yy}^{c}$ on the volume fraction $\phi$ is predicated on the dominance of the pairwise interactions validated a posteriori by the simulation results (also see Malipeddi & Sarkar, Reference Malipeddi and Sarkar2019a , Reference Malipeddi and Sarkar2021). We have shown that $f_{2}$ does not change with changing initial compact layer volume fraction (in the range of 25 %–43 %). Equation (2.5) for a fixed number of particles spreading by diffusion yields a self-similar parabolic volume fraction (Grandchamp et al. Reference Grandchamp, Coupier, Srivastav, Minetti and Podgorski2013)

(2.6) \begin{equation} \psi \left(\eta \right)=\left(f_{2}t'\right)^{1/3}\phi =\left(b-\eta ^{2}/6\right),\begin{array}{ll} & \end{array}\eta =y'\!/\!\left(f_{2}t'\right)^{1/3}, \end{equation}

with b as a free parameter, showing a 1/3rd scaling of the half-width $\underline{w}$ of the profile with time

(2.7) \begin{equation} \underline{w}^{3}-\underline{w}_{o}^{3}=K\hspace{0pt}t',\begin{array}{ll} & \end{array}K=9f_{2}N_{0}/\left(4\sqrt{2}\right), \begin{array}{ll} & \end{array} N_{0}=\int \phi \left(y',t'\right){\textrm d}y', \end{equation}

allowing determination of $f_{2}$ . Here, $\underline{w}_{o}$ is the initial width and $N_{0}$ is a conserved quantity representing the number of drops diffusing out from a layer, calculated as $N_{0}=NV/aL_{x}L_{y}$ . Malipeddi & Sarkar (Reference Malipeddi and Sarkar2019b ) have identified an alternative half-width $w$ computed from the drop positions $y_{i}^{\prime}, i=1,\ldots, N$

(2.8) \begin{equation} w=\sqrt{\frac{1}{N}\sum _{i=1}^{N}\left(y_{i}^{\prime}-{ \mu} \right)^{2}},\quad { \mu} =\frac{1}{N}\sum _{i=1}^{N}y_{i}^{\prime}. \end{equation}

It can be related to the first moment of the analytical solution (2.6) leading to

(2.9) \begin{equation} w^{3}-w_{0}^{3}=K't',\qquad K'=9f_{2}N_{0}/\big(10\sqrt{5}\big) .\end{equation}

Using either expression for finding $f_{2}$ led to identical results within estimation uncertainties (Malipeddi & Sarkar, Reference Malipeddi and Sarkar2019b ). However, the second method avoids the intermediate step of obtaining the coarse-grained volume fraction from the drop location. We plot $w^{3}-w_{0}^{3}$ vs. $t^{\prime}(=t\dot{\gamma })$ and calculate the slope to find the collective diffusivity. The simulations are run up to $t^{\prime}$ ∼ 180 inverse shear unit times and $t^{\prime}\leq t_{0}^{\prime}=20$ is discarded to make sure the drops are deformed and a linear region of $w^{3}-w_{0}^{3}$ vs. $t^{\prime}$ is reached. Here, $w_{0}$ represents drop layer width at $t^{\prime}=t_{0}^{\prime}$ . Similar to Malipeddi & Sarkar (Reference Malipeddi and Sarkar2019b ), the time evolution data excluding the initial transient are divided into smaller overlapping time intervals (∼80 non-dimensional times) with $f_{2}$ computed in each of them; their average is reported as an ensemble average with corresponding standard deviation as the uncertainty in its estimation. Previously, we have shown that $f_{2}$ does not depend on the initial volume fraction (in the range of ∼25 %–43 %) of the drop layer (Malipeddi & Sarkar, Reference Malipeddi and Sarkar2019b ).

Shear-induced particle diffusivity, both self and gradient (or collective), can also be obtained from the dynamic structure factor (Rallison & Hinch Reference Rallison and Hinch1986; Morris & Brady Reference Morris and Brady1996; Marchioro & Acrivos Reference Marchioro and Acrivos2001; Leshansky & Brady Reference Leshansky and Brady2005, Leshansky et al. Reference Leshansky, Morris and Brady2008). The dynamic structure factor has been useful for particle sizing in dynamic light scattering techniques, where the scattering of a monochromatic laser from the suspension volume is analysed to obtain the fluctuation autocorrelation. The exponential decay time of the autocorrelation is related to the diffusivity, which in turn is related to the particle size. Measurement of gradient diffusivity using the dynamic structure factor has previously been performed in homogeneous suspensions where the decay of spontaneously appearing fluctuation underlies the wavenumber or scale-dependent diffusivities. However, we have convincingly demonstrated that it can be used in a non-homogeneous system such as the centrally packed layer considered here to get useful information. It offered an excellent match with the diffusivity obtained by the continuum diffusion equation method (Malipeddi & Sarkar, Reference Malipeddi and Sarkar2019a , Reference Malipeddi and Sarkarb , Reference Malipeddi and Sarkar2021). A detailed description of the method, including a discussion of the various issues arising from inhomogeneity of the system, has been given in Malipeddi and Sarkar (Malipeddi & Sarkar, Reference Malipeddi and Sarkar2019a ). Briefly, the method is based on the particle autocorrelation function in the Fourier ( $\boldsymbol{k}$ wavenumber) space of $N$ drops located at $\boldsymbol{x}_{\alpha }^{\prime}(t'),\alpha =1,2,\ldots, N$

(2.10) \begin{equation} \begin{array}{l} F\big(\boldsymbol{k},t'\big)=\dfrac{1}{N}\big\langle \sum _{\alpha ,\beta =1}^{N}e^{i\boldsymbol{k}\boldsymbol{\cdot }({\boldsymbol{x}_{\alpha }^{\prime}}(t')-{\boldsymbol{x}'_{\beta }}(0))}\big\rangle ,\\[8pt] n\left(\boldsymbol{x}',t'\right)=\sum _{\alpha =1}^{N}\delta \left(\boldsymbol{x}'-\boldsymbol{x}_{\alpha }^{\prime}\right),\begin{array}{ll} & \end{array}\hat{n}\left(\boldsymbol{k},t'\right)=\sum _{\alpha =1}^{N}e^{i\boldsymbol{k}\boldsymbol{\cdot }{\boldsymbol{x}_{\alpha }^{\prime}}}. \end{array} \end{equation}

It can be shown that the number density $n(\boldsymbol{x}',t')$ satisfies an advection-diffusion equation (Leshansky & Brady Reference Leshansky and Brady2005)

(2.11) \begin{equation} \frac{\partial n}{\partial t}+\left(\boldsymbol{U}+\dot{\boldsymbol{\varGamma }}\boldsymbol{\cdot }\boldsymbol{x}\right)\boldsymbol{\cdot }\boldsymbol{\nabla }n=\boldsymbol{D}^{c}\boldsymbol{\nabla} ^{2}n \end{equation}

in a shear flow $\boldsymbol{U}+\dot{\boldsymbol{\varGamma }}\boldsymbol{\cdot }\boldsymbol{x}$ ( $\boldsymbol{U}$ is the average flow and $\dot{\boldsymbol{\varGamma }}$ is the velocity gradient tensor), with the diffusivity in the gradient direction being related to $F(\boldsymbol{k},t')$ as

(2.12) \begin{equation} D_{yy}^{c}=-\frac{1}{k^{2}}\frac{{\textrm d}\left(\ln F\right)}{{\textrm d}t'} .\end{equation}

Note that, as noted in our previous publication (Malipeddi & Sarkar Reference Malipeddi and Sarkar2019b ). the advection term in (2.11) does not affect the dynamics in a simple shear due to the orthogonality of the $\boldsymbol{k}$ ( $=k\hat{\boldsymbol{y}}$ ) vector to the velocity field. We showed before that, even in this initially non-homogeneous system, the method can be applied to obtain a meaningful diffusivity $D_{yy}^{c}$ in the limit of $k\rightarrow 0$ which can be matched with the results obtained from the continuum method described before. Note that, similar to the two different methods to compute $f_{2}$ described before, the autocorrelation can also be computed either directly from the particle positions (2.10) or using the coarse-grained volume fraction $\phi (y,t)$ satisfying the continuum (2.5). Both were shown to yield the same $D_{yy}^{c}$ (Malipeddi & Sarkar Reference Malipeddi and Sarkar2019b ). To compute $F(\boldsymbol{k},t')$ , signals for $t'\lt t_{0}$ are discarded and the rest are used to calculate the autocorrelation in 80 % overlapping intervals. They are averaged with the corresponding standard deviation used as an error of estimation. We will compare the gradient diffusivity represented by $f_{2}$ and $D_{yy}^{c}$ .