3.1 Hydrodynamic dispersion and velocity statistics

For hydrodynamic dispersion, the mean-square displacements (MSD) versus time are calculated to quantify the hydrodynamic diffusion of a suspension undergoing DIP. Initially, the MSD curve exhibits a quadratic growth with time. After a few particle–particle interactions, a transition to the diffusive regime takes place for the particle motion with linear growth of the mean-square displacement. In the inset of figure 2(a), the log–log plot of mean-square displacements at $\unicode[STIX]{x1D719}=50\,\%$ is presented to clearly show the diffusive behaviour in concentrated regimes, where a transition from a ballistic regime ( ${\sim}t^{2}$ ) to a diffusive regime ( ${\sim}t$ ) is observed. The average slopes of the MSD curves in the diffusive regime indicates an effective hydrodynamic diffusivity tensor $\unicode[STIX]{x1D63F}$ . It was previously reported that the hydrodynamic diffusion is primarily a consequence of ICEP interactions causing rearrangement of the particle configuration with no long-lasting clusters, while DEP interactions result in the formation of stable chains and trapping the particles in larger structures, leading to a subdiffusion (Park & Saintillan Reference Park and Saintillan2011a ,Reference Park and Saintillan b ).

Figure 2(a) shows the dependence of the diffusivity on the volume fraction from dilute to concentrated regimes. For the error bars on the plot, the block-averaging method is used to calculate the standard errors of the time-averaged quantity (Flyvbjerg & Petersen Reference Flyvbjerg and Petersen1989). Notably, a non-trivial behaviour is observed in concentrated regimes, where a secondary peak arises at $\unicode[STIX]{x1D719}\approx 50\,\%$ . As shown in Park & Saintillan (Reference Park and Saintillan2010), for a dilute suspension, the hydrodynamic diffusivity increases with the volume fraction before reaching a (primary) maximum, and then significantly decreases in the semi-dilute regime up to $\unicode[STIX]{x1D719}=20\,\%$ . The current simulation results show that it continues to decrease and then reaches the local minimum at a volume fraction of about 35 %. Interestingly, the hydrodynamic diffusivity starts to increase again in a concentrated regime until it reaches a local (secondary) maximum at a volume fraction of about 50 %, and then drastically decreases on approaching random close packing. It is worth noting that the tracer diffusion (not shown) is found to show a trend qualitatively similar to that of figure 2(a).

Figure 2.

(a) Hydrodynamic diffusivities in the transverse direction ( $x$ ) and the field direction ( $z$ ) versus volume fraction, with the error bars based on the block-averaging method (Flyvbjerg & Petersen Reference Flyvbjerg and Petersen1989). Inset: mean-square displacement (MSD) curves in the transverse ( $x$ , $y$ ) and field direction ( $z$ ) at $\unicode[STIX]{x1D719}=50\,\%$ as functions of time on a log–log scale. (b) Standard deviations of particle velocities as functions of the volume fraction in the $x$ and $z$ directions.

In order to further characterize the suspension dynamics, a statistical analysis is carried out for particle velocities. Since the velocity of the suspended particles only comes from the particle–particle interactions due to DIP, there is no net particle motion – the mean particle velocity is exactly zero at every time step. Therefore, the particle velocity statistics can be readily characterized by the standard deviation of velocity distributions. Figure 2(b) presents the standard deviation of particle velocity distributions in the transverse $(x)$ and field $(z)$ directions as a function of the volume fraction. As seen in the figure, the magnitude of the velocity fluctuations increases with the volume fraction, reaching a maximum in the dilute regime ( $\unicode[STIX]{x1D719}\approx 8\,\%$ ), and then decreases in semi-dilute suspensions. At $\unicode[STIX]{x1D719}\approx 30\,\%$ , however, the velocity fluctuation starts to increase again, reaching a local peak at $\unicode[STIX]{x1D719}\approx 45\,\%$ , and then decreases with the volume fraction. The velocity fluctuation indeed exhibits a similar trend to the diffusivity (figure 2 a). In addition, it is seen that the fluctuation in the transverse direction starts to increase earlier (at $\unicode[STIX]{x1D719}\approx 25\,\%$ ) than that in the field direction (at $\unicode[STIX]{x1D719}\approx 30\,\%$ ). This observation suggests that suspension dynamics starts to change with regard to the preferred direction of particle motions at $\unicode[STIX]{x1D719}=25\,\%\sim 30\,\%$ .

The secondary peaks in the hydrodynamic diffusivity and velocity fluctuation are counter-intuitive since it is expected that the mobility of particles tends to decrease with concentration due to a consequence of strong excluded volume interactions. This distinct dynamics is most easily seen in supplementary movies 2 and 3 for $\unicode[STIX]{x1D719}=25\,\%$ and 45 %, where the latter displays faster dynamics than the former. This non-trivial behaviour can result from several reasons. Clarification of the mechanisms contributing to the secondary peak will be given as the paper proceeds.

3.2 Suspension microstructure

The change in the suspension dynamics can be investigated quantitatively by calculating the pair distribution function $p(r,z)$ , revealing the local microstructure of the suspension. This function provides the probability of finding the particle at a location $(r,z)$ in cylindrical coordinates ( $r^{2}=x^{2}+y^{2}$ ) with respect to a probe particle located at the origin. The functions are shown in figure 3 for various volume fractions. While dilute and semi-dilute suspensions ( $\unicode[STIX]{x1D719}=10\,\%$ and 20 %) show the maximum probability of particles paired up in the field direction (near the pole of the particles), which has also been observed in suspensions of spherical and rod-like particles (Saintillan et al. Reference Saintillan, Darve and Shaqfeh2006; Rose et al. Reference Rose, Hoffman, Saintillan, Shaqfeh and Santiago2009; Park & Saintillan Reference Park and Saintillan2010), the peak location shifts from the pole to the equator in the concentrated regime ( $\unicode[STIX]{x1D719}>30\,\%$ ). The maximum probability region then seems to emanate and propagate towards the pole on further increasing the volume fraction. Finally, the peak region turns to entirely cover the particle surface at $\unicode[STIX]{x1D719}=59\,\%$ . The suspension microstructure undergoes a drastic change for $\unicode[STIX]{x1D719}>30\,\%$ , which might infer the non-trivial behaviours observed in the diffusivities and velocity fluctuations.

Figure 3.

Pair distribution functions in the suspensions undergoing dipolophoresis at six different volume fractions in cylindrical coordinates $(r,z)$ , where $r^{2}=x^{2}+y^{2}$ . A probe particle is located at the origin.

In addition to the pair distribution function, particle occupancy statistics can be presented to characterize the number-density fluctuations for non-uniform microstructure on a larger length scale. To calculate the particle occupancy distribution, a cubic cell with a fixed volume $V$ , which is calculated based on the fixed average number of the particles $\langle N\rangle =\unicode[STIX]{x1D719}V/V_{p}$ , where $V_{p}$ is the volume of one sphere, is placed arbitrarily inside the simulation box. The cell will then contain $N$ number of particles that can differ from the expected value $\langle N\rangle$ . The probability or distribution of the number of particles $P(N)$ in the fixed cell volume presents the number-density fluctuation of the suspension. Figure 4(a) shows the particle occupancy distribution for a cell containing $\langle N\rangle =15$ particles on average. The curves include the occupancy distributions at steady state for different volume fractions, along with a Poisson distribution given by

(3.1) $$\begin{eqnarray}P(N)=\frac{\langle N\rangle ^{N}\text{e}^{-\langle N\rangle }}{N!}.\end{eqnarray}$$

It is shown that none of the distribution curves for $\unicode[STIX]{x1D719}>15\,\%$ are captured by the Poisson distribution, owing to strong excluded volume interactions in semi-dilute and concentrated suspensions as opposed to a dilute suspension at $\unicode[STIX]{x1D719}=2.5\,\%$ . To make a quantitative comparison, the standard deviation $\unicode[STIX]{x1D70E}_{N}$ of the particle occupancy distribution can be calculated for a range of volume fractions. As is the case in non-Poisson statistics (Bergougnoux & Guazzelli Reference Bergougnoux and Guazzelli2009), a power law of $\unicode[STIX]{x1D70E}_{N}$ as a function of $\langle N\rangle$ can effectively capture a number-density deviation from the Poisson distribution $\unicode[STIX]{x1D70E}_{N}=\langle N\rangle ^{1/2}$ . It is found that all curves for $\unicode[STIX]{x1D719}>15\,\%$ are well-fitted by a power law $\unicode[STIX]{x1D70E}_{N}\approx \langle N\rangle ^{0.36}$ , which is a quite substantial deviation from the Poisson law. It suggests a significant number-density deviation from a random distribution at larger scales.

Figure 4.

(a) Particle occupancy distributions for $\langle N\rangle =15$ at different volume fractions along with a Poisson distribution. (b) The standard deviations of number-density fluctuations normalized with their maximum value as functions of volume fraction for three expected values: $\langle N\rangle =5$ , 10 and 20.

As seen in inset of figure 4(a), the distributions at $\unicode[STIX]{x1D719}=15\,\%$ and 45 % are broader compared to those at $\unicode[STIX]{x1D719}=25\,\%$ and 59 %. Specifically, it indicates that the number-density fluctuation initially decreases, then increases, and finally decreases with the volume fraction. To clearly illuminate this trend, figure 4(b) shows the normalized standard deviation $\widetilde{\unicode[STIX]{x1D70E}}_{N}$ of number-density fluctuations as a function of volume fraction for three average numbers $\langle N\rangle =5$ , 10 and 20. For all three cases, they share a trend similar to that of the diffusivities in figure 2(b). After a local minimum in the semi-dilute regime ( $\unicode[STIX]{x1D719}\approx 25\,\%$ ), the standard deviation of the number-density fluctuation increases with the volume fraction in concentrated suspensions, and then decreases on approaching random close packing. This similar behaviour indicates that the fluctuation of the number of adjacent particles is strongly related to the hydrodynamic diffusion of the DIP suspension – the larger the number-density fluctuations, the stronger the hydrodynamic diffusion in the suspension.

3.3 Dynamics under dipolophoresis in concentrated regimes

We now attempt to illuminate the mechanisms behind the non-trivial behaviours in concentrated regimes. It has been known that in dense suspensions, particle interactions undergoing a lubricated-to-frictional transition provide a coherent mechanistic basis for the shear thickening in these suspensions (Morris Reference Morris2018). Mari et al. (Reference Mari, Seto, Morris and Denn2014) performed simulations for non-Brownian suspensions under shear, and showed that a growing number of frictional contacts and its intermittent occurrence contribute to discontinuous shear thickening. Similarly, for concentrated DIP suspensions, particle contacts also appear to provide a coherent mechanistic basis for the non-trivial behaviours observed in concentrated regimes.

Figure 5.

(a) Radial components of the total DIP relative velocity between two spheres as functions of the angle $\unicode[STIX]{x1D6E9}$ between the direction of the electric field $\boldsymbol{E}_{0}$ and separation vector $\boldsymbol{R}$ for different separation distances. Velocities are scaled by $U_{0}=(\unicode[STIX]{x1D700}a/\unicode[STIX]{x1D702})E_{0}^{2}$ . Inset: the magnitude ratio of the radial relative velocities in the transverse direction ( $\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2$ ) and in the field direction ( $\unicode[STIX]{x1D6E9}=0$ ) as a function of separation distance $R/a$ – the transverse velocity becomes dominant at $R/a\leqslant 2.23$ . (b) The average minimum separation distance of particles in a suspension as a function of volume fraction, along with a dashed line of $R/a=2.23$ .

The suspension dynamics undergoing both DEP and ICEP or DIP is originated from pairing dynamics through electric and hydrodynamic interactions, contributing to particle contacts. For the case of a two-particle interaction, a radial component of the relative velocity $\boldsymbol{U}=\boldsymbol{U}_{2}-\boldsymbol{U}_{1}$ can provide the nature of particle contacts. Figure 5(a) shows the radial relative velocity $U_{r}=\boldsymbol{U}\boldsymbol{\cdot }\hat{\boldsymbol{R}}$ for different separation distances between two particles as a function of the angle $\unicode[STIX]{x1D6E9}$ between the axis of the two particles and the direction of the external field. The particles are attracted ( $U_{r}<0$ ) when aligned with the direction of the electric field, while they are repulsive ( $U_{r}>0$ ) when aligned in the transverse direction. At a large separation distance (e.g., $R/a=2.5$ ), the field-direction velocity $U_{r}(\unicode[STIX]{x1D6E9}=0)$ is greater than the transverse-direction velocity $U_{r}(\unicode[STIX]{x1D6E9}=\unicode[STIX]{x03C0}/2)$ and a large portion of the curve is on the attractive side. These results suggest that attractive contacts are likely to happen between two particles. As the particles come much closer (e.g., $R/a<2.2$ ), however, the transverse-direction velocity is much greater the field-direction velocity, and most of the curve is on the repulsive side. In such a case, the particles are strongly separated in the traverse direction and repulsive contacts are likely to happen if another particle is nearby – this contact is highly expected to happen for concentrated suspensions. In the inset of figure 5(a), the ratio of the transverse-direction velocity to the field-direction velocity is shown, where the transverse velocity becomes predominant at $R/a\leqslant 2.23$ . At $R/a=2.01$ , the repulsive velocity is almost four times the attractive velocity. To relate the changeover distance $R/a=2.23$ to a characteristic separation distance of the suspension, the time-averaged minimum separation distance of all particles in a suspension is calculated at different volume fractions, as seen in figure 5(b). At $\unicode[STIX]{x1D719}\approx 21\,\%$ , the average minimum distance is about $R/a=2.23$ , at which the dominant direction of the pairing dynamics is expected to change from the field direction to the transverse direction. At higher volume fractions ( $\unicode[STIX]{x1D719}>45\,\%$ ), the average minimum distance of a suspension is almost $R/a=2.01$ , where repulsive forces are quite strong in the transverse directions. These forces are expected to strongly push particles into each other along the transverse directions to make them contact nearby neighbours very often in concentrated suspensions.

Figure 6.

(a) The time-averaged number of total contacts, attractive contacts and repulsive contacts, as functions of the volume fraction. (b) The contact ratio $f$ , which is the fraction of the number of repulsive contacts to total contacts, as a function of the volume fraction.

Finally, we can make the link between particle contacts (excluded volume interactions) and the non-trivial behaviours in concentrated regimes. Given the contact direction, the particle contacts can be classified into attractive and repulsive classes with respect to an angle $\unicode[STIX]{x1D6E9}$ at which a contact happens. We used the critical angle $\unicode[STIX]{x1D6E9}_{c}=0.11\unicode[STIX]{x03C0}$ at which $R/a=2.075$ – the issue of sensitivity to the chosen angle was tested using different angles at shorter distances, giving essentially identical results. The contact is called a repulsive contact if $\unicode[STIX]{x1D6E9}>\unicode[STIX]{x1D6E9}_{c}$ or an attractive contact if $\unicode[STIX]{x1D6E9}<\unicode[STIX]{x1D6E9}_{c}$ . Figure 6(a) shows the number of total, attractive and repulsive contacts as a function of volume fraction. In a dilute regime ( $\unicode[STIX]{x1D719}<10\,\%$ ), there are more attractive contacts than repulsive contacts – attractive interactions are the dominant mechanism for pairing dynamics. For $\unicode[STIX]{x1D719}>16\,\%$ , however, most of the contacts are from repulsive ones, indicating that repulsive interactions are the dominant mechanism of pairing dynamics that causes pushing particles to contact. To further manifest the dominance of repulsive contacts, we calculate the fraction of repulsive contacts $f(\unicode[STIX]{x1D719})$ , which is defined as the ratio of the number of repulsive contacts divided by the total number of contacts. Figure 6(b) shows that there is a sharp increase in $f$ from 0.26 ( $\unicode[STIX]{x1D719}=10\,\%$ ) to 0.90 ( $\unicode[STIX]{x1D719}=25\,\%$ ). Then, almost 95 % of total contacts are repulsive in the range $30\,\%<\unicode[STIX]{x1D719}<40\,\%$ . These strong, massive repulsive contacts are direct evidence as to the non-trivial behaviours in concentrated regimes – particles tend to move very fast and meet others in the repulsive zone. The number of repulsive contacts is still around 85 % of the total contacts up to random close packing. These results are indeed related to the pair distribution functions, where the sharp peak probability region is located at the equator region and eventually covers the particle surface in concentrated regimes.