2 Model description

Numerous numerical techniques have been used to study particle suspensions. Among them we find dissipative particle dynamics (Hoogerbrugge & Koelman Reference Hoogerbrugge and Koelman1992), the lattice Boltzmann method (Ladd & Verberg Reference Ladd and Verberg2001) as well as the Lagrange multiplier fictitious domain method (Glowinski et al. Reference Glowinski, Pan, Hesla and Joseph1999) and the finite element method (Maury Reference Maury1999). The most widely used technique for suspensions at zero Reynolds number is Stokesian dynamics (Brady & Bossis Reference Brady and Bossis1988; Sierou & Brady Reference Sierou and Brady2001). In this paper, we use a two-dimensional finite element method for rigid non-Brownian fully lubricated particles in Stokes flow (i.e. zero Reynolds number and infinite Péclet number). The implementation is based on MILAMIN (Dabrowski, Krotkiewski & Schmid Reference Dabrowski, Krotkiewski and Schmid2008), which is an efficient open-source implementation of a finite element model Stokes solver in MATLAB; we have used this model previously to investigate the transient nature of suspensions and the statistics of close particle encounters (Thøgersen, Dabrowski & Malthe-Sørenssen Reference Thøgersen, Dabrowski and Malthe-Sørenssen2016).

We solve the incompressible Stokes equations. Conservation of mass yields

where $\boldsymbol{v}$ is the velocity field. Conservation of momentum yields

where $\mathbf{0}$ ensures a neutrally buoyant suspension. Here $\unicode[STIX]{x1D748}$ is the stress tensor,

where $p$ is the pressure, $\unicode[STIX]{x1D644}$ is the identity matrix and

is the deviatoric stress tensor.

We introduce an approximation space through a set of basis functions for velocity and pressure, which are usually called shape functions. We use a mixed formulation finite element method and the seven-node Crouzeix–Raviart triangular element (Crouzeix & Raviart Reference Crouzeix and Raviart1973). We use $N_{i}$ to denote continuous shape functions for velocity, and $\unicode[STIX]{x1D6F1}_{i}$ for linear discontinuous shape functions for pressure. Approximations of velocities and pressure marked with $\tilde{~}$ are then

and

where $nnod$ and $np$ are the number of velocity and pressure degrees of freedom, respectively. In Voigt notation, we define the following quantities at the element level:

The weak formulation of Stokes equations can then be written through the following matrices:

and

where the integrals over the domain $\unicode[STIX]{x1D6FA}$ are carried out at an element level. Using isoparametric elements, the elements are mapped onto a fixed reference frame in the local coordinates $[\unicode[STIX]{x1D712},\unicode[STIX]{x1D701}]$ . The derivatives of the shape functions of global coordinates $x$ and $y$ are found from the derivatives of the shape functions with respect to the local coordinates:

where $\unicode[STIX]{x1D645}$ is the Jacobian,

The volume integrals over $\unicode[STIX]{x1D6FA}$ are carried out using elementwise quadrature, which yields sums over $nip$ integration points located at $[\unicode[STIX]{x1D712}_{k},\unicode[STIX]{x1D701}_{k}]$ :

where $|\unicode[STIX]{x1D645}|$ is the Jacobian determinant, and $W_{k}$ are integration point specific weights. Using a weighted residual Galerkin method yields a symmetric, albeit indefinite, matrix equation. The full matrix assembly yields the global matrix equation, where we use the augmented Lagrangian method and solve

Here, we have introduced the general term $[\boldsymbol{g},\boldsymbol{h}]$ for the right-hand side after boundary conditions have been applied. In practice,

is solved iteratively for $p$ using a preconditioned conjugate gradient method, with $\unicode[STIX]{x1D648}$ as the preconditioner. Since the pressures shape functions are discontinuous, $\unicode[STIX]{x1D648}^{-1}$ can be obtained at the element level, which improves performance significantly. For the augmented matrix $\unicode[STIX]{x1D63C}+\unicode[STIX]{x1D706}\unicode[STIX]{x1D64C}\unicode[STIX]{x1D648}^{-1}\unicode[STIX]{x1D64C}^{\text{T}}$ , we use Cholesky decomposition, which only has to be carried out once, so that the incremental steps can be calculated efficiently using the method of backward–forward substitution. The resulting incompressible velocity field can be recovered from

For a more detailed discussion about implementation and optimization procedures, we refer to Dabrowski et al. (Reference Dabrowski, Krotkiewski and Schmid2008).

The rigidity of particles is enforced directly on the matrix level through a matrix transformation. Using an operator $\unicode[STIX]{x1D64B}$ , we replace $\unicode[STIX]{x1D63C}$ , $\unicode[STIX]{x1D64C}$ and $\boldsymbol{g}$ in (2.20) with

and

respectively. The resulting matrix equation is still symmetric and positive definite, and ensures that the net forces and torques on the particles vanish. The particle motion is fully coupled to the velocity field of the fluid phase, as we solve for all velocities simultaneously. Note that there are no constraints on the particle shape, as long as the particle centre of mass is given. The operator $\unicode[STIX]{x1D64B}$ can be found in appendix A. Validation of the implementation can be found in appendix B.

3 Discretization and system set-up

Direct modelling of particle suspensions is numerically challenging. In particular, in finite element models with unstructured body-fitting meshes, particle overlaps is a problem if we wish to run systems to very large strains. Increasing the resolution between particles is crucial for accuracy, but, since particles can be arbitrarily close, this may lead to unconstrained computational costs, and ill-conditioned systems of equations due to numerous small elements in between close particle pairs. Still, mesh refinement in particle apertures as well as a small enough time step make it possible to run the system to large enough strains to address shear-induced self-diffusion for particle volume fractions up to $\unicode[STIX]{x1D719}=0.4$ , without introducing any artificial repulsive forces between the particles.

In this paper we study the purely hydrodynamic limit without repulsive forces between the particles. In order to achieve this, special care has to be taken over spatial discretization as well as over time integration. We generate an unstructured triangular mesh at every time step using ‘Triangle’ (Shewchuk Reference Shewchuk1996), and use adaptive mesh refinement to ensure that there are at least two mesh elements between all particle pairs at any time during the simulation. In this way we ensure that the velocity field is accurately computed between close particles (quadratic velocity profiles can be represented within a single triangular element). For time integration we use the second-order Runge–Kutta method with time step $\dot{\unicode[STIX]{x1D6FE}}\,\text{d}t=0.02$ . If particle overlap is detected, the simulation is terminated, which for dense suspensions sets limits on the maximum strains that we can reach. For statistical measures, it would be beneficial to have one very long run rather than ensemble averaging over multiple shorter ones, but due to computational challenges the latter is used throughout the paper.

A strength of the finite element method is that it allows access to the velocity field in the fluid at all time steps. This means that, in addition to measuring the shear-induced self-diffusion of the solid phase, we can also measure the shear-induced self-diffusivity of the fluid phase using the fluid velocity field. While we acknowledge that particle shape can play a role in migration of particles (Rusconi & Stone Reference Rusconi and Stone2008), we limit this study to monodisperse discs.

The system set-up is demonstrated in figure 1. We use a simulation box with dimensions $L\times L$ . The left and right boundaries are periodic (the rows and columns of the corresponding periodic degrees of freedom are replaced by their sums in the system of equations (2.20)), and we set up a shear rate $\dot{\unicode[STIX]{x1D6FE}}$ in the $x$ -direction by applying Dirichlet boundary conditions at the top and bottom boundaries: constant $v_{x}$ and zero $v_{y}$ . The strain $\unicode[STIX]{x1D6FE}=\dot{\unicode[STIX]{x1D6FE}}\,\text{d}t$ is then the dimensionless time in the system. We initialize the system of particles using random sequential adsorption (Widom Reference Widom1966), with a small initial particle radius. We then increase the radius of the particles and iterate until no overlaps are present using a steepest descent algorithm.

Figure 2 shows a snapshot of a simulation of $N=4096$ particles at volume fraction $\unicode[STIX]{x1D719}=0.3$ in simple shear. The presence of particles perturbs the velocity field and sets up a non-zero velocity field in the $y$ -direction, both for the particles and for the fluid. In figure 2(c) one can see a low-velocity layer ( $y$ -component) close to the boundary that is approximately two particle diameters thick. In bounded domains, one can also find swapping trajectories (Zurita-Gotor, Bławzdziewicz & Wajnryb Reference Zurita-Gotor, Bławzdziewicz and Wajnryb2007) as well as layered hexagonal structures close to the boundaries (Gallier et al. Reference Gallier, Lemaire, Lobry and Peters2016). Layers of low velocity close to the boundary introduce boundary effects on measurements of mean square displacement as well as shear-induced self-diffusion that scale with $1/\sqrt{N}$ . Throughout this paper we will use the scaling with $1/\sqrt{N}$ when we look at measurements towards infinite system sizes. An alternative approach that was taken by Gallier et al. (Reference Gallier, Lemaire, Peters and Lobry2014) is to perform measurements in the core region (mid 50 %) in order to remove boundary effects. Figure 3 shows the probability density of the fluid velocity field for systems of $N=1024$ and $\unicode[STIX]{x1D719}\in [0.1,0.4]$ equilibrated for $\dot{\unicode[STIX]{x1D6FE}}t=20$ from simulations of $4\times 10^{4}$ fluid tracer particles integrated forward in time with a second-order Adams–Bashforth method with time step $\dot{\unicode[STIX]{x1D6FE}}\,\text{d}t=0.02$ . The boundary effects have been removed by removing tracers in a band of width $4r$ near the top and bottom boundaries.

Figure 1.

Sketch of the system set-up. The left and right boundaries are periodic, while the top and bottom boundaries have velocities in opposite directions to give a background shear rate $\dot{\unicode[STIX]{x1D6FE}}$ . The particle configuration is a snapshot from a simulation of $N=1024$ particles at volume fraction $\unicode[STIX]{x1D719}=0.4$ at $\unicode[STIX]{x1D6FE}=30$ .

Figure 2.

Snapshot of simulation at $\unicode[STIX]{x1D6FE}=30$ of 4096 particles at volume fraction $\unicode[STIX]{x1D719}=0.3$ in simple shear. The top and bottom boundaries have fixed velocities in the $x$ -direction, while the left and right boundaries are periodic. (a) The $x$ -component of the velocity field. (b) Close-up of the $x$ -component of the velocity field with subtracted background shear rate. (c) Close-up of the $y$ -component of the velocity field. The presence of particles sets up a non-zero $y$ -component of the velocity field, leading to shear-induced self-diffusion of particles as well as of the surrounding fluid. (d) Close-up of the curl field $\boldsymbol{e}_{z}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{v}$ , where $\boldsymbol{e}_{z}$ is the unit vector out of the plane.

Figure 3.

Probability density of the fluid velocity field for $N=1024$ and $\unicode[STIX]{x1D719}\in [0.1,0.4]$ . (a) The $y$ -component of the velocity field scaled with the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ and the particle radius $r$ . (b) The $x$ -component of the velocity field scaled with the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ and the system size $L$ . (c) The $x$ -component of the velocity, where the background velocity field $\dot{\unicode[STIX]{x1D6FE}}y$ has been subtracted, scaled with the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ and the particle radius $r$ . (The breaking of symmetry in panel (c) is a system size effect, which becomes evident when the particle concentration is not perfectly constant in the $y$ -direction. Since shear-induced self-diffusion is slow, the average shear rate is not always identical to $\dot{\unicode[STIX]{x1D6FE}}y$ for finite strains.)

Figure 4.

Fluid mixing in a system with $N=512$ , $\unicode[STIX]{x1D719}=0.3$ at $Pe=10^{4}$ . The concentration profile is estimated using the diffusive strip method. The dimensions of the window are $(0.2\times 0.2)L$ in the central part of the simulation box.

Figure 5.

Fluid tracer concentration field calculated using the diffusive strip method. The concentration profiles are snapshots of a simulation at $\unicode[STIX]{x1D6FE}=10$ for $N=512$ , $\unicode[STIX]{x1D719}=0.3$ for $Pe=10^{3}$ , $10^{4}$ , $10^{5}$ and $10^{6}$ . The dimensions of the window are $(0.6\times 0.2)L$ in the central part of the simulation box.

Figure 6.

Fluid tracer concentration field calculated using the diffusive strip method. The concentration profiles are snapshots of a simulation at $\unicode[STIX]{x1D6FE}=10$ for $N=512$ , $Pe=10^{4}$ for $\unicode[STIX]{x1D719}=0.1$ , $0.2$ , $0.3$ and $0.4$ . The dimensions of the window are $(0.6\times 0.2)L$ in the central part of the simulation box.

4 Structure of mixing in the fluid phase

While we assume that the particle Péclet number is infinite, we assume that diffusion of mass is present in the fluid phase. Recent experiments by Souzy et al. (Reference Souzy, Yin, Villermaux, Abid and Metzger2015) have shown complex structures in the fluid phase when a strip of dye is placed near the boundary in a sheared suspension. Here, we introduce finite fluid tracer Péclet numbers $Pe=\dot{\unicode[STIX]{x1D6FE}}r^{2}/D$ , where $D$ is the fluid tracer diffusivity, $r$ is the particle radius and $\dot{\unicode[STIX]{x1D6FE}}$ is the shear rate. We take advantage of the velocity fields that we calculate using finite elements, and use the diffusive strip method developed by Meunier & Villermaux (Reference Meunier and Villermaux2010), which allows us to study filaments of fluid tracer particles at finite Péclet numbers with strain.

5 Mixing at small strains – concentration distributions from the diffusive strip method

For the strains reached using the diffusive strip method, the shear-induced self-diffusive regime has not been established (see e.g. the velocity autocorrelation function in figure 11). In order to quantify mixing at these strains, we measure the distribution of tracer concentration fields from the previous section as a function of $\dot{\unicode[STIX]{x1D6FE}}t$ for various $\unicode[STIX]{x1D719}$ and $Pe$ .

Figure 9 shows the concentration distribution $c_{0}/c$ in the fluid phase for $\dot{\unicode[STIX]{x1D6FE}}t\in [0.02,5]$ , $\unicode[STIX]{x1D719}\in [0.1,0.4]$ and $Pe\in [10^{3},10^{6}]$ . Figure 9(a) shows the concentration distribution for $\unicode[STIX]{x1D719}=0.3$ at $Pe=10^{4}$ as $\dot{\unicode[STIX]{x1D6FE}}t$ increases. The initial distribution is U-shaped, with a peak at $c/c_{0}=1$ (as was also described by Meunier & Villermaux (Reference Meunier and Villermaux2010)), which quickly shifts to lower concentrations, and the distribution becomes dominated by the low concentrations between the strips. To highlight instead the changes close to the strip, we also give the probability density of the inverse of the concentration $c_{0}/c-1$ (figure 9 b). Then one can clearly see the decay and shift of the peak concentration as a function of $\dot{\unicode[STIX]{x1D6FE}}t$ . Figure 9(c,d) shows the dependence on the concentration probability density at $\dot{\unicode[STIX]{x1D6FE}}t$ on $\unicode[STIX]{x1D719}$ and $Pe$ , respectively. Increasing the particle volume fraction shifts the peak concentration and widens the inverse concentration distribution. The same applies on decreasing the Péclet number, which demonstrates a strong dependence on the concentration distribution of fluid tracers even at large $Pe$ . This is also clearly visible in figure 5. Interestingly, increasing strain, increasing particle volume fraction and decreasing Péclet number all result in very similar changes to the inverse fluid tracer concentration distribution.

We also measure the average concentration averaged over the $x$ -direction $\langle c/c_{0}\rangle _{x}$ as a function of $\dot{\unicode[STIX]{x1D6FE}}t$ . The shape of the distribution for $Pe=10^{4}$ , $\unicode[STIX]{x1D719}=0.3$ and $N=512$ is shown in figure 10(a). This distribution widens as the strain increases. In order to investigate the effect of $\unicode[STIX]{x1D719}$ and $Pe$ , we measure the width of the average concentration as a function of $\dot{\unicode[STIX]{x1D6FE}}t$ , which we define as the standard deviation of the normalized concentration distribution,

where $\unicode[STIX]{x1D6E5}_{y}$ increases with $\dot{\unicode[STIX]{x1D6FE}}t$ and with increasing $\unicode[STIX]{x1D719}$ . However, it is more or less independent of $Pe$ for $Pe>10^{3}$ . This demonstrates that, although the probability density of concentration has a clear dependence on $Pe$ for $Pe\in [10^{3},10^{6}]$ , the spatial distribution of concentration measured through $\unicode[STIX]{x0394}_{y}$ is largely dominated by advection for $Pe\geqslant 10^{3}$ .

Figure 10.

(a) Concentration averaged over $x$ as a function of $y/r$ averaged over nine diffusive strips shifted to $y/r=0$ at $\dot{\unicode[STIX]{x1D6FE}}t=0$ , for $Pe=10^{4}$ , $\unicode[STIX]{x1D719}=0.3$ and $N=512$ . (b) Width of the average concentration $\unicode[STIX]{x1D6E5}_{y}$ as a function of $\dot{\unicode[STIX]{x1D6FE}}t$ . The width is taken to be the standard deviation of the normalized average concentration distribution.

6 Mixing at large strains – shear-induced self-diffusion

At large strains, the mixing of particles in a suspension close to equilibrium is described by the self-diffusion coefficient. In this section, we measure the self-diffusion coefficient for both the fluid and the solid phase. Before we measure the diffusion coefficient, we equilibrate the system for $\dot{\unicode[STIX]{x1D6FE}}t=20$ , which is chosen because of the decay of the velocity autocorrelation function in the $y$ -direction (figure 11 b). This resulting mixing of the particles can be described using the self-diffusion coefficient $D_{s}$ . Note that the self-diffusion coefficient is actually a tensor of rank 2 in two-dimensional systems that are not isotropic (see e.g. Breedveld et al. Reference Breedveld, Van Den Ende, Bosscher, Jongschaap and Mellema2002). In this work we focus on the translation of particles normal to the direction of shear, and in the following, $D_{s}$ will denote this component of the self-diffusion tensor (often denoted $D_{s,\langle yy\rangle }$ in the literature). Also note that self-diffusion and collective diffusion are in general two separate quantities. In this paper we discuss self-diffusion normal to the direction of shear only. The interested reader is referred to Leighton & Acrivos (Reference Leighton and Acrivos1987b ), Phillips et al. (Reference Phillips, Armstrong, Brown, Graham and Abbott1992), Marchioro & Acrivos (Reference Marchioro and Acrivos2001), Mauri (Reference Mauri2003), Leshansky & Brady (Reference Leshansky and Brady2005) and Leshansky, Morris & Brady (Reference Leshansky, Morris and Brady2008), for example, for discussions on the collective shear-induced diffusivity.

Three different ways of measuring self-diffusivity are reported in the literature. The simplest one is to determine self-diffusion from the mean square displacement using the Einstein–Smoluchowski relation in one dimension,

where the brackets $\langle ~\rangle$ denote an ensemble average over available time intervals of length $t$ . Equivalently, self-diffusivity can be probed through the integral over the velocity autocorrelation function $\text{VACF}=(1/N)\sum _{i=1}^{N}\langle v_{i}^{y}(t^{\prime }),v_{i}^{y}(0)\rangle$ :

A third approach that has been used for measuring shear-induced diffusion in suspensions involves the self-dynamic structure factor (Leshansky & Brady Reference Leshansky and Brady2005; Leshansky et al. Reference Leshansky, Morris and Brady2008).

Figure 11 demonstrates the various ways of measuring the self-diffusivity for a system consisting of 512 particles at volume fraction $\unicode[STIX]{x1D719}=0.1$ . In order to remove system size effects, we measure the self-diffusivity using the integral over the velocity autocorrelation function for particle numbers $N$ from $64$ to $1024$ for fractions from $\unicode[STIX]{x1D719}=0.1$ , $0.2$ and $0.4$ , and up to $4096$ for $\unicode[STIX]{x1D719}=0.3$ . Assuming that there is a thin layer at the boundary that has a thickness proportional to the particle radius, and that does not contribute to the self-diffusion coefficient because the $y$ -component of velocity is small, we expect that the system size effect in the diffusion coefficient will scale as $N^{-1/2}$ . This allows us to extrapolate towards a large particle number as shown in figure 11(d). Figure 12 shows $D_{s}(\unicode[STIX]{x1D719})$ extrapolated towards infinite resolution for volume fractions up to $\unicode[STIX]{x1D719}=0.4$ . The best fit for a power law yields $D_{s}(\unicode[STIX]{x1D719})\sim \unicode[STIX]{x1D719}^{1.38\pm 0.50}$ . We do not expect this relation to hold for volume fractions that approach the jamming threshold where previous studies suggested that $D_{s}$ will flatten out at high volume fractions (Breedveld et al. Reference Breedveld, Van Den Ende, Bosscher, Jongschaap and Mellema2002; Metzger et al. Reference Metzger, Rahli and Yin2013b ), and, as the system jams, vertical motion is inhibited. However, simulating large volume fractions using finite elements for strains large enough to assess self-diffusion remains a very challenging task numerically. In addition we expect that particle interactions such as friction become increasingly important as the volume fraction increases (Boyer, Guazzelli & Pouliquen Reference Boyer, Guazzelli and Pouliquen2011; Seto et al. Reference Seto, Mari, Morris and Denn2013).

Table 1.

Shear-induced self-diffusivity compared to three-dimensional simulations from Sierou & Brady (Reference Sierou and Brady2004) where they used accelerated Stokesian dynamics (ASD). The results from Sierou & Brady match experimental measurements fairly well, although there are large variations of self-diffusivities measured from experiments. Our measured finite element method (FEM) results are systematically higher, which is expected for a two-dimensional system.

In table 1, $D_{s}$ obtained from our simulations are shown along with one example from Sierou & Brady (Reference Sierou and Brady2004) in three dimensions using accelerated Stokesian dynamics simulations. Qualitatively, our results show similar $\unicode[STIX]{x1D719}$ dependence as their results, although the self-diffusion coefficients from two-dimensional simulations are significantly larger than those obtained from three-dimensional simulations. While the self-diffusivity is larger in two dimensions than it is in three dimensions, the underlying mechanisms are the same.

Figure 11.

(a) Mean square displacement in the vertical direction of $N=512$ particles in simple shear at volume fraction $\unicode[STIX]{x1D719}=0.1$ . Initially, the displacement is ballistic and goes as $(\dot{\unicode[STIX]{x1D6FE}}t)^{2}$ , before it becomes diffusive at $\dot{\unicode[STIX]{x1D6FE}}t\sim 20$ and goes as $\dot{\unicode[STIX]{x1D6FE}}t$ . The self-diffusivity is determined from the slopes of either of the two curves. (b) Velocity autocorrelation function ( $v_{y}$ ) for 512 particles at volume fraction $\unicode[STIX]{x1D719}=0.1$ : grey lines, different simulations; blue (darker) line, ensemble average. (c) Integral of velocity autocorrelation function ( $v_{y}$ ) for $N=512$ particles at volume fraction $\unicode[STIX]{x1D719}=0.1$ : grey lines, different simulations; blue (darker) line, ensemble average. (d) Extrapolation $N\rightarrow \infty$ to remove system size effects measured with the integral over the velocity autocorrelation function. For $N\rightarrow \infty$ , $D_{s}/\dot{\unicode[STIX]{x1D6FE}}r^{2}$ is $0.031\pm 0.008$ . All data in this figure are for the solid phase, and are scaled upon the radius $r$ and the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ for dimensionless quantities.

Figure 12.

The self-diffusion coefficient $D_{s}$ scaled with the shear rate $\dot{\unicode[STIX]{x1D6FE}}$ and the particle radius $r$ as a function of particle volume fraction $\unicode[STIX]{x1D719}$ extrapolated towards $N\rightarrow \infty$ . The black broken lines show best fits with $D_{s}$ scaling as $\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FC}}$ . Error bars show 95 % confidence bounds

7 Discussion and conclusion

In this paper we introduced a two-dimensional adaptive unstructured finite element method for linear incompressible flow of suspensions at zero Reynolds number, and used it to study mixing in sheared suspensions. Direct numerical modelling of particle suspensions using the finite element model is challenging, and, in particular, close particle pairs are a fundamental problem for several reasons. First, particles can get arbitrarily close without touching, which poses issues about the resolution between the particles. Second, large gradients in both velocities and stresses are often associated with close particle pairs, so that it is important to obtain a decent resolution in these strategic areas. Particle overlaps pose a big challenge when modelling dense systems. And even though mesh refinement in strategic areas helps, it does not solve the problem, since infinitesimal distances between particles are not inhibited by Stokes equations. Still, mesh refinement allows us to run fairly dense systems to large enough strains to perform statistical analyses when combined with ensemble averaging.

If a dye is introduced into a sheared suspension, complex structures form in the fluid phase. Starting from thin filaments placed in the fluid phase, we characterized these structures using the diffusive strip method, and the structures show close resemblance to recent structures observed in experiments (Souzy et al. Reference Souzy, Yin, Villermaux, Abid and Metzger2015). As demonstrated in figures 4 and 8, from our simulations, we observe complicated structures developing with increased strain. The fractal dimension of these structures quickly increase towards $(F-1)/(1-\unicode[STIX]{x1D719})=1$ , and is governed by the mean square displacement in the fluid phase as well as the particle volume fraction.

At infinite Péclet number, we found that the shear-induced self-diffusion coefficient of the fluid phase is larger than for the solid phase. This observation has recently been observed in experiments (Souzy et al. Reference Souzy, Pham and Metzger2016), and we hypothesize that this difference is related to the rotational degree of freedom of the particles. One possible way to investigate further the relation between particle rotations and enhanced diffusivity in the fluid phase would be through simulations with full slip boundary conditions at the particle rims, which should remove effects of particle rotation in the motion of the fluid phase.

In reality, the shear-induced self-diffusion is a form of forced convection and has to be combined with the thermal or mass diffusion in order to assess the transport properties of the sheared suspension. At finite Péclet number, the differences between the fluid and the solid phase are further amplified. The differences in shear-induced self-diffusivity of the solid and the fluid phases has consequences for the way we interpret the relation between enhanced heat and mass transport of sheared suspensions and shear-induced self-diffusion as was done by Metzger et al. (Reference Metzger, Rahli and Yin2013b ). Although the shear-induced self-diffusion coefficients of the solid and the fluid phases are strongly coupled, our results suggest that the velocity field in the fluid phase needs to be determined in order to accurately correlate shear-induced self-diffusion with enhanced thermal and mass diffusivity. While particle trajectories are more readily available in experiments, particle image velocimetry measurements giving the fluid velocity field as was done by Souzy et al. (Reference Souzy, Yin, Villermaux, Abid and Metzger2015) to characterize fluid transport close to the boundary look promising.

To compare this to the self-diffusion coefficients of the fluid phase from figure 12, we need some constraint for $\dot{\unicode[STIX]{x1D6FE}}$ and $r$ . Typical numbers reported for slowly sheared systems in the literature are $\dot{\unicode[STIX]{x1D6FE}}\approx 1~\text{s}^{-1}$ (Metzger et al. Reference Metzger, Rahli and Yin2013b ), while particle sizes can be in the range $r\in [1~\unicode[STIX]{x03BC}\text{m},~1~\text{mm}]$ . While we do expect the self-diffusivity to be lower in three dimensions, these assumptions would result in $D_{s}\in [10^{-13}~\text{m}^{2}~\text{s}^{-1},~10^{-7}~\text{m}^{2}~\text{s}^{-1}]$ , which for millimetre-sized particles is approximately of the same order as the thermal diffusivity of water (the thermal diffusivity of water at room temperature is $D_{\text{H}_{2}\text{O}}\simeq 1.5\times 10^{-7}~\text{m}^{2}~\text{s}^{-1}$ (James Reference James1968)). For enhanced mass transport, we can compare the shear-induced self-diffusivity to the diffusion coefficients of, for example, methane and carbon dioxide in water, which are $D_{\text{CH}_{4}}\simeq 1.5\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$ and $D_{\text{CO}_{2}}\simeq 2\times 10^{-9}~\text{m}^{2}~\text{s}^{-1}$ , respectively (Cussler Reference Cussler2009). For rhodamine dye in water, the diffusion coefficient is $D_{\text{rhodamine}}\simeq 4.14\times 10^{-10}~\text{m}^{2}~\text{s}^{-1}$ (Culbertson, Jacobson & Ramsey Reference Culbertson, Jacobson and Ramsey2002), and for rhodamine dye in the mixture used by Souzy et al. (Reference Souzy, Yin, Villermaux, Abid and Metzger2015), the diffusion coefficient is as low as $D\simeq 10^{-13}~\text{m}^{2}~\text{s}^{-1}$ . This means that, to enhance mass transport in sheared suspensions, one can use even lower shear rates or smaller particles than one would need in order to significantly enhance thermal transport.

Here, we have found that the enhanced transport fits the functional form $D_{s}/D=1+\unicode[STIX]{x1D6FD}\unicode[STIX]{x1D719}^{\unicode[STIX]{x1D6FC}}Pe^{\unicode[STIX]{x1D701}}$ with $\unicode[STIX]{x1D6FD}=2.98\pm 0.39$ , $\unicode[STIX]{x1D6FC}=1.61\pm 0.26$ and $\unicode[STIX]{x1D701}=0.900\pm 0.031$ for $\unicode[STIX]{x1D719}\in [0.1,0.4]$ and $Pe\in [5,10^{4}]$ . The functional dependence that we find is similar to recent results for enhanced heat transport in suspensions of spheres (Metzger et al. Reference Metzger, Rahli and Yin2013b ), as well as experimental results for enhanced mass diffusivity in suspensions of spheres for large $Pe$ that have reported $D_{s}/D\sim Pe^{0.89}$ (Wang & Keller Reference Wang and Keller1985).

In microfluidic experiments, the Reynolds number is usually very small, and mixing can be difficult (Stroock et al. Reference Stroock, Dertinger, Ajdari, Mezić, Stone and Whitesides2002; Lee et al. Reference Lee, Chang, Wang and Fu2011). Using particles as mixers is therefore a possible option to enhance mass and thermal transport in microfluidic devices. Workamp, Saggiomo & Dijksman (Reference Workamp, Saggiomo and Dijksman2015) used the concept of enhanced transport properties in sheared suspensions to create a simple microfluidic mixer where the pressure drop across the mixer is low. This microfluidic mixer consists of a cylindrical chamber filled with particles, and a shear motion is set up using a moving magnet. As demonstrated in this paper, the enhancement of mass transport in such systems can be huge.

An additional aspect has to be taken into account if the mixing occurs in a system with a non-homogeneous $\dot{\unicode[STIX]{x1D6FE}}$ , as in Poiseuille flow, for example. In systems with a non-zero shear-induced collective diffusion coefficient, which scales with $\dot{\unicode[STIX]{x1D6FE}}$ , particles would migrate from regions of high shear stress towards regions of lower shear stress (Leighton & Acrivos Reference Leighton and Acrivos1987b ). In the case of Poiseuille flow, this would result in migration towards the centre of the channel for a suspension initially distributed uniformly in space (Boschan, Aguirre & Gauthier Reference Boschan, Aguirre and Gauthier2015). This change of particle density could potentially lower the shear-induced self-diffusion coefficient close to the boundary because the particle concentration decreases there.

In this paper we have studied the fully hydrodynamic limit. This means that the Reynolds number is zero, the particle Péclet number is infinite and there are no additional interactions between the particles such as normal contact, friction or repulsive forces. While we have simulated the purely hydrodynamic limit, several works have demonstrated the importance of particle contacts (e.g. Blanc, Peters & Lemaire Reference Blanc, Peters and Lemaire2011; Seto et al. Reference Seto, Mari, Morris and Denn2013; Gallier et al. Reference Gallier, Lemaire, Peters and Lobry2014; Pham et al. Reference Pham, Metzger and Butler2015, Reference Pham, Butler and Metzger2016). Not only do particle contacts alter the suspension microstructure (Blanc et al. Reference Blanc, Peters and Lemaire2011), for periodically sheared suspensions there is an onset of irreversible particle trajectories that occurs at a critical strain amplitude which correlates strongly with particle roughness Pham et al. (Reference Pham, Butler and Metzger2016). In periodically sheared suspensions, neither lubrication interactions Metzger & Butler (Reference Metzger and Butler2010) nor long-range hydrodynamic interactions Metzger, Pham & Butler (Reference Metzger, Pham and Butler2013a ) are the source of irreversibility. The irreversible behaviour is dominated by particle contact. It is worth noting that, while there is no direct source of irreversibility in our model, since the system is chaotic and all numerical models are subject to numerical errors, numerical error is bound to be a source of irreversibility at large strains.

Frictional interactions, repulsive forces, finite particle Péclet numbers and non-zero Reynolds numbers can all break the fore–aft symmetry that is present in the purely hydrodynamic limit (Bossis & Brady Reference Bossis and Brady1984; Brady & Morris Reference Brady and Morris1997; Metzger et al. Reference Metzger, Pham and Butler2013a ; Haddadi & Morris Reference Haddadi and Morris2014; Pham et al. Reference Pham, Metzger and Butler2015). Breaking of fore–aft symmetry will probably affect the motion of the particles, increasingly so as the particle volume fraction increases and the frequency of particle collisions increases. We also find it likely that any of these additional particle interactions could influence the shear-induced self-diffusion coefficient that we measure, both in the solid and in the fluid phase. We leave the investigation of the effect of additional particle interactions on shear-induced self-diffusion for future studies.