The rheological behaviors of suspension of ideally conductive particles in an electric field are studied using large-scale numerical simulations in the limit of the zero-shear-rate flow. Under the action of an electric field, the particles undergo the nonlinear electrokinetic phenomenon termed dipolophoresis, which is the combination of dielectrophoresis and induced-charge electrophoresis. For ideally conductive particles, the dynamics of the suspension are primarily controlled by induced-charge electrophoresis. To characterize the rheological properties of the suspension, the particle stress tensor and particle pressure are calculated in a range of volume fraction up to almost random close packing. The particle normal stress and particle pressure are shown to behave nonmonotonically with volume fraction, especially in concentrated regimes. In particular, the particle pressure is positive for volume fraction up to 30%, after which it becomes negative, indicating a change in the nature of the particle pressure. The microstructure expressed by the pair distribution function and suspension entropy is also evaluated. Visible variations in the local microstructure seem to correlate with the nonmonotonic variation in the particle normal stresses and particle pressure. These nonmonotonic behaviors are also correlated with the change in the dominant mechanism of particle pairing dynamics observed in our recent study [Mirfendereski and Park, J. Fluid Mech. 875, R3 (2019)]. Finally, the effects of confinement on the particle stress and particle pressure are investigated. It is found that the particle pressure changes its nature very quickly at high volume fractions as the level of confinement increases. This study should motivate control strategies to fully exploit the distinct changing nature of the pressure for rheological manipulation of such a suspension system.

An electric-field-driven suspension of particles in a viscous fluid has been widely studied in a number of fields including material science, microfluidics and nanofluidics, and bioengineering [1–4]. In general, there are two classes in such systems. The most notable class is the so-called electrorheological (ER) suspension in which nonconductive but electrically active particles are suspended in an electrically insulating fluid. A wide spectrum of industrial applications for ER suspensions has been enabled by a rich rheological property [5]. The other class is comparatively new, where conductive or polarizable particles are suspended in an electrically conductive fluid or electrolyte. Recently, such suspension has gained a growing interest in the additive manufacturing of 3D printing technology and electrochemical energy storage technology [6–8]. As highly concentrated suspensions are commonly used for these technologies, understanding the rheological behaviors of such suspensions is of practical interest. However, the studies on the rheology of such suspensions remain limited [9,10] to which this work is intended. The dynamics of such suspensions is governed by the so-called induced-charge electrophoresis (ICEP) for ideally polarizable particles [11–13]. Recently, such suspension displays nontrivial behaviors in concentrated regimes [11].

In the class of ER suspensions, electronic control of stress transfer tends to establish the concept of distinct fluid types. These suspensions are sometimes denoted as smart fluids, leading to various applications, including active shock absorbers, clutches, brakes, dampers, actuators, and artificial joints [2,14]. The ER fluids under the action of an electric field are known to undergo dielectrophoresis (DEP) and exhibit a dramatic viscosity enhancement, which is reversible and can be controlled by the electric field [2,15,16]. This increase in viscosity, potentially leading to a transition to solid-state, is known as a consequence of a rapid formation of particle chains and columns along the applied field direction due to dipolar interactions between particles [17–20]. These fibrillated chain structures are regarded as a typical feature of ER fluids and were first found by Winslow [21]. Above a limiting volume fraction, these fibrillated structures, which span the electrode gap in quiescent suspensions, tend to become larger with particle concentration and eventually degrade with increasing shear rate in a shear flow [22]. It has also been well known that the unique rheological properties of the ER fluids are attributed to the resistance to deformation of the field-induced structures [22,23].

In a similar fashion, the presence of external magnetic fields in a suspension of highly magnetizable particles can also lead to a significant change in apparent viscosity and a reversible transition from liquid to nearly solid-state [24,25]. Such suspension is also termed a magnetorheological (MR) fluid. The viscosity of MR fluids, which contain the particles with a typical diameter of 1–10 μm, can also be tuned with the field strength. As opposed to MR fluids, ferrofluids are composed of much smaller magnetic particles and exhibit a very small change in viscosity in the presence of an external field due to the considerable Brownian effects, disrupting the formation of the structures [26].

Rheological properties due to dipolar interactions have received much attention in both fields of electrorheology and magnetorheology [27,28]. Sim et al. studied the large amplitude oscillatory shear (LAOS) behavior of ER fluids, relating the nonlinear LAOS behavior to the microstructural change of the suspension [29]. They also found that the cluster formation along with a slight rearrangement within a cluster results in the strain overshoot phenomenon, often seen in the system of complex fluids. Bonnecaze and Brady observed different types of rearrangements at small Mason numbers (MNs), the ratio of viscous to electrostatic forces [23]. They correlated the macroscopic rheology to the dynamics of the suspension for a range of MNs. They also showed that the decrease in MN, which is associated with an increase in electrostatic force relative to the viscous force, leads to an increase in effective viscosity [23,30]. This MN is originally used for MR fluids, where it represents the ratio between the viscous and magnetic forces. Similar to ER fluids, the apparent shear viscosity of an MR suspension was seen to collapse onto a single function of MN for a range of conditions [31,32]. It is worth noting that Anderson claimed that the polarizability of particles could significantly affect the ER behavior of the suspension as conductive particles can acquire additional charges and, in turn, fail to participate in the formation of the structure as opposed to dielectric particles [33].

As for the other class, which is the case of interest for this study, conductive particles are suspended in a conducting fluid. The dynamics of the system becomes more complex and far apart from that of the conventional ER suspensions governed by DEP interactions. In this system, ICEP [34,35] arises as the particles acquire additional surface charges resulting in a nonlinear fluid flow around particle surfaces. This phenomenon was first analyzed in the Russian colloidal literature [36,37], which was reviewed by Murtsovkin [38], and later rediscovered and coined by Squires and Bazant [34]. It has been observed that ICEP is predominant over DEP for a suspension of ideally conductive particles [12], but this predominance can be modulated by surface contamination [39]. Recent numerical simulations were performed to investigate the effects of ICEP on the dynamics of such suspensions in a range of volume fractions up to a maximum of about 60%, where nontrivial behaviors in large-scale dynamics were discovered in concentrated regimes (volume fraction of 45–50%) and explained by the nature of contact mechanisms [11]. Yet, no rheological studies of such system have been done in the literature, which motivates the current study.

It should be noted that owing to qualitatively similar far-field fluid disturbances governed by a Stokes dipole, the system studied shares similarities with active suspensions such as a suspension of spherical squirmers driven with a prescribed axisymmetric tangential velocity [40,41]. Hydrodynamic interactions in active suspensions induced by permanent swimming dipoles result in complex dynamics and mixing and diffusive behaviors [42], which are qualitatively similar to ones of the current system governed by ICEP. However, the expected differences are the magnitude and orientation of the surface (or slip) velocity, which may modify the relative importance of attractive and repulsive interactions between particles [43,44]. In the context of the rheology of active suspensions, theoretical, numerical, and experimental studies have shown that particle activity causes an increase in intrinsic viscosity of the suspension of pullers, while reducing the intrinsic viscosity of the suspension of pushers when the particles are strongly active relative to the weak external flow (low Péclet numbers) [42,45–47]. It is also consistent with the observations that the suspension of living algae (pullers) is more viscous than the suspension of dead ones, while the suspension of both E. coli bacteria and sperms (pushers) is less viscous when alive than dead [48].

In this paper, large-scale numerical simulations are used to investigate the rheology of the suspension of noncolloidal, ideally conductive spheres in a uniform electric field for volume fractions up to random close packing in the limit of the zero-shear-rate flow. The suspensions under these conditions are known to undergo both ICEP and DEP. A combination of these two phenomena is also known as dipolophoresis (DIP) [49]. The particle-particle interactions arising from DIP are generally governed by ICEP in a suspension of ideally conductive particles since the leading-order contributions to ICEP and DEP interactions are a Stokes dipole and a potential quadrupole, respectively, where the former has a slower decay with separation distance R as O(R2) than the later as O(R4) [11–13]. Under these conditions, the particles undergo random chaotic motions, which result in the constant rearrangement of the particle configurations but do not lead to the formation of chains as in the case of DEP only. The details of the simulation method adopted to facilitate the study of highly concentrated suspensions are reported in Sec. II. In Sec. III, the simulation results for rheology and microstructure in the suspension as functions of volume fraction along with the system relaxation are presented and discussed. The effect of wall confinement on the rheology is also evaluated in Sec. III E.

We consider a suspension of N identical neutrally buoyant ideally polarizable spheres of radius a in a viscous electrolyte with the permittivity ε and viscosity η. A cubic periodic domain is used to simulate an unbounded, infinite suspension. An external uniform electric field E0=E0z^ is applied in the z direction. The particles are assumed to carry no net charge, so the linear electrophoresis is not expected to occur. We also assume weak electric fields, thin Debye layers, and zero Dukhin number for no surface conduction [34]. It is also assumed that the particle size is large enough so that the Brownian motion is negligible. Under these conditions and assumptions, the suspension dynamics results entirely from the effect of DIP.

As a uniform electric field is applied, electric and hydrodynamic interactions between particles arise as a result of DEP and ICEP, which may lead to relative motions of particles. For ICEP, each spherical particle polarizes and forms a nonuniform surface charge distribution, which then attracts counterions in an electrolyte. The migration and accumulation of these counterions near the polarized surface result in the formation of a nonuniform Debye layer. The charging time of this nonuniform Debye layer is very small on the order of τc=λDa/D, where λD is the Debye layer screening length, a is the particle radius, and D is the characteristic diffusivity of ions in solution. In a typical experiment (a10μm, λD10 nm, D105cm2 s1), the charging time τc104 s, which is much smaller than the diffusion time across the particle τa=a2/D101 s. Therefore, the nonuniform Debye layer can be assumed to remain at equilibrium. The effect of the electric field on the nonuniform Debye layer drives disturbance flows around the particle surface, which may lead to relative motions due to ICEP. The current calculation for ICEP interactions is based on the standard model of induced-charged electroosmosis (ICEO) around an ideally conductive sphere [34], which has limitations compared to experimental observations, such as electrostatic correlations, electroviscous effects, dielectric decrement, among others [50]. Thus, it should be noted that, in general, the various neglected conditions are likely to modify the magnitude of ICEP interactions and shift their frequency response compared to the standard model used in the current study. In addition, given that in the current model, Debye layers are sufficiently thin (λD/a103) and the minimum distance between particles in a suspension is larger than 2λD even at high volume fractions, it can also be assumed that the Debye layer overlapping is negligible. Hence, the capacitance or polarizability of particles remains at equilibrium even at high volume fractions. However, in the case of highly concentrated suspensions, the Debye layer overlapping might happen due to strong excluded volume interactions, leading to the modification of the Debye layer capacity [51], which is not considered in the current study. For DEP, the presence of other particles like a suspension causes disturbances to the local electric field around particles, resulting in a nonuniform Maxwell stress tensor O(E2) in the fluid. This tensor can yield a nonzero DEP force on surrounding particles, leading to relative motions in particle suspensions due to DEP [12].

In a uniform electric field, both ICEP and DEP do not lead to the motion of a single sphere owning to fore-aft symmetry. Indeed, the presence of other particles in a suspension leads to the symmetry breaking, resulting in the relative motion between the particles due to ICEP and DEP. The detailed description of paring dynamics associated with DIP for ideally conductive spheres in a uniform electric field was presented in the previous studies [12,13].

For electric and hydrodynamic interactions in a suspension, the method of simulation is based on the previous work of Park and Saintillan [12], for which a new approach was introduced to efficiently prevent particle overlaps in our recent work [11]. The outline of the method is as follows. Based on pair interactions due to DIP [13], the translational velocity x˙α of a given particle α in a suspension can be expressed by

x˙α=εaE02ηβ=1N[MDEP(Rαβ/a)+MICEP(Rαβ/a)]:z^z^+M0IFpα=1,,N,
(1)

where Rαβ=xβxα is the separation vector between particle α and particle β, and MICEP and MDEP are third-order dimensionless tensors accounting for the ICEP and DEP interactions, respectively. M0 is a mobility constant of a single sphere, I is an identity tensor, and Fp is the interparticle force for excluded volume interactions. It is shown that these two tensors MICEP and MDEP can be entirely determined by the scalar functions of the dimensionless inverse separation distance λ=2a/|R|. Specifically, for both DEP and ICEP, M is calculated as follows:

M={Mp(R/a)R/a4,Mp(R/a)MFF(R/a)+MTM(R/a)R/a<4,
(2)

where Mp denotes the periodic version of the far-field tensors MFF, which are given by

MFFDEP(R/a)=112T(R/a)+O(λ5),
(3)
MFFICEP(R/a)=98S(R/a)1124T(R/a)+O(λ5),
(4)

where the two tensors S and T=2S are the Green’s functions for a Stokes dipole and for a potential quadrupole, respectively [52]. These two third-order tensors are also given in index notations as follows:

Sijk(R)=1R3(δijRkδikRjδjkRi)3RiRjRkR5,
(5)
Tijk(R)=6R5(δijRk+δikRj+δjkRi)+30RiRjRkR7.
(6)

These far-field tensors are asymptotically valid to order O(λ4) for any pair of particles, and their use is consequently justified when they are sufficiently far apart (i.e., their separation distance is greater than 4a). The far-field interactions (λ1) can be readily computed using the method of reflections. These tensors reconfirm the dominance of hydrodynamic interactions due to ICEP over ones due to DEP in a suspension of ideally conductive particles. However, the near-field corrections are necessary as the method of reflections becomes inaccurate when particles are close to each other (typically |Rαβ|<4a), for instance, during pairing events. This is achieved by correcting the far-field tensor MFF with a more accurate version MTM calculated using the method of twin multiple expansions [13]. This method is very accurate down to separation distances on the order of |Rαβ|2.005a [12]. In addition, this method captures the strong modification of the local electric and hydrodynamic interactions between the two particles very accurately even at high volume fractions.

To account for hydrodynamic interactions between particles α and β and all its periodic images, the tensors of Eqs. (3) and (4) can be expressed as the periodic version of the far-field tensors, which are valid to order O(Rαβ4). Since a high-order computation O(N2) is required to direct calculation of the sums in Eq. (1), the smooth particle mesh Ewald (SPME) algorithm based on the Ewald summation formula of Hasimoto [53] and on fast Fourier transforms is used to accelerate the calculation of the sums to O(NlogN) operations [12]. The details of the SPME algorithm on calculations of particle velocities can be found in the  Appendix.

Once all the particle velocities x˙α (α=1,,N) are calculated, particle positions are advanced in time using a second-order Adams–Bashforth time-marching scheme, with an explicit Euler scheme for the first time step. A fixed time step Δt is used and is chosen so as to ensure that particles only travel a fraction of the mean interparticle distance during one integration step. In order to prevent particle overlaps that occur due to the use of finite time steps in simulations, the application of a repulsive interparticle force is necessary. For this purpose, an effective algorithm, functionally identical to the potential-free algorithm [54], is implemented to prevent excessive particle overlaps as a result of DIP interactions, where particles are moved almost exactly, within roundoff errors (2.005a), to contact. The form of the repulsive potential for excluded volume (EV) interactions is

UEV=12k(2aR)2,
(7)

where k is the time step-dependent prefactor, which can be expressed as k=3πηa/Δt [55]. A particle velocity driven by a short-range repulsive force, which is the negative gradient of the potential with respect to the coordinates of the particle, can be obtained by the Stokes drag law. This repulsive force corresponds to the interparticle force in Eq. (1). The resulting velocity contributes to the displacement along the direction connecting the center of the two spheres at the points of closest approach. We tested the robustness with respect to the excluded volume interactions by using different values of the prefactor in Eq. (7) and no excluded volume interactions, where almost identical results were produced. Another troublesome factor for suspension simulations in concentrated regimes is to provide the initial configuration of random particle distributions. To this end, the initial random configurations were generated using a similar procedure to ones suggested for dense hard-sphere systems [56–59]. In the present simulations, all runs were started with hard-sphere equilibrium configurations, but the first 100 time steps were discarded for better steady-state configurations beginning to compute averages.

In the remainder of the paper, all variables are made dimensionless using the following characteristic length and velocity scales: lc=a and uc=εaE02/η.

We performed large-scale simulations in a cubic periodic domain (Lx×Ly×Lz=203) for a range of volume fraction ϕ up to almost random close packing (ϕ64%). In our previous study, the nontrivial suspension dynamics was observed in concentrated regimes, where the velocity fluctuation, hydrodynamic diffusivity, and number density fluctuation tend to increase with volume fraction at ϕ=35 before reaching a local maximum at ϕ45% and then drop as approaching to the random close packing [11]. We attributed this nonmonotonic behavior to the change in the dominant direction of particle-particle contacts from the field direction to the transverse direction. In this study, the local microstructure, which is correlated with the suspension dynamics, is evaluated by a pair distribution function. The rheology of the system is characterized by computing suspension bulk stress. Particle extra stress tensor Σp, which provides the essential explanation of the mean particle effect on the flow, is calculated based on the Batchelor calculation of the average stress tensor in particle suspension [60]. The particle pressure as an isotropic part of the stress is evaluated for a range of volume fraction [61].

The mean-square displacement (MSD) of particles versus time, which generally has been served as a start of the pathway to quantify the collective dynamics of suspended particles, is presented for four different volume fractions in Figs. 1(a) and 1(b) for the z direction (field direction) and the x direction (transverse direction), respectively. After a few particle-particle interactions, the initial quadratic growth of the MSD curve is followed by a transition to the diffusive regime in which the MSD curve linearly grows with time in both x and z directions. To compute a time required to reach the diffusive regime due to loss of memory, denoted as a relaxation time (also known as crossover time), the autocorrelation functions Cuu of particle velocities in the x and z directions are calculated, as seen in Fig. 2. The relaxation time can be approximated as the time when the function reaches its first global minimum. In both directions, the function decorrelates very quickly with increasing volume fraction.

FIG. 1.

The mean-square displacement (MSD) in (a) the z direction (field direction) and (b) the x direction (transverse direction) as a function of time on a log–log scale at four different volume fractions.

FIG. 1.

The mean-square displacement (MSD) in (a) the z direction (field direction) and (b) the x direction (transverse direction) as a function of time on a log–log scale at four different volume fractions.

Close modal
FIG. 2.

The autocorrelation function of the particle velocity components in (a) the z (field) direction and in (b) the x (transverse) direction as a function of time for four different volume fractions.

FIG. 2.

The autocorrelation function of the particle velocity components in (a) the z (field) direction and in (b) the x (transverse) direction as a function of time for four different volume fractions.

Close modal

Figure 3 shows the relaxation time τ as a function of volume fraction ϕ on a log-log scale. In the dilute regime (ϕ<5%), the relaxation time tends to decrease with volume fraction and scale as τϕ0.5. Subsequently, in the semidilute regime (ϕ=535%), the relaxation time decreases faster than for the dilute regime, approximately scaling as τϕ1. The decreasing trend in the relaxation time with volume fraction is easily explained by an increase in the magnitude of particle-particle interactions with increasing volume fraction. Interestingly, the relaxation time appears to increase at ϕ35% up to ϕ47.5% and then decrease again as approaching random close packing. It is this range of volume fraction that the nontrivial behavior of increasing hydrodynamic diffusivity and velocity fluctuation was observed in our previous study for the same suspension [11].

FIG. 3.

The relaxation time in the transverse (x) and field (z) directions as a function of volume fraction on a log–log scale. Different zones associated with the different slopes are seen on a log–log scale for ϕ<35%, and a nontrivial variation in relaxation time is seen for ϕ>35%.

FIG. 3.

The relaxation time in the transverse (x) and field (z) directions as a function of volume fraction on a log–log scale. Different zones associated with the different slopes are seen on a log–log scale for ϕ<35%, and a nontrivial variation in relaxation time is seen for ϕ>35%.

Close modal

To further distinguish different behaviors in different ranges of volume fraction, the mean interparticle gap h¯=R¯min2a is calculated as a function of volume fraction, as shown in Fig. 4. Four different zones can be readily distinguished by different slopes in the log–log plot, namely, dilute, semidilute, concentrated, and very concentrated (close to random close packing) regimes. In dilute regime (first zone), the average interparticle gap is proportional to h¯ϕ1/3, which has also been observed in many dilute suspensions [62–65]. The exponent then becomes larger at ϕ5%, providing the second slope of 0.75 (second zone). In concentrated regime (third zone), it decreases dramatically with a much higher exponent as a volume fraction increases. Finally, the slope in the very high concentrated regime (fourth zone) decreases slowly compared to the third zone because the average interparticle gap is approaching the value of designated prefactor corresponding to excluded volume interaction in Eq. (7).

FIG. 4.

The average interparticle gap h¯/a=R¯min/a2 as a function of volume fraction on a log–log scale. Four different zones associated with different slopes in a power law of the volume fraction are identified.

FIG. 4.

The average interparticle gap h¯/a=R¯min/a2 as a function of volume fraction on a log–log scale. Four different zones associated with different slopes in a power law of the volume fraction are identified.

Close modal

The suspension microstructure is known to provide important information in complex flows, especially its implications to rheology [66–68]. An appropriate way to characterize the suspension microstructure, which is the spatial arrangement of particles, is to calculate the pair distribution function. This function provides information about the probability of finding a particle with respect to a probe one at the origin. Figure 5 shows the effect of volume fraction on this function. As seen in the figure, the maximum probability is first located near the particle poles at ϕ=10% and ϕ=20%, as previously observed in both suspension of spherical [11,12] and rodlike [69] particles. In the concentrated regime (ϕ>30%), the maximum shifts from the poles to equators due to the change of the dominant mechanism of pairing dynamics [11]. With a further increase in volume fraction up to ϕ=59%, this shifted maximum probability region seems to emanate and propagate over the particle surface toward the poles where the peak region finally turns to entirely cover the particle surface, forming seemingly microstructural isotropy. Specifically, starting at ϕ20%, the dominant mechanism and direction of particle parings change from attractive in the field direction to repulsive in the transverse direction due to strong ICEO flows in the lateral direction. This change in the dominant mechanism leads to pushing the particles to contact in the transverse direction, resulting in the appearance of the second high probability region at the equator [11]. The further increase in the repulsive interactions with volume fraction eventually causes nearly 95% of particle contacts to occur in the lateral directions, leading to a clear transition of the high probability region from the polar to the equator at ϕ35%. As a volume fraction is further increased, the increase in excluded volume interactions results in increasing particle contacts in all directions by which the high probability region starts to entirely cover particle surface for 45ϕ59%. Again, it should be noted that the local microstructure of the current system is primarily governed by ICEP as the pair distribution function of the suspensions undergoing only ICEP (not shown) is almost identical to that of Fig. 5. Finally, the function seems to show a crystal structure at random close packing (ϕ=64%), where the multiple high probability spots indicates the typical crystal structural information [70–72].

FIG. 5.

Pair distribution functions in suspensions undergoing DIP at six different volume fractions in coordinates (sgn(x)r,z), where r2=x2+y2. A probe particle is located at the origin (white areas are indeed an excluded volume).

FIG. 5.

Pair distribution functions in suspensions undergoing DIP at six different volume fractions in coordinates (sgn(x)r,z), where r2=x2+y2. A probe particle is located at the origin (white areas are indeed an excluded volume).

Close modal

Having determined the pair distribution functions, the suspension microstructure can be linked to the suspension entropy. Typically, the Shannon entropy S is calculated [73] and is given by

S(ϕ)=i=1MP(xi)log2P(xi),
(8)

where P is the pair distribution function at a volume fraction ϕ, xi is the ith location in P, and M is the total number of the locations in P. This entropy is also known to quantify the order of a suspension. Figure 6 shows the deviation of the suspension entropy normalized by the maximum entropy, 1S/Smax, as a function of volume fraction. The maximum entropy, Smax, results from an equiprobable pair distribution function, meaning that the probability of finding a particle is constant at all locations. It is found that the normalized suspension entropy (1S/Smax) decreases as ϕ0.75 from ϕ=0.25% to ϕ=20% at which the Shannon entropy S eventually reaches a maximum value, meaning that the suspension becomes the least ordered. The normalized entropy starts to increase and reach the maximum value at the random close packing, where the crystal structure is likely to be formed. At the volume fraction of the minimum normalized entropy (ϕ=20%), the secondary high probability region starts to appear at the equators in a pair distribution function (see Fig. 5). In addition, it is worth noting that at this volume fraction, the nature of pairing dynamics is changed from repulsive to attractive, where particle contacts along the transverse direction become predominant over the attractive contacts along the field direction [11].

FIG. 6.

The normalized deviation of the suspension entropy from the maximum entropy (1S(ϕ)/Smax) as a function of volume fraction on a log–log scale, where S(ϕ) is the entropy of suspension at a certain volume fraction ϕ and Smax is the maximum entropy resulted from the equiprobable pair distribution function.

FIG. 6.

The normalized deviation of the suspension entropy from the maximum entropy (1S(ϕ)/Smax) as a function of volume fraction on a log–log scale, where S(ϕ) is the entropy of suspension at a certain volume fraction ϕ and Smax is the maximum entropy resulted from the equiprobable pair distribution function.

Close modal

For defining the rheological properties of the suspension, the calculation of the bulk stress Σ is needed, where the angle bracket denotes an ensemble average over the particles. Here, we use the Batchelor calculation for the average stress tensor in a suspension of force-free particles [60]. The average of bulk stress is then determined by the average of the Cauchy stress tensor σ over a characteristic volume V,

Σ=1VVσ(x)dV=pfI+2ηE+Σp,
(9)

where pf=pf is the (constant) pressure in the fluid, I is the identity matrix, and 2ηE is the deviatoric contribution from an incompressible fluid containing rigid particles. These first two terms are the contributions of the fluid to the stress, comprising a Newtonian behavior [70]. The non-Newtonian behavior is captured by the particle contribution to the stress Σp, which is given by (in index notation)

Σijp=1Vα=1NA0[σikxjnkαη(uinjα+ujniα)]dA=1Vα=1NSijα+12Vα=1NϵijkLkα,
(10)

where A0 is the surface of particle α, ui is the velocity component in an ambient fluid, and niα is the component of the normal vector pointing outward from the particle surface to the fluid. Sα and Lαis the stresslet and torque on particle α, respectively. For torque-free rigid particles, in which case the present suspension is, the stresslet is the only contribution to the particle stress and is simply computed by

Sijα=A0[12(σikxj+σjkxi)nkα]dA,
(11)

where the position vector x can be taken from an arbitrary origin as long as the total hydrodynamic force on N particles vanishes. The particle stress tensor is then readily obtained by

Σijp=nSij,
(12)

where n=N/V is the number density of the suspension. Specifically, the particle stress is the number density times the average symmetric force dipole per particle.

For the suspension considered, which includes non-Brownian particles in a quiescent flow, the stresslet can be further decomposed into two stress-generating mechanisms, namely, ICEP effects and hydrodynamic interactions

Σijp=nSijICEP+nSijH.
(13)

The first term is associated with an induced nonlinear slip velocity over particle surface due purely to electric interactions, while the second term corresponds to the effect of hydrodynamic interactions between particles. Note that as mentioned above, the effect of electrostatic forces (DEP effect) is negligible compared with ICEP for the system studied; therefore, the contribution of DEP to the stress would not be included.

As the slip velocity on the particle surface drives the flow in a suspending fluid, the velocity field around a given particle due to the slip velocity can be found by solving the Stokes equations [13,34]. To leading order of O(R3), the contribution of the slip velocity to the particle stress SijICEP is now given by

SijICEP=916ϵE02(aR)3(AE^0+E^0A)+O(R4),
(14)

where A(I3R^R^)E^0, R^=R/a, and E^0=E0/E0. Subsequently, the velocity field generated by the slip velocity causes the motion of other particles, which can be obtained by Faxen’s law [52]. This particle motion corresponds to the first effect of hydrodynamic interactions. The contribution of this hydrodynamic interaction to the particle stress SijH is then captured by the disturbance velocity field due to the particle motion when placed in the velocity field generated by the slip velocity. In other words, the contribution to SijH is related to the hydrodynamic interactions captured by the reflection of a velocity disturbance from one particle to another. To leading order of O(R3) in the first effect of hydrodynamic interactions, the contribution to SijH is given by

SijH=4532ϵE02(aR)3[(I3R^R^)3(E^0R^)(E^0R^+R^E^0)3(E^0R^)2(I5R^R^)]+O(R4).
(15)

Given these two contributions to the particle stress, we turn to the dependence of the particle stress tensor on a volume fraction. Figure 7 shows the diagonal entities (i.e., normal stresses) in the particle stress tensor. It was found that the diagonal entities are several orders higher in magnitude than nondiagonal ones (i.e., shear stresses), suggesting that the rheological behavior of the current system is essentially governed by the normal stresses.

FIG. 7.

Dependence of the particle normal stresses, namely, ICEP SICEP and hydrodynamic SH contributions, on volume fraction in (a) the transverse (x,y) directions (note that Σxxp=nSxxH and Σyyp=nSyyH) and (b) the field (z) direction. The ICEP effect only contributes to the particle normal stress in the field direction. Due to symmetry with respect to the field direction, Σxxp=Σyyp.

FIG. 7.

Dependence of the particle normal stresses, namely, ICEP SICEP and hydrodynamic SH contributions, on volume fraction in (a) the transverse (x,y) directions (note that Σxxp=nSxxH and Σyyp=nSyyH) and (b) the field (z) direction. The ICEP effect only contributes to the particle normal stress in the field direction. Due to symmetry with respect to the field direction, Σxxp=Σyyp.

Close modal

Figure 7(a) shows the normal stresses in the transverse directions (x and y directions). To leading order, the ICEP effect appears to provide no contribution to the normal stresses in both x and y directions, essentially SxxICEP=SyyICEP=0. However, the hydrodynamic interactions contribute to the particle normal stress in the transverse directions, which are almost the same in both x and y directions, that is, SxxHSyyH. Furthermore, nontrivial behavior is observed in terms of sign and trend. The transverse hydrodynamic normal stresses are positive in dilute regimes and increase with volume fraction before reaching a local maximum at ϕ=25%. They start to decrease and become negative at ϕ32.5%, reaching a minimum at ϕ=40%. They start to increase again and become positive at ϕ=47.5%. Finally, they increase rapidly as approaching random close packing.

For the field direction, the particle normal stress is presented in Fig. 7(b). As opposed to the transverse directions, the ICEP effect contributes to the particle normal stress in the field direction. Nontrivial behaviors are also observed in the field direction. However, in this case, both ICEP and hydrodynamic contributions are all negative in dilute regimes and become more negative with volume fraction up to ϕ=25%. They then start to increase and become positive at ϕ32.5%, reaching a maximum at ϕ=45% and ϕ=40% for ICEP and hydrodynamic contributions, respectively. While the ICEP contribution stays positive for the rest of volume fraction with a slight decrease toward random close packing, the hydrodynamic contribution rapidly decreases and becomes negative again at ϕ=47.5%, reaching a minimum at random close packing. This nonmonotonic behavior of the particle normal stresses with volume fraction could imply distinctive rheological behaviors of the system studied.

Interestingly, for the hydrodynamic contributions, the normal stresses in the transverse directions, SxxH and SyyH, are exactly opposite to that in the field direction, SzzH.

This can be explained by the stress-generating response of the particle pairing originated due to DIP or ICEP. Specifically, for each pair of particles undergoing ICEP, SzzH=(SxxH+SyyH) always holds. Considering SyyHSxxH for the current system with an electric field applied in the z direction, then SzzH2SxxH. Note that the particle pairing dynamics, where particles tend to be attracted along the field direction and repulsive along the transverse direction, essentially results in exactly opposite flow fields around particles as the particles move in the transverse and field directions, resulting in the opposite signs for stresslets [11–13].

The sign of all the contributions in the particle normal stress seems to be correlated with the local microstructure seen in Fig. 5. In dilute regimes, the high probability region in the pair distribution function is located near the particle poles, which corresponds to a negative sign in SzzH<0 and SzzICEP<0, indicating that the particles are mostly paired up along the field direction (attractive pairings). Note that there is a positive sign in SxxH>0, but its effect is negligible on the microstructure. Their signs eventually change with the increase in volume fraction at ϕ32.5%, which corresponds to transition in local microstructure that the high probability regions now shift to the equators (repulsive pairings). With further increase in volume fraction, the second transition in particle normal stress associated with a change in the sign of hydrodynamic contributions is also confirmed by the local microstructure, where the high probability region in the equators propagates toward the poles. Therefore, it is suggested that the nonmonotonic variation of the particle stress system ties into the variation of the microstructure due to a change in the dominant mechanism and direction of particle parings [11].

It is worth mentioning that the contribution of the short-range repulsive force to the bulk stress can be negligible because it decays very quickly to zero within a very short distance ha [74]. To validate that the contribution of the short-range repulsive force to the bulk stress is negligible for the current system, we have calculated the hard-sphere interparticle stresslet SP=nxFP [75] in the range of volume fraction (not shown). It is seen that the maximum contribution of the short-range repulsive force to the particle normal stress is of O(104), which is almost two orders of magnitude smaller than that of ICEP and hydrodynamic contributions shown in Fig. 7.

The particle pressure represents the isotropic contribution of particles in bulk stress and is mechanically defined as the negative mean normal stress exerted by the particles in a viscous fluid, i.e., Π=1/3Σiip [61]. The particle pressure is also referred to as the nonequilibrium continuation of the osmotic pressure [76,77]. This particle pressure has been employed in many studies to characterize the rheological properties of the suspension of hard particles under the shear flow [70,76,78]. Although the whole suspension, which is the mixture of hard particles and a viscous fluid, is typically incompressible, the entire collection of the suspended particles is compressible from a macroscopic point of view [78]. With that said, the particle pressure has been quantified based upon the tendency of particle phase to expand [77]. As a way to evidencing the particle pressure, a U-shaped tube was sheared in a Couette device with a semipermeable membrane placed in the middle to separate a pure liquid on one side and a suspension on another (i.e., the particles are restricted to pass the membrane) [77]. The pure liquid has been observed to be sucked into the other side of a tube containing particles as the device is sheared. This observed suction indicates that the particle phase tends to expand, which means placing the suspending fluid in tension relatively, leading to the positive particle pressure [76]. In this regard, we investigate the sign of the particle pressure in the suspension studied.

Figure 8 shows the particle pressure as a function of volume fraction. To leading order of O(R3), it turns out that there is no hydrodynamic contribution to the particle pressure because SzzH=(SxxH+SyyH). The only contribution to the particle pressure results from the ICEP effect. At ϕ30%, the particle pressure is positive (Π>0) similar to that of the hard-sphere suspensions in a simple shear flow, which is always positive [76,77]. Again, it implies that the particle pressure places the suspending fluid in tension relative to the pure fluid state as the particles exhibit volumetric expansion regarding pf=Π. Interestingly, the particle pressure becomes negative at ϕ>30%, which indicates the change in its nature in a way that the suspending fluid is now relatively placed in compression as particles show volumetric contraction. Eventually, the particle pressure reaches its minimum and also maximum in magnitude at ϕ=45%. It is at this volume fraction that a nontrivial local maximum was observed in suspension dynamics, for instance, hydrodynamic diffusivity and velocity fluctuation [11]. It suggests a clear relationship between the particle pressure and the suspension dynamics in concentrated regimes, where the negative particle pressure might be the manifestation of the change in the predominance of particle pairings from mild attractive contacts to massive, strong repulsive contacts [11].

FIG. 8.

The particle pressure Π as a function of volume fraction. The nature of the particle pressure is changed at ϕ30% as it turns into negative from positive.

FIG. 8.

The particle pressure Π as a function of volume fraction. The nature of the particle pressure is changed at ϕ30% as it turns into negative from positive.

Close modal

It has been well-known that the hydrodynamics and rheology of suspensions are strongly affected by confinement [79–81]. We investigate the effects of confinement on the particle stress and pressure in the current system. To this end, we modify the boundary conditions in the z direction from periodic to wall-confined conditions as followed by our previous study [20], where only short-range interactions with the boundaries were captured. Specifically, to introduce the rigid boundaries in a domain, a short-range repulsive force is imposed between the particles and the boundaries by implementing the algorithm similar to the one used to prevent excessive particle overlaps. However, it should be noted that long-range interactions with the boundaries are likely to have an effect on particle dynamics, particularly in the direct vicinity of the boundaries, which will be included in future work. We introduce the confinement factor χ=2a/Lz to quantify the level of confinement, where Lz is the electrode spacing. The larger χ, the stronger the level of confinement. Typically, χ<0.05 corresponds to a weakly confined regime [82].

Figure 9 shows the dependence of the particle normal stresses on volume fraction for periodic suspensions and confined suspensions (χ=0.1, where Lz = 20). The confinement level χ=0.1 can be regarded as a moderately confined regime. The general shapes of the particle normal stress for a confined suspension are similar to ones for a periodic suspension. Notably, the confinement appears to cause the opposite effects on the transverse and field directions. Compared to the periodic suspension, the transverse normal stresses (Σxxp,Σyyp) are lower, while the field normal stress (Σzzp) is higher at all volume fractions considered. A difference of the particle stress between the periodic and confined suspensions gets larger as a volume fraction is increased and becomes constant in concentrated regimes.

FIG. 9.

The particle normal stresses in the transverse (x) and field (z) directions for periodic and confined suspensions. The confinement is placed in the z direction and the confinement level is χ=0.1, where χ=2a/Lz, Lz=20, and Lz is the electrode spacing. For both periodic and confined suspensions, ΣxxpΣyyp.

FIG. 9.

The particle normal stresses in the transverse (x) and field (z) directions for periodic and confined suspensions. The confinement is placed in the z direction and the confinement level is χ=0.1, where χ=2a/Lz, Lz=20, and Lz is the electrode spacing. For both periodic and confined suspensions, ΣxxpΣyyp.

Close modal

Figure 10 compares the particle pressure between periodic and confined suspensions (χ=0.1). A difference of the particle pressure between the periodic and confined suspensions keeps larger as a volume fraction is increased. The noticeable change is that the confinement makes the particle pressure decrease and become negative earlier starting at ϕ20%. It implies that the confinement essentially results in augmenting the volumetric particle contraction. The short-range repulsive interactions between particles and boundaries could be responsible for this change, for which detailed investigations will be included in future work.

FIG. 10.

The particle pressure Π for periodic and confined suspensions. The confinement level is χ=0.1, where the confinement factor χ=2a/Lz, Lz=20, and Lz is the electrode spacing.

FIG. 10.

The particle pressure Π for periodic and confined suspensions. The confinement level is χ=0.1, where the confinement factor χ=2a/Lz, Lz=20, and Lz is the electrode spacing.

Close modal

Lastly, the effects of the confinement factor χ on the particle stress and pressure are considered. Figure 11(a) shows the particle normal stresses of suspensions at ϕ=10% from weakly confined regime χ=0.05 to strongly confined regime χ=0.25 (note that χ=0.05 and 0.25 correspond to 20 and 4 particle diameters, respectively). As the confinement level gets stronger (i.e., increasing the confinement factor), the noticeable rise of the normal stress in the field direction is observed, while a relatively small decrease is observed for the transverse directions. Thus, it appears that the change in the particle normal stress in the z direction due to the increase in the level of confinement essentially governs the particle pressure reduction as no significant change is observed for the normal stresses in the x and y directions in the range of χ. More interestingly, Σxxp and Σyyp are almost the same up to χ=0.13 but become separated from each other and almost opposite in sign beyond χ=0.13. This separation might result from the instability induced by the strong confinement. Finally, the dependence of the particle pressure on the level of confinement is presented in Fig. 11(b). Five different volume fractions up to ϕ=25% are considered because the corresponding particle pressures are all positive for a periodic domain, as seen in Fig. 8. As the level of confinement increases, the particle pressure decreases with the confinement factor and then eventually becomes negative for all volume fractions. Moreover, the changeover at which the particle pressure turns into negative from positive arises earlier as a volume fraction increases. It is suggested that the strong confinement increases the volumetric contraction of the particle phase in a suspension.

FIG. 11.

The effect of the level of confinement on (a) the particle normal stresses at ϕ=10% and (b) the particle pressure at different volume fractions. The level of confinement is captured by the confinement factor χ=2a/Lz, where Lz is the electrode spacing.

FIG. 11.

The effect of the level of confinement on (a) the particle normal stresses at ϕ=10% and (b) the particle pressure at different volume fractions. The level of confinement is captured by the confinement factor χ=2a/Lz, where Lz is the electrode spacing.

Close modal

We have performed large-scale Stokesian dynamics simulations to study suspensions of ideally conductive spheres undergoing DIP, the combination of DEP and ICEP. In the current system, it is found that ICEP dominates DEP as a result of the slower decay of hydrodynamic interactions—a Stokes dipole for ICEP as O(R2) and a potential quadrupole for DEP as O(R4), where R is the separation distance between two particles. The suspension dynamics was investigated by the relaxation time, where it appears to relax very quickly with increasing volume fraction up to ϕ=35%, but the relaxation time increases up to ϕ47.5% and decreases again as approaching random close packing.

The particle stress tensor is computed to characterize the suspension rheology and shows that the normal stresses are predominant over the shear stresses. The particle normal stress is primarily a sum of two contributions, the ICEP effect, which only contributes to the normal stress in the field direction, and hydrodynamic interactions, which contribute to both transverse and field directions.

The nonmonotonic behaviors of both the particle normal stress (Σxxp,Σyyp,Σzzp) and the particle pressure (Π=1/3Σiip) with the volume fraction were observed and shown to tie into the suspension dynamics and microstructure. At ϕ<30%, the particle pressure Π>0, which is similar to a hard-sphere suspension in a simple shear flow, implying that the particle pressure places the suspending fluid in tension relatively. Interestingly, the nature of the particle pressure is changed at ϕ35% at which it turns negative, indicating that the suspending fluid is now placed in compression. The negative particle pressure eventually becomes maximum at ϕ=45%, where the hydrodynamic diffusivity and the velocity fluctuation reach a local maximum, as observed in our previous study [11]. It suggests that the negative particle pressure is connected to the observed nontrivial behaviors in the system.

Lastly, the wall confinement in the field direction is shown to strongly affect both the particle normal stress and particle pressure. It enhances the negative particle pressure drastically as a result of the augmentation of the volumetric particle contraction as a result of the short-range repulsive force between walls and particles. The confinement effect is further investigated by the confinement factor (χ=2a/Lz, where Lz is the electrode spacing). It is observed that the field normal stress is more affected by the confinement than the transverse ones, and the particle pressure turns negative much earlier with increasing confinement factor. For a confined suspension in microfluidic geometries, where the confinement is likely to introduce environmental heterogeneity, the conductive particles can display even more complex interactions with the walls, depending on the direction and frequency of the applied electric field, wall conductivity, and particle asymmetry [83–85]. Such interactions have yet to be fully explained and will be a subject of interesting future work.

Moving forward, the investigation of suspensions of ideally conductive particles in the shear flow is necessary for a complete study of the rheology of such suspension, which will be included in future work. Another interesting future work is to control the rheology of such suspension via surface treatments such as surface coating or ion surface absorption. It was shown that a thin dielectric coating over ideally conductive particles leads to the change in the suspension dynamics from ICEP-dominated to DEP-dominated, where the local aggregation of particles is observed [39]. Thus, the surface treatments can be considered as a tuning factor to control dynamics and rheology of such suspension, which may provide a useful avenue toward new engineering applications. Finally, these distinct changing nature of rheological properties could suggest that different flow control techniques could be used to delay or promote the changing nature in such suspension systems.

The authors gratefully acknowledge financial support from the Collaboration Initiative and Interdisciplinary Research Grants at the University of Nebraska and, in part, from the National Science Foundation through Grant No. CBET-1936065 (Particulate and Multiphase Processes program). This work was completed utilizing the Holland Computing Center of the University of Nebraska, which receives support from the Nebraska Research Initiative.

For large-scale dynamics driven by DEP and ICEP, we wish to calculate sums of form

us(xα)=β=1NSP(xβxα):z^z^,
(A1)
ut(xα)=β=1NTP(xβxα):z^z^,
(A2)

where α=1,,N and Sp and Tp denote the periodic versions of Green’s function of Stoke dipole and a potential quadrupole. By making use of the known Ewald summation formula [53,86], Eqs. (A1) and (A2) can be recast into the following Ewald summations:

us(xα)=pβ=1NAs(ξ,xβxα+p):z^z^+k0e2πikxαS(k)Bs(ξ,k):z^z^,
(A3)
ut(xα)=pβ=1NAt(ξ,xβxα+p):z^z^+k0e2πikxαS(k)Bt(ξ,k):z^z^,
(A4)

where ξ, called the Ewald coefficient, which determines the relative importance of the real and Fourier sums. This coefficient is user-defined and is chosen to minimize the overall cost of the algorithm. The first sums (real sums) in Eqs. (A3) and (A4) are over all particle positions xβ and their periodic images (which are denoted by the lattice vectors p), and the second sums (Fourier sums) are over wave vector k. The structure factor S(k) in Eqs. (A3) and (A4) is obtained for the suspension as follows:

S(k)=β=1Ne2πikxβ.
(A5)

The convolution kernels As, At, Bs, and Bt are third-order tensors and can be calculated analytically (see Park and Saintillan [12]). These tenors decay exponentially, which consequently results in exponential convergence of the sums in Eqs. (A3) and (A4). The details of evaluating these tenors in the SPME algorithm are provided in Saintillan et al. [86]. An O(N) cost for the evaluation of the real sums at all particle positions is obtainable by choosing Ewald coefficient ξ so as to exclude all the terms beyond a fixed cutoff distance rc from the reals sums, which allows truncation of these sums after a finite number of terms independent of the system size. For the Fourier sums (second sums), the particle distribution is transformed to Fourier space using the fast Fourier transform algorithm, after assigning to Cartesian gird by B-spline interpolation [87]. It yields structure factor S(k), which is then multiplied by the convolution kernels Bs and Bt. The inverse Fourier transform is applied, and the value of Fourier sums is determined at the particle locations by interpolation. The computation cost of Fourier sums is limited by the fast Fourier transform algorithm, scaling as O(KlogK) with the number K of grid points (or Fourier modes). This number is typically chosen to be proportional to the number of particles N; therefore, the overall cost for the velocities evaluation is O(NlogN).

1.
Bazant
,
M. Z.
, and
T. M.
Squires
, “
Induced-charge electrokinetic phenomena
,”
Curr. Opin. Colloid Interface Sci.
15
,
203
213
(
2010
).
2.
Sheng
,
P.
, and
W.
Wen
, “
Electrorheological fluids: Mechanisms, dynamics, and microfluidics applications
,”
Annu. Rev. Fluid Mech.
44
,
143
174
(
2012
).
3.
Nikonenko
,
V. V.
,
A. V.
Kovalenko
,
M. K.
Urtenov
,
N. D.
Pismenskaya
,
J.
Han
,
P.
Sistat
, and
G.
Pourcelly
, “
Desalination at overlimiting currents: State-of-the-art and perspectives
,”
Desalination
342
,
85
106
(
2014
).
4.
Li
,
J.
,
X.
Liang
,
F.
Liou
, and
J.
Park
, “
Macro-/micro-controlled 3D lithium-ion batteries via additive manufacturing and electric field processing
,”
Sci. Rep.
8
,
1
11
(
2018
).
5.
Wang
,
X.-B.
,
J.
Yang
,
Y.
Huang
,
J.
Vykoukal
,
F. F.
Becker
, and
P. R.
Gascoyne
, “
Cell separation by dielectrophoretic field-flow-fractionation
,”
Anal. Chem.
72
,
832
839
(
2000
).
6.
Presser
,
V.
,
C. R.
Dennison
,
J.
Campos
,
K. W.
Knehr
,
E. C.
Kumbur
, and
Y.
Gogotsi
, “
The electrochemical flow capacitor: A new concept for rapid energy storage and recovery
,”
Adv. Energy Mater.
2
,
895
902
(
2012
).
7.
Soloveichik
,
G. L.
, “
Flow batteries: Current status and trends
,”
Chem. Rev.
115
,
11533
11558
(
2015
).
8.
Tan
,
H. W.
,
J.
An
,
C. K.
Chua
, and
T.
Tran
, “
Metallic nanoparticle inks for 3D printing of electronics
,”
Adv. Electronic Mater.
5
,
1800831
(
2019
).
9.
Bauer
,
W.
, and
D.
Nötzel
, “
Rheological properties and stability of NMP based cathode slurries for lithium ion batteries
,”
Ceram. Int.
40
,
4591
4598
(
2014
).
10.
Akuzum
,
B.
,
K.
Maleski
,
B.
Anasori
,
P.
Lelyukh
,
N. J.
Alvarez
,
E. C.
Kumbur
, and
Y.
Gogotsi
, “
Rheological characteristics of 2D titanium carbide (MXene) dispersions: A guide for processing MXenes
,”
ACS Nano.
12
,
2685
2694
(
2018
).
11.
Mirfendereski
,
S.
, and
J. S.
Park
, “
Dipolophoresis in concentrated suspensions of ideally polarizable spheres
,”
J. Fluid Mech.
875
,
R3
(
2019
).
12.
Park
,
J. S.
, and
D.
Saintillan
, “
Dipolophoresis in large-scale suspensions of ideally polarizable spheres
,”
J. Fluid Mech.
662
,
66
90
(
2010
).
13.
Saintillan
,
D.
, “
Nonlinear interactions in electrophoresis of ideally polarizable particles
,”
Phys. Fluids
20
,
067104
(
2008
).
14.
Hao
,
T.
, “
Electrorheological fluids
,”
Adv. Mater.
13
,
1847
1857
(
2001
).
15.
Pohl
,
H. A.
,
Dielectrophoresis: The Behavior of Neutral Matter in Nonuniform Electric Fields
, Cambridge Monographs on Physics (
Cambridge University
,
Cambridge, UK
,
1978
).
16.
von Pfeil
,
K.
,
M. D.
Graham
,
D. J.
Klingenberg
, and
J. F.
Morris
, “
Pattern formation in flowing electrorheological fluids
,”
Phys. Rev. Lett.
88
,
188301
(
2002
).
17.
Halsey
,
T. C.
, “
Electrorheological fluids
,”
Science
258
,
761
766
(
1992
).
18.
Halsey
,
T. C.
, and
W.
Toor
, “
Structure of electrorheological fluids
,”
Phys. Rev. Lett.
65
,
2820
2823
(
1990
).
19.
Fossum
,
J.
,
Y.
Méheust
,
K.
Parmar
,
K.
Knudsen
,
K.
Måløy
, and
D. M.
Fonseca
, “
Intercalation-enhanced electric polarization and chain formation of nano-layered particles
,”
Europhys. Lett.
74
,
438
444
(
2006
).
20.
Park
,
J. S.
, and
D.
Saintillan
, “
Electric-field-induced ordering and pattern formation in colloidal suspensions
,”
Phys. Rev. E
83
,
041409
(
2011
).
21.
Winslow
,
W. M.
, “Method and means for translating electrical impulses into mechanical force,” U.S. Patent 2,417,850 (1947).
22.
Parthasarathy
,
M.
, and
D. J.
Klingenberg
, “
Electrorheology: Mechanisms and models
,”
Mater. Sci. Eng. R. Rep.
17
,
57
103
(
1996
).
23.
Bonnecaze
,
R. T.
, and
J.
Brady
, “
Dynamic simulation of an electrorheological fluid
,”
J. Chem. Phys.
96
,
2183
2202
(
1992
).
24.
De Vicente
,
J.
,
D. J.
Klingenberg
, and
R.
Hidalgo-Alvarez
, “
Magnetorheological fluids: A review
,”
Soft Matter
7
,
3701
3710
(
2011
).
25.
Bossis
,
G.
,
E.
Lemaire
,
O.
Volkova
, and
H.
Clercx
, “
Yield stress in magnetorheological and electrorheological fluids: A comparison between microscopic and macroscopic structural models
,”
J. Rheol.
41
,
687
704
(
1997
).
26.
Klingenberg
,
D. J.
, “
Magnetorheology: Applications and challenges
,”
AIChE J.
47
,
246
249
(
2001
).
27.
Martin
,
J. E.
,
R. A.
Anderson
, and
C. P.
Tigges
, “
Simulation of the athermal coarsening of composites structured by a uniaxial field
,”
J. Chem. Phys.
108
,
3765
3787
(
1998
).
28.
Adriani
,
P. M.
, and
A. P.
Gast
, “
A microscopic model of electrorheology
,”
Phys. Fluids
31
,
2757
2768
(
1988
).
29.
Sim
,
H. G.
,
K. H.
Ahn
, and
S. J.
Lee
, “
Three-dimensional dynamics simulation of electrorheological fluids under large amplitude oscillatory shear flow
,”
J. Rheol.
47
,
879
895
(
2003
).
30.
Bonnecaze
,
R. T.
, and
J.
Brady
, “
Yield stresses in electrorheological fluids
,”
J. Rheol.
36
,
73
115
(
1992
).
31.
Domínguez-García
,
P.
,
S.
Melle
,
O. G.
Calderón
, and
M.
Rubio
, “
Doublet dynamics of magnetizable particles under frequency modulated rotating fields
,”
Colloids Surf. A Physicochem. Eng. Aspects
270
,
270
276
(
2005
).
32.
Klingenberg
,
D. J.
,
J. C.
Ulicny
, and
M. A.
Golden
, “
Mason numbers for magnetorheology
,”
J. Rheol.
51
,
883
893
(
2007
).
33.
Anderson
,
R.
, Effects of finite conductivity in electrorheological fluids, Technical Report, Institution Sandia National, Albuquerque, NM, 1991.
34.
Squires
,
T. M.
, and
M. Z.
Bazant
, “
Induced-charge electro-osmosis
,”
J. Fluid Mech.
509
,
217
252
(
2004
).
35.
Squires
,
T. M.
, and
M. Z.
Bazant
, “
Breaking symmetries in induced-charge electro-osmosis and electrophoresis
,”
J. Fluid Mech.
560
,
65
101
(
2006
).
36.
Gamayunov
,
N.
,
V.
Murtsovkin
, and
A.
Dukhin
, “
Pair interaction of particles in electric field. 1. Features of hydrodynamic interaction of polarized particles
,”
Colloid J. USSR
48
,
233
239
(
1986
).
37.
Dukhin
,
A.
, and
V.
Murtsovkin
, “
Pair interaction of particles in electric field. 2. Influence of polarization of double layer of dielectric particles on their hydrodynamic interaction in a stationary electric field
,”
Colloid J. USSR
48
,
240
247
(
1986
).
38.
Murtsovkin
,
V.
, “
Nonlinear flows near polarized disperse particles
,”
Colloid J. Russ. Acad. Sci.
58
,
341
349
(
1996
).
39.
Park
,
J. S.
, and
D.
Saintillan
, “
From diffusive motion to local aggregation: Effect of surface contamination in dipolophoresis
,”
Soft Matter
7
,
10720
10727
(
2011
).
40.
Lauga
,
E.
, and
T. R.
Powers
, “
The hydrodynamics of swimming microorganisms
,”
Rep. Prog. Phys.
72
,
096601
(
2009
).
41.
Saintillan
,
D.
, “
Rheology of active fluids
,”
Annu. Rev. Fluid Mech.
50
,
563
592
(
2018
).
42.
Saintillan
,
D.
, “
Extensional rheology of active suspensions
,”
Phys. Rev. E
81
,
056307
(
2010
).
43.
Ishikawa
,
T.
,
M.
Simmonds
, and
T. J.
Pedley
, “
Hydrodynamic interaction of two swimming model micro-organisms
,”
J. Fluid Mech.
568
,
119
160
(
2006
).
44.
Evans
,
A. A.
,
T.
Ishikawa
,
T.
Yamaguchi
, and
E.
Lauga
, “
Orientational order in concentrated suspensions of spherical microswimmers
,”
Phys. Fluids
23
,
111702
(
2011
).
45.
Hatwalne
,
Y.
,
S.
Ramaswamy
,
M.
Rao
, and
R. A.
Simha
, “
Rheology of active-particle suspensions
,”
Phys. Rev. Lett.
92
,
118101
(
2004
).
46.
Haines
,
B. M.
,
I. S.
Aranson
,
L.
Berlyand
, and
D. A.
Karpeev
, “
Effective viscosity of dilute bacterial suspensions: A two-dimensional model
,”
Phys. Biol.
5
,
046003
(
2008
).
47.
Haines
,
B. M.
,
A.
Sokolov
,
I. S.
Aranson
,
L.
Berlyand
, and
D. A.
Karpeev
, “
Three-dimensional model for the effective viscosity of bacterial suspensions
,”
Phys. Rev. E
80
,
041922
(
2009
).
48.
McDonnell
,
A. G.
,
T. C.
Gopesh
,
J.
Lo
,
M.
O’Bryan
,
L. Y.
Yeo
,
J. R.
Friend
, and
R.
Prabhakar
, “
Motility induced changes in viscosity of suspensions of swimming microbes in extensional flows
,”
Soft Matter
11
,
4658
4668
(
2015
).
49.
Shilov
,
V.
, and
T.
Simonova
, “
Polarization of electric double-layer of disperse particles and dipolophoresis in a steady (DC) field
,”
Colloid J. USSR
43
,
90
96
(
1981
).
50.
Bazant
,
M. Z.
,
M. S.
Kilic
,
B. D.
Storey
, and
A.
Ajdari
, “
Towards an understanding of induced-charge electrokinetics at large applied voltages in concentrated solutions
,”
Adv. Colloid Interface Sci.
152
,
48
88
(
2009
).
51.
Jones
,
T. B.
, and
T. B.
Jones
,
Electromechanics of Particles
(
Cambridge University
,
Cambridge
,
2005
).
52.
Kim
,
S.
, and
S. J.
Karrila
,
Microhydrodynamics: Principles and Selected Applications
(
Courier Corporation
,
Mineola, NY
,
2013
).
53.
Hasimoto
,
H.
, “
On the periodic fundamental solutions of the Stokes equations and their application to viscous flow past a cubic array of spheres
,”
J. Fluid Mech.
5
,
317
328
(
1959
).
54.
Melrose
,
J.
, and
D.
Heyes
, “
Simulations of electrorheological and particle mixture suspensions: Agglomerate and layer structures
,”
J. Chem. Phys.
98
,
5873
5886
(
1993
).
55.
Sherman
,
Z. M.
,
D.
Ghosh
, and
J. W.
Swan
, “
Field-directed self-assembly of mutually polarizable nanoparticles
,”
Langmuir
34
,
7117
7134
(
2018
).
56.
Rintoul
,
M.
, and
S.
Torquato
, “
Computer simulations of dense hard-sphere systems
,”
J. Chem. Phys.
105
,
9258
9265
(
1996
).
57.
Clarke
,
A.
, and
J.
Wiley
, “
Numerical simulation of the dense random packing of a binary mixture of hard spheres: Amorphous metals
,”
Phys. Rev. B
35
,
7350
7356
(
1987
).
58.
Stillinger Jr
,
F.
,
E.
DiMarzio
, and
R.
Kornegay
, “
Systematic approach to explanation of the rigid disk phase transition
,”
J. Chem. Phys.
40
,
1564
1576
(
1964
).
59.
Jodrey
,
W.
, and
E.
Tory
, “
Computer simulation of close random packing of equal spheres
,”
Phys. Rev. A
32
,
2347
2351
(
1985
).
60.
Batchelor
,
G.
, “
The stress system in a suspension of force-free particles
,”
J. Fluid Mech.
41
,
545
570
(
1970
).
61.
Jeffrey
,
D.
,
J.
Morris
, and
J. F.
Brady
, “
The pressure moments for two rigid spheres in low-Reynolds-number flow
,”
Phys. Fluids A Fluid Dyn.
5
,
2317
2325
(
1993
).
62.
Hao
,
T.
,
Electrorheological Fluids: The Non-aqueous Suspensions
(
Elsevier
,
Amsterdam
,
2011
).
63.
Michot
,
L. J.
,
I.
Bihannic
,
S.
Maddi
,
S. S.
Funari
,
C.
Baravian
,
P.
Levitz
, and
P.
Davidson
, “
Liquid–crystalline aqueous clay suspensions
,”
Proc. Natl. Acad. Sci. U.S.A.
103
,
16101
16104
(
2006
).
64.
Jiang
,
W.
,
C.-h.
Liu
,
Z.-g.
Wang
,
L.-j.
An
,
H.-j.
Liang
,
B.-z.
Jiang
,
X.-h.
Wang
, and
H.-x.
Zhang
, “
Brittle-tough transition in PP/EPDM blends: Effects of interparticle distance and temperature
,”
Polymer
39
,
3285
3288
(
1998
).
65.
Borggreve
,
R.
,
R.
Gaymans
,
J.
Schuijer
, and
J. I.
Housz
, “
Brittle-tough transition in nylon-rubber blends: Effect of rubber concentration and particle size
,”
Polymer
28
,
1489
1496
(
1987
).
66.
Stickel
,
J. J.
, and
R. L.
Powell
, “
Fluid mechanics and rheology of dense suspensions
,”
Annu. Rev. Fluid Mech.
37
,
129
149
(
2005
).
67.
Morris
,
J. F.
, “
A review of microstructure in concentrated suspensions and its implications for rheology and bulk flow
,”
Rheol. Acta
48
,
909
923
(
2009
).
68.
Butler
,
J. E.
, and
B.
Snook
, “
Microstructural dynamics and rheology of suspensions of rigid fibers
,”
Annu. Rev. Fluid Mech.
50
,
299
318
(
2018
).
69.
Saintillan
,
D.
,
E.
Darve
, and
E. S.
Shaqfeh
, “
Hydrodynamic interactions in the induced-charge electrophoresis of colloidal rod dispersions
,”
J. Fluid Mech.
563
,
223
259
(
2006
).
70.
Sierou
,
A.
, and
J.
Brady
, “
Rheology and microstructure in concentrated noncolloidal suspensions
,”
J. Rheol.
46
,
1031
1056
(
2002
).
71.
Foss
,
D. R.
, and
J. F.
Brady
, “
Structure, diffusion and rheology of Brownian suspensions by Stokesian dynamics simulation
,”
J. Fluid Mech.
407
,
167
200
(
2000
).
72.
Gasser
,
U.
, and
A.
Fernandez-Nieves
, “
Crystal structure of highly concentrated, ionic microgel suspensions studied by small-angle x-ray scattering
,”
Phys. Rev. E
81
,
052401
(
2010
).
73.
Shannon
,
C. E.
, “
A mathematical theory of communication
,”
Bell Syst. Tech. J.
27
,
379
423
(
1948
).
74.
Nazockdast
,
E.
, and
J. F.
Morris
, “
Microstructural theory and the rheology of concentrated colloidal suspensions
,”
J. Fluid Mech.
713
,
420
452
(
2012
).
75.
Singh
,
A.
, and
P. R.
Nott
, “
Normal stresses and microstructure in bounded sheared suspensions via Stokesian dynamics simulations
,”
J. Fluid Mech.
412
,
279
301
(
2000
).
76.
Yurkovetsky
,
Y.
, and
J. F.
Morris
, “
Particle pressure in sheared Brownian suspensions
,”
J. Rheol.
52
,
141
164
(
2008
).
77.
Deboeuf
,
A.
,
G.
Gauthier
,
J.
Martin
,
Y.
Yurkovetsky
, and
J. F.
Morris
, “
Particle pressure in a sheared suspension: A bridge from osmosis to granular dilatancy
,”
Phys. Rev. Lett.
102
,
108301
(
2009
).
78.
Guazzelli
,
É.
, and
O.
Pouliquen
, “
Rheology of dense granular suspensions
,”
J. Fluid Mech.
852
,
P1
(
2018
).
79.
Swan
,
J. W.
, and
J. F.
Brady
, “
The hydrodynamics of confined dispersions
,”
J. Fluid Mech.
687
,
254
299
(
2011
).
80.
Jeanneret
,
R.
,
D. O.
Pushkin
, and
M.
Polin
, “
Confinement enhances the diversity of microbial flow fields
,”
Phys. Rev. Lett.
123
,
248102
(
2019
).
81.
Ramaswamy
,
M.
,
N. Y.
Lin
,
B. D.
Leahy
,
C.
Ness
,
A. M.
Fiore
,
J. W.
Swan
, and
I.
Cohen
, “
How confinement-induced structures alter the contribution of hydrodynamic and short-ranged repulsion forces to the viscosity of colloidal suspensions
,”
Phys. Rev. X
7
,
041005
(
2017
).
82.
Graham
,
M. D.
, “
Fluid dynamics of dissolved polymer molecules in confined geometries
,”
Annu. Rev. Fluid Mech.
43
,
273
298
(
2011
).
83.
Kilic
,
M. S.
, and
M. Z.
Bazant
, “
Induced-charge electrophoresis near a wall
,”
Electrophoresis
32
,
614
628
(
2011
).
84.
Gangwal
,
S.
,
O. J.
Cayre
,
M. Z.
Bazant
, and
O. D.
Velev
, “
Induced-charge electrophoresis of metallodielectric particles
,”
Phys. Rev. Lett.
100
,
058302
(
2008
).
85.
Wu
,
Z.
,
Y.
Gao
, and
D.
Li
, “
Electrophoretic motion of ideally polarizable particles in a microchannel
,”
Electrophoresis
30
,
773
781
(
2009
).
86.
Saintillan
,
D.
,
E.
Darve
, and
E. S.
Shaqfeh
, “
A smooth particle-mesh Ewald algorithm for Stokes suspension simulations: The sedimentation of fibers
,”
Phys. Fluids
17
,
033301
(
2005
).
87.
Essmann
,
U.
,
L.
Perera
,
M. L.
Berkowitz
,
T.
Darden
,
H.
Lee
, and
L. G.
Pedersen
, “
A smooth particle mesh Ewald method
,”
J. Chem. Phys.
103
,
8577
8593
(
1995
).