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.

## I. INTRODUCTION

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 $\mu $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(R\u22122)$ than the later as $O(R\u22124)$ [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.

## II. GOVERNING EQUATION AND SIMULATION METHOD

We consider a suspension of $N$ identical neutrally buoyant ideally polarizable spheres of radius $a$ in a viscous electrolyte with the permittivity $\epsilon $ and viscosity $\eta $. 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 $\tau c=\lambda Da/D$, where $\lambda 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 ($a\u223c10$ $\mu $m, $\lambda D\u223c10$ nm, $D\u223c10\u22125$ $cm2$ s$\u22121$), the charging time $\tau c\u223c10\u22124$ s, which is much smaller than the diffusion time across the particle $\tau a=a2/D\u223c10\u22121$ 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 ($\lambda D/a\u224810\u22123$) and the minimum distance between particles in a suspension is larger than $2\lambda 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\u02d9\alpha $ of a given particle $\alpha $ in a suspension can be expressed by

where $R\alpha \beta =x\beta \u2212x\alpha $ is the separation vector between particle $\alpha $ and particle $\beta $, 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 $\lambda =2a/|R|$. Specifically, for both DEP and ICEP, $M$ is calculated as follows:

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

where the two tensors $S$ and $T=\u22072S$ 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:

These far-field tensors are asymptotically valid to order $O(\lambda 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 ($\lambda \u226a1$) 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\alpha \beta |<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\alpha \beta |\u22482.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 $\alpha $ and $\beta $ 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\alpha \beta \u22124)$. 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(Nlog\u2061N)$ 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\u02d9\alpha $ ($\alpha =1,\u2026,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 $\Delta 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 $(\u223c2.005a)$, to contact. The form of the repulsive potential for excluded volume ($EV$) interactions is

where $k$ is the time step-dependent prefactor, which can be expressed as $k=3\pi \eta a/\Delta 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=\epsilon aE02/\eta $.

## III. RESULTS AND DISCUSSION

We performed large-scale simulations in a cubic periodic domain ($Lx\xd7Ly\xd7Lz=203$) for a range of volume fraction $\varphi $ up to almost random close packing ($\varphi \u224864%$). 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 $\varphi =35$ before reaching a local maximum at $\varphi \u224845%$ 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 $\Sigma 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].

### A. Suspension dynamics: Relaxation time and interparticle distance

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.

Figure 3 shows the relaxation time $\tau $ as a function of volume fraction $\varphi $ on a log-log scale. In the dilute regime $(\varphi <5%)$, the relaxation time tends to decrease with volume fraction and scale as $\tau \u223c\varphi \u22120.5$. Subsequently, in the semidilute regime ($\varphi =5\u2212\u221235%$), the relaxation time decreases faster than for the dilute regime, approximately scaling as $\tau \u223c\varphi \u22121$. 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 $\varphi \u224835%$ up to $\varphi \u224847.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].

To further distinguish different behaviors in different ranges of volume fraction, the mean interparticle gap $h\xaf=R\xafmin\u22122a$ 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\xaf\u223c\varphi \u22121/3$, which has also been observed in many dilute suspensions [62–65]. The exponent then becomes larger at $\varphi \u22485%$, providing the second slope of $\u22120.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).

### B. Suspension microstructure and entropy

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 $\varphi =10%$ and $\varphi =20%$, as previously observed in both suspension of spherical [11,12] and rodlike [69] particles. In the concentrated regime ($\varphi >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 $\varphi =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 $\varphi \u224820%$, 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 $\varphi \u224835%$. 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\u2a7d\varphi \u2a7d59%$. 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 ($\varphi =64%$), where the multiple high probability spots indicates the typical crystal structural information [70–72].

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

where $P$ is the pair distribution function at a volume fraction $\varphi $, $xi$ is the $i$th 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, $1\u2212S/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 ($1\u2212S/Smax$) decreases as $\u223c\varphi \u22120.75$ from $\varphi =0.25%$ to $\varphi =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 ($\varphi =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].

### C. Particle stress

For defining the rheological properties of the suspension, the calculation of the bulk stress $\u27e8\Sigma \u27e9$ 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 $\sigma $ over a characteristic volume $V$,

where $pf=\u27e8p\u27e9f$ is the (constant) pressure in the fluid, $I$ is the identity matrix, and $2\eta \u27e8E\u221e\u27e9$ 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 $\Sigma p$, which is given by (in index notation)

where $A0$ is the surface of particle $\alpha $, $ui$ is the velocity component in an ambient fluid, and $ni\alpha $ is the component of the normal vector pointing outward from the particle surface to the fluid. $S\alpha $ and $L\alpha $is the stresslet and torque on particle $\alpha $, 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

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

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

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(R\u22123)$, the contribution of the slip velocity to the particle stress $SijICEP$ is now given by

where $A\u2261(I\u22123R^R^)\u22c5E^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(R\u22123)$ in the first effect of hydrodynamic interactions, the contribution to $SijH$ is given by

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.

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, $SxxH\u2248SyyH$. 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 $\varphi =25%$. They start to decrease and become negative at $\varphi \u224832.5%$, reaching a minimum at $\varphi =40%$. They start to increase again and become positive at $\varphi =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 $\varphi =25%$. They then start to increase and become positive at $\varphi \u224832.5%$, reaching a maximum at $\varphi =45%$ and $\varphi =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 $\varphi =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, $\u27e8SxxH\u27e9$ and $\u27e8SyyH\u27e9$, are exactly opposite to that in the field direction, $\u27e8SzzH\u27e9$.

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=\u2212(SxxH+SyyH)$ always holds. Considering $\u27e8SyyH\u27e9\u2248\u27e8SxxH\u27e9$ for the current system with an electric field applied in the $z$ direction, then $\u27e8SzzH\u27e9\u2248\u22122\u27e8SxxH\u27e9$. 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 $\u27e8SzzH\u27e9<0$ and $\u27e8SzzICEP\u27e9<0$, indicating that the particles are mostly paired up along the field direction (attractive pairings). Note that there is a positive sign in $\u27e8SxxH\u27e9>0$, but its effect is negligible on the microstructure. Their signs eventually change with the increase in volume fraction at $\varphi \u224832.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 $h\u226aa$ [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 $\u27e8SP\u27e9=\u2212n\u27e8xFP\u27e9$ [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(10\u22124)$, which is almost two orders of magnitude smaller than that of ICEP and hydrodynamic contributions shown in Fig. 7.

### D. Particle pressure

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., $\Pi =\u22121/3\Sigma 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(R\u22123)$, it turns out that there is no hydrodynamic contribution to the particle pressure because $\u27e8SzzH\u27e9=\u2212(\u27e8SxxH\u27e9+\u27e8SyyH\u27e9)$. The only contribution to the particle pressure results from the ICEP effect. At $\varphi \u2a7d30%$, the particle pressure is positive ($\Pi >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=\u2212\Pi $. Interestingly, the particle pressure becomes negative at $\varphi >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 $\varphi =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].

### E. Confinement effect

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 $\chi =2a/Lz$ to quantify the level of confinement, where $Lz$ is the electrode spacing. The larger $\chi $, the stronger the level of confinement. Typically, $\chi <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 ($\chi =0.1$, where $Lz$ = 20). The confinement level $\chi =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 ($\Sigma xxp,\Sigma yyp$) are lower, while the field normal stress ($\Sigma 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.

Figure 10 compares the particle pressure between periodic and confined suspensions ($\chi =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 $\varphi \u224820%$. 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.

Lastly, the effects of the confinement factor $\chi $ on the particle stress and pressure are considered. Figure 11(a) shows the particle normal stresses of suspensions at $\varphi =10%$ from weakly confined regime $\chi =0.05$ to strongly confined regime $\chi =0.25$ (note that $\chi =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 $\chi $. More interestingly, $\Sigma xxp$ and $\Sigma yyp$ are almost the same up to $\chi =0.13$ but become separated from each other and almost opposite in sign beyond $\chi =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 $\varphi =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.

## IV. CONCLUDING REMARKS

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(R\u22122)$ and a potential quadrupole for DEP as $O(R\u22124)$, 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 $\varphi =35%$, but the relaxation time increases up to $\varphi \u224847.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 ($\Sigma xxp,\Sigma yyp,\Sigma zzp$) and the particle pressure ($\Pi =\u22121/3\Sigma iip$) with the volume fraction were observed and shown to tie into the suspension dynamics and microstructure. At $\varphi <30%$, the particle pressure $\Pi >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 $\varphi \u224835%$ at which it turns negative, indicating that the suspending fluid is now placed in compression. The negative particle pressure eventually becomes maximum at $\varphi =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 ($\chi =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.

## ACKNOWLEDGMENTS

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.

### APPENDIX: SPME ALGORITHM

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

where $\alpha =1,\u2026,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:

where $\xi $, 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\beta $ 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:

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 $\xi $ 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(Klog\u2061K)$ 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(Nlog\u2061N)$.