Effects of electrostatic interaction on clustering and collision of bidispersed inertial particles in homogeneous and isotropic turbulence

Abstract In sandstorms and thunderclouds, turbulence-induced collisions between solid particles and ice crystals lead to inevitable triboelectrification. The charge segregation is usually size dependent, with small particles charged negatively and large particles charged positively. In this work, we perform numerical simulations to study the influence of charge segregation on the dynamics of bidispersed inertial particles in turbulence. Direct numerical simulations of homogeneous isotropic turbulence are performed with the Taylor Reynolds number ${Re}_{\lambda }=147.5$, while particles are subjected to both electrostatic interactions and fluid drag, with Stokes numbers of 1 and 10 for small and large particles, respectively. Coulomb repulsion/attraction is shown to effectively inhibit/enhance particle clustering within a short range. Besides, the mean relative velocity between same-size particles is found to rise as the particle charge increases because of the exclusion of low-velocity pairs, while the relative velocity between different-size particles is almost unaffected, emphasizing the dominant roles of differential inertia. The mean Coulomb-turbulence parameter, ${Ct}_0$, is then defined to characterize the competition between the Coulomb potential energy and the mean relative kinetic energy. In addition, a model is proposed to quantify the rate at which charged particles approach each other and to capture the transition of the particle relative motion from the turbulence-dominated regime to the electrostatic-dominated regime. Finally, the probability distribution function of the approach rate between particle pairs is examined, and its dependence on the Coulomb force is further discussed using the extended Coulomb-turbulence parameter.


Introduction
Particle-laden turbulent flows are widespread in both nature and industry.One of the most important features of turbulence is its ability to effectively transport suspended particles and create highly inhomogeneous particle distributions, which eventually lead to frequent particle collisions or bubble/droplet coalescence (Falkovich et al. 2001;Toschi & Bodenschatz 2009;Pumir & Wilkinson 2016).Typical natural processes can be found as rain formation (Shaw 2003;Grabowski & Wang 2013), sandstorms (Zhang & Zhou 2020), and the formation of marine snow in the ocean (Arguedas-Leiva et al. 2022), while in industrial applications turbulence-induced collisions are also shown essential in dusty particle removal (Jaworek et al. 2018) and flocculation (Zhao et al. 2021).
It has been widely shown that the spatial distribution of inertial particles is highly nonuniform and intermittent, which is known as the preferential concentration (Maxey 1987;Squires & Eaton 1991).Preferential concentration could significantly enhance local particle number density and facilitate interparticle collisions (Reade & Collins 2000).At the dissipative range, clustering is driven by Kolmogorov-scale eddies that eject heavy particles out of high-vorticity regions while entrap light bubbles (Balkovsky et al. 2001;Calzavarini et al. 2008;Mathai et al. 2016).Moreover, later studies show that in addition to clustering at the dissipative scale, the coexistence of multiscale vortices in turbulence also affect the dynamics of heavy particles with relaxation times much larger than the Kolmogorov time scale, and drive self-similar clustering at larger scales (Bec et al. 2007;Goto & Vassilicos 2006;Yoshimoto & Goto 2007;Baker et al. 2017).Besides, the inertial-range clustering of heavy particles can be alternatively explained by the sweep-stick mechanism, where inertial particles are assumed to stick to and move with fluid elements with zero acceleration (Goto & Vassilicos 2008).
Apart from preferential concentration, the sling effect or the formation of caustics becomes more significant with the growth of particles inertia or turbulence intensity (Falkovich et al. 2002;Wilkinson & Mehlig 2005;Wilkinson et al. 2006;Falkovich & Pumir 2007).When sling happens, inertial particles encountering intermittent fluctuations get ejected out of local flow.Then clouds of particles with different flow histories interpenetrate each other with large relative velocities (Bewley et al. 2013).As a result, the large relative velocity gives rise to an abrupt growth of the collision frequency, and the sling effect becomes the dominant collision mechanism for large-inertia particles (Voßkuhle et al. 2014).
By means of direct numerical simulations (DNS), the contributions of both preferential concentration and the sling effect (or the turbulence transport effect) to the collision rate could be quantified separately.A framework based on the radial distribution function and the mean relative velocity was proposed to predict the collision kernel function of neutral monodisperse particles with arbitrary inertia (Sundaram & Collins 1997;Wang et al. 2000).Based on this framework, the impacts of turbulence Reynolds number, nonlinear drag, gravity, and hydrodynamic interactions were systematically investigated (Ayala et al. 2008b,a;Ireland et al. 2016a,b;Ababaei et al. 2021;Bragg et al. 2022).In addition to monodisperse particles, this framework was also extended to bidispersed particle systems.Compared to monodisperse system, the clustering of different-size particles is found weaker, while the turbulent transport effect is more pronounced due to particles' different responses to turbulent fluctuations (Zhou et al. 2001).
In addition to neutral particles, particles in nature could easily get charged after getting through ionized air (Marshall & Li 2014;van Minderhout et al. 2021) or encountering collisions (Lee et al. 2015).The resulting electrostatic interaction could significantly modulate the fundamental particle dynamics in turbulence, such as clustering, dispersion and agglomeration (Karnik & Shrimpton 2012;Yao & Capecelatro 2018;Ruan et al. 2021;Boutsikakis et al. 2022;Ruan & Li 2022;Zhang et al. 2023).For identically charged particles, it has been shown that the Coulomb repulsion could effectively repel close particle pairs and mitigate particle clustering.By assuming that the drift velocity caused by Coulomb force can be superposed to the turbulence drift velocity, the radial distribution function of particles with negligible/finite inertia can be derived, which shows good agreement with both experimental and numerical results (Lu et al. 2010a,b;Lu & Shaw 2015).
Different from the like-charged monodisperse system, the real charged particle-laden flows may be more complicated because: (1) the size of suspended particles is generally polydispersed (Zhang & Zhou 2020), and (2) the charge polarity varies with the particle size because of charge segregation.That is, after triboelectrification, smaller particles tend to accumulate negative charges, while large ones are preferentially positively charged (Forward et al. 2009;Lee et al. 2015).In the previous work by Di Renzo & Urzay (2018), DNS of the bidispersed suspensions of oppositely charged particles was performed.It was shown that, the aerodynamic responses of different-size particles to turbulent fluctuations varies significantly.Although the system is overall neutral, different clustering behaviors lead to a spatial separation of positive and negative charge generating mesoscale electric fields in the space.
However, even though the charge segregation is often correlated with the size difference, few studies have been devoted to the dynamics of bidispersed particles with both charge and size difference.In this system, three kinds of electrostatic interactions are introduced simultaneously, i.e., repulsion between small-small/large-large particles and attraction between small-large particles.How would these electrostatic forces affects the clustering of bidispersed particles is less understood.Moreover, the collision frequency of charged particles is of great importance to understand the physics of rain initiation (Harrison et al. 2020) and planet formation (Steinpilz et al. 2020), and the influence of charge segregation on these processes are still unclear.
To address the above issues, in this study we perform direct number simulations to investigate the dynamics of bidispersed oppositely charged particles in homogeneous isotropic turbulence.The numerical methods are introduced in Sec.2.In Sec.3, the statistics of the radial distribution function and the mean relative velocity are first presented to show the effects of electrostatic force on the clustering and the relative motion between different kinds of particle pairs.Then the modulation of collision frequency is quantified using the particle flux, and the mean Coulombturbulence parameter is proposed to characterize the relative importance of electrostatic force compared to particle-turbulence interaction.Finally, the particle flux density in the relative velocity space is discussed, which illustrates in detail how electrostatic force affects the collision of charged particle pairs with different relative velocities.

Simulation conditions
The transport of bidispersed solid particles in turbulence is numerically investigated using the Eulerian-Lagrangian framework.The flow field is homogeneous isotropic turbulence with Re  = 147.5.Solid particles are treated as discrete points and their motions are updated based on Maxey-Riley equation (Maxey & Riley 1983).
Table 1 lists simulation parameters used in this study.Particle properties are chosen based on typical values of silica particles suspended in gaseous turbulence.Here, the particle density is much larger than that of the air ( p / f ∼  (10 3 )), and the particle size is smaller than the Kolmogorov length (i.e.,  p,S / = 0.1,  p,L / ≈ 0.3).The particle Stokes , defined as the ratio of the particle relaxation time  p to the Kolmogorov time scale   , is given by (2.1) In previous field measurement (Zhang & Zhou 2020), the Stokes number of sand particles mainly lies with the range of  = 1 − 10 in strong sandstorms.We thus choose  S = 1 and  L = 10 as the typical Stokes numbers for small and large particles, respectively.In addition, since the measured particle mass concentration in Zhang & Zhou (2020) is low, the particle volume fraction is set small ( par = 6.87 × 10 −6 ) in our simulation, and the force feedback from the dilute particle phase to the fluid phase is omitted (Balachandar & Eaton 2010).
When choosing the particle charge, we assume that particles carry a certain amount of charge as a result of triboelectrification.Both small and large particles are assigned the same amount of charge  but with the opposite polarities.Because of charge segregation, small particles ( S = 1) are negatively charged while large particles ( L = 10) are positive.In previous measurements of tribocharged particles, the surface charge density could reach  p ≈ 10 −5 C/m 2 (Waitukaitis et al. 2014;Lee et al. 2015;Harper et al. 2021).For particle size of  p ∼ 10 m, the order of particle charge can be estimated as  =  2 p  p ∼ 10 −14 C.Moreover, since the considered particle concentration is low, particle collision is negligible.So the charge transfer is not included and the prescribed charge on each individual particle remains constant.

Fluid phase
The direct numerical simulation of the incompressible homogeneous isotropic turbulence is performed using the pseudo-spectral method.The computation domain is a triply periodic cubic box with the domain size  3 = (2 × 10 −2 m) 3 and the grid number  3 grid = 192 3 .The governing equations of the fluid phase are given by and Here, (  ) is the fluid velocity at the grid node   ,  is the pressure,  f is the fluid density, and  f is the fluid kinematic viscosity. F is the forcing term that injects kinetic energy from large scales and keeps the turbulence steady.In the wavenumber space, the forcing term is given by Here  is the constant forcing coefficient, the critical wavenumber is  crit = 5 0 with  0 = 2/ the smallest wavenumber.When evolving the turbulence, the second-order Adam-Bashforth scheme is used for the nonlinear term while the viscous term is integrated exactly (Vincent & Meneguzzi 1991).The fluid time step is Δ F = 1 × 10 −5 s.Table 1 lists typical parameters of the homogeneous and isotropic turbulence.The fluctuation velocity  ′ and the dissipation rate  can be obtained from the energy spectrum  () as 3 2 (2.5) The Kolmogorov length and time scales are given by  = ( f 3 /) 1/4 and   = ( f /) is  e = / ′ , and the Taylor Reynolds number is Re  =  ′ / f with  =  ′ (15 f /) 1/2 the Taylor microscale.

Particle motion
The movements of charged particles are integrated using the second-order Adam-Bashforth time stepping.The particle time step is Δ p = 2.5 × 10 −6 s.This time step is much smaller than the particle relaxation time to well resolve particles' aerodynamic response to the turbulence (Δ p / p,S = 0.23%).Besides, within each step, particles only move a short distance compared to the Kolmogorov length ( ′ Δ p / = 1.4%), so that the variation of the electrostatic force is well captured when particles are getting close.The governing equations of particle motions are d   =   . (2.6b) Here,   and   are the velocity and the location of particle , respectively.  =  p  3 p, /6 is the particle mass, with  p the density and  p, the diameter. F  and  E  are the fluid force and the electrostatic force acting on particle .
Since particles are small ( p < ) and heavy ( p / f ≫ 1), the dominant fluid force comes from the Stokes drag (Maxey & Riley 1983) (2.7) Here  f is the dynamic viscosity of the fluid, and (  ) is the fluid velocity at the particle position.
For two consecutive particle time steps within the same fluid time step, the flow field remains unchanged, and the fluid velocity at the particle location is interpolated at each particle time step using the second-order tri-linear interpolation scheme (Marshall 2009;Qian et al. 2022).The effect of fluid inertia is neglected by assuming that the particle Reynolds number Re p is much smaller than unity.In this dilute system, the average separation between charged particles is large compared to particle size ( l ≫  p ) indicating insignificant effect of particle polarization (Ruan et al. 2022).Therefore, the electrostatic force acting on each particle  can be obtained by summing the pairwise Coulomb force from other particles , which reads (2.8) Here,  0 = 8.8542 × 10 −12 F/m is the vacuum permittivity,   is the point charge located at the centroid of particle .
To properly apply the periodic boundary condition of the long-range electrostatic force, several layers of image domains are added around the original computation domain (Di Renzo & Urzay 2018; Boutsikakis et al. 2023).When calculating the electrostatic force acting on the th particle located at   , the contributions from other particles  at   as well as their images at   + (  +   + ) are all considered: where , ,  = − per , ...,  per ., ,  are unit vectors along the , , and  directions.According to Eq. 2.9, when computing the electrostatic force acting on the target particle , the cost is  ( s ) =  (( p − 1) × (2 per + 1) 3 ), where  ( s ) is the number of source particles to each target particle.The total cost of electrostatic computation then scales with  ( 2 p  3 per ).Considering the large particle number and the addition of image domains, direct summating Eq. 2.9 is extremely expensive.Therefore, we employ the fast multipole method (FMM) to accelerate the calculation.In FMM, the field strength generated by neighboring particles are directly computed using Eq.2.9, while the contribution from clouds of particles far from the target particle is approximated using fast multipole expansion, which reduces the required computation cost to  ( p log( p  3 per )) (Greengard & Rokhlin 1987, 1997).In this study the opensource library FMMLIB3D is employed to conduct fast electrostatic computation (Gimbutas & Greengard 2015).The accuracy of FMM with various layers of image domains is compared in Appendix A. Based on the analysis, two layers of image boxes ( per = 2) are added, which is sufficient to ensure the convergence of the electrostatic force acting on all particles in the original domain.
Moreover, the Coulomb force is capped if the distance between two particles is smaller than a preset critical distance to avoid nonphysically strong force: (2.10) Here the critical distance is  cap =  p,SS /2 or  cap = /20.This value is chosen so that even for a pair of small particles approaching each other, the Coulomb force between them is still accurate at the collision distance (|  −   | =  p,SS >  cap ).Moreover, since  cap is set small compared to the Kolmogorov length scale , capping the electrostatic force does not affect the statistics at larger scales presented below.
In each run, the single-phase turbulence is first evolved until reaching the steady state.Particles are then injected randomly in the domain with their initial velocities equal to the fluid velocity at particle locations.It takes roughly 3.5 e for particle dynamics to reach equilibrium, statistics are then taken over another 5 e .Three parallel runs with different initial particle locations are performed for each case, and their results are averaged and presented below.

Clustering and relative velocity of charged particles
Figure 1 compares the spatial distribution of bidispersed particles within a thin slice || ⩽ 20.Although oppositely-charged bidispersed particles are simulated in each case, we show small and large particles in the left and the right panels separately for better illustration.In the neutral case (Fig. 1(a) and (b)), particle behavior is solely determined by the particle-turbulence interaction.Small particles ( = 1) are responsive to fluctuations at the Kolmogorov length scale.As a result, their spatial distribution is highly nonuniform, and small-scale particle clusters can be observed We employ the radial distribution function (RDF) to quantify how particles from the same (or different) size groups cluster.The radial distribution functions of different kinds of particle pairs, i.e., small-small (SS), large-large (LL), and small-large (SL) paris, are defined as In the definition of  SS (), Δ SS () is the number of SS pairs with their separation distances lie within the range of  ± Δ/2, and Δ () = 4[( + Δ/2) 3 − ( − Δ/2) 3 ]/3 is the shell volume in the separation distance bin. p,S ( p,S − 1)/2/ 3 is the average number of SS pairs in the whole domain.Therefore,  SS () > 1 means there is a higher density of SS pairs at a given separation distance of  compared to the average pair density.The definitions of  LL () and  SL () are similar and thus omitted.
As shown in Fig. 2(a),  SS () for neutral pairs is significantly larger than unity at  ⩽ , indicating strong clustering at the small scale. SS () also follows a clear power law  SS () ∼ (/) −c 1 with the fitted exponent c 1 = 0.69.The value of c 1 is close to those reported in previous studies with similar Re  , e.g., c 1 = 0.67 for Re  = 143 (Saw et al. 2012) and c 1 = 0.68 for Re  = 140 (Ireland et al. 2016a).In comparison,  LL () for neutral pairs is only slightly larger than one, which means the clustering of large particles are relatively weak.However,  LL () decreases slowly with the increase of  and stays above unity until  ≈ 100, while  SS () drops rapidly and approaches unity around 20 − 30.This is because the clustering of inertial particles with the relaxation time  p are driven by vortices whose timescales are comparable to particle's relaxation time (i.e., () ∼  p ).In the inertial range, the time scale satisfy () ∼  −1/3  2/3 , so the relationship can be given as  ∼  1/3  3/2 p (Yoshimoto & Goto 2007;Bec et al. 2011).With the increase of particle inertia, the length and the time scales of vortices that could affect particle behaviors also become larger.As a result, the size of particle clusters also grow larger, leading to the large correlation length in  LL () (Yoshimoto & Goto 2007;Liu et al. 2020).Besides, when comparing those particles with large St difference, because their relaxation times differ by a factor of ten, they respond to flows of very different scales.Therefore, there is little spatial correlation between small and large particles, and  SL () is close to unity at all separation distance  when particles are neutral.
When particles are charged, since charge segregation correlates with size, the same-size particles will repel each other when they get close.Since the amount of particle charge is the same, the effect of repulsion is more drastic between SS pairs rather than between LL pairs.As a result, the order of  SS () at small  drops by several decades as  increases.The same decreasing trend is observed for  LL (), but the influence is less significant because of the large particle inertia.This can also be shown by comparing the top and the bottom panels in Fig. 1.In the charged case ( = 2 × 10 −14 C), the small particles become less concentrated (Fig. 1(c)), while no obvious difference is seen in the spatial distribution of large particles (Fig. 1(d)).It is worth noting that, since the Coulomb force decays rapidly with , its effect is only significant within a relatively short range ( ≈ 1 − 10).Beyond this range, both  SS () and  LL () in Fig. 2 (a) and (b) recover to their neutral values again, so clustering could still be observed at large .As for  SL , because of the Coulomb attraction, particles of different sizes are more likely to stay close.More specifically, considering the large mass difference between SL pairs ( L / S = ( L / S ) 1.5 ≈ 30), small particles will always get attracted towards the large ones, giving rise to the rapid growth of  SL () at small separation ( ⩽ ).Again, the opposite-sign attraction decays with increasing , so  SL () approaches one when  is sufficiently large.
Apart from spatial correlation, it is also of interest to study the relative velocity between particle pairs and focus on the modulation caused by the electrostatic force.particle charge .For neutral SS pairs, when  is in the inertial range ( ≪  ≪  int ),  p,S is much smaller than the characteristic time scales of turbulent fluctuations at this length scale.The radial relative velocity between SS pairs, ⟨| r,SS |⟩, thus follows that of fluid tracers  r .Here, the fluid relative velocity is obtained from the the second order longitudinal structure function as  r =  2 1/2  1/3  1/3 and the constant is  2 = 2.13 (Yeung & Zhou 1997) (shown as dash-dotted lines in Fig. 3).If  is within the dissipative range, the timescale of local fluctuations becomes comparable to  p,S .Small neutral particles cannot perfectly follow the background flow at this separation, so ⟨| r,SS |⟩ deviates from the fluid relative velocity  r = (/15) 1/2  (dashed lines in Fig. 3).
When particles are charged, since the influence of Coulomb force is negligible at large , curves for different  collapse.In contrast, ⟨| r,SS |⟩ is found to rise significantly as particles get close to each other.When  gets larger, the increase of ⟨| r,SS |⟩ becomes more pronounced and the deviation from the neutral result also occurs at a larger .This seems counter-intuitive because the Coulomb repulsion between SS pairs should always slow approaching particles down and reduce the relative velocity.One explanation for this change is as follows.When  is large, the strong Coulomb repulsion causes biased sampling at small : only those pairs with large inward relative velocity could overcome the energy barrier and get close.Meanwhile, since the electrical potential energy is conserved in the approach-then-depart process, the pairs approaching each other at a high velocity will also be repelled away at a high velocity as the electrostatic force pushing them apart.As a result, pairs with a large | r,SS | take a large proportion as  increase, so the mean radial relative velocity between both approaching and departing pairs, ⟨ r,  ⟩, shows the increasing trend in Fig. 3(a).This increasing trend with  will be further discussed in Sec 3.3.
Compared to small particles with  = 1, large particles with  = 10 are very inertial and insensitive to fluctuations even at large separations.The curves of ⟨| r,LL |⟩ are therefore much flatter.Besides, the constant relative velocity between LL pairs at small  is also significantly larger, because their relative velocity comes from the energetic large-scale motions.When particles are charged, the same increasing trend with  is also observed but less significant, which can be attributed to the large kinetic energy compared to the electrostatic potential energy.
As for SL pairs, ⟨| r,SL |⟩ is larger than both ⟨| r,SS |⟩ and ⟨| r,LL |⟩ in the neutral cases.This is caused by the different responses of small/large particles to turbulent fluctuations, which is termed as the differential inertia effect (Zhou et al. 2001).Interestingly, there is no obvious change when comparing curves between different charge  in Fig. 3(c).As shown in Fig. 2(c), Coulomb force attracts small particles towards large ones and enhances their spatial correlation.Nevertheless, even though charged SL pairs experience a more similar flow history than neutral SL pairs, they will still develop a large relative velocities over time because their response to local fluctuations is very different.Therefore, the influence of charge on ⟨| r,SL |⟩ is weak.

Effect of Coulomb force on particle flux
Once the radial distribution function and the radial relative velocity are known, we could further investigate the collision frequency between charged particles.For a steady system, the collision frequency can be measured by the kinematic collision kernel function as (Sundaram & Collins 1997;Wang et al. 2000) Even though other particle parameters (e.g., the Stokes number) are set the same, different choices of  C will affect the final outcome of Γ   by changing the collision geometry.Therefore, instead of directly using Γ   , we define the particle flux Φ   as the ratio of the collision kernel Γ   () to the area of the collision sphere 4 2 at  as: Φ   can be understood as the number of particles crossing the collision sphere per area per unit time, which is independent of the prescribed  C .Apart from Eq.3.3 that uses information of both approaching and departing pairs, the particle flux can also be defined using only the approaching or departing pairs.In real situations, it is the approaching pairs that lead to collisions, so a natural way to define particle flux is based on the inward flux as Here,  ( r,  < 0|) is the fraction of particles moving inwards at , and the mean inward relative velocity is ⟨ − r,  ⟩ = − ∫ 0 −∞  r,     ( r,  )d r,  .After the system reaches the equilibrium, RDFs between particle pairs no longer change, suggesting that the inward flux should balance the outward flux at all .Here, the outward particle flux is with the mean outward relative velocity given by ⟨ + r,  ⟩ = ∫ ∞ 0  r,     ( r,  )d r,  .The definition of Φ   is based on a pair of particles  and  without specifying their sizes.If the average is taken over all SS pairs, the result is the small-small particle flux denoted by Φ SS .The flux between LL pairs (Φ LL ) and SL pairs (Φ SL ) can also be obtained by taking the average over corresponding particle pairs.
The SS fluxes Φ SS defined by Eqs.3.3, 3.4 and 3.5 are first compared and show good agreement with each other, indicating that the flux-balance condition is valid.Since Eq.3.3 uses the information of both approaching and departing particle pairs, it is adopted in the following sections for better statistics.
Fig. 4(a) compares the SS fluxes with various .The flux for neutral pairs remain constant when  is comparable to .This indicates that, although both  SS (Fig. 2(a)) and ⟨| r,SS |⟩ (Fig. 3(a)) keep changing as  decreases, their product is almost constant for small .For inertial particle pairs with small , it has been shown that () ∝   2 −3 with  2 the correlation dimension (Bec et al. 2007), while the relative velocity follows ⟨ r,SS ⟩ ∝  3− 2 (Gustavsson & Mehlig 2014).Therefore, the dependence of Φ SS on  is canceled out.Besides, in practical situations the collision diameter  C between micron particles/droplets is smaller than , so the constant flux at such separation distance leads to a quadratic dependence of collision kernel on  C as When particles are charged, the SS flux is found to decrease rapidly as  drops.To characterize the effect of Coulomb force on the reduction of Φ SS , we need to quantify the competition between the driving force and the resistance of particle collisions.For a pair of same-sign particles  and  separating by , the driving force can be evaluated by the relative kinetic energy  Kin () =    ⟨| r,  |⟩ 2 /2.Here,    = ( −1  +  −1  ) −1 is the effective mass, ⟨| r,  |⟩() is the mean radial relative velocity between neutral pairs with a separation of .The resistance is the electrical energy barrier  Coul =     /4 0 .The energy ratio is then given as (3.6)Note that in previous studies regarding the clustering of charged particles with weak inertia (Lu et al. 2010a,b), the approaching velocity between particle pairs can be directly derived through perturbation expansion in Stokes number (Chun et al. 2005).In this work, however, particles are very inertial, so the mean relative velocity ⟨| r,  |⟩ between neutral pairs is used instead.By using the corresponding effective mass (e.g.,  SS =  S /2,  LL =  L /2, and  SL =  S  L /( S +  L )) and the mean relative velocity (⟨| r,SS |⟩, ⟨| r,LL |⟩, and ⟨| r,SL |⟩), Eq. 3.6 can quantify the Coulomb-turbulence competition for different kinds of particle pairs.Since Ct 0 measures the relative importance of the Coulomb force to the mean relative kinetic energy, it is called the mean Coulomb-turbulence parameter hereinafter.
According to the definition in Eq.3.6, the value of Ct 0 can be varied as the particle charge or the separation distance changes.Therefore, for each data point in Fig. 4 (a), we plot the ratio of Φ SS () to its corresponding value in the neutral case Φ SS ( = 0) as a function of Ct 0 in Fig. 5(a).The flux ratios for different  are found to follow the same trend.The particle flux of charged particles is almost unaffected when Ct 0 is small, and a significant decay is observed when Ct 0 exceeds unity.The LL fluxes Φ LL in Fig. 4(b) are analyzed in a similar way, and the results are plotted in Fig. 5(b).Because of the large kinematic energy between LL pairs, the Coulomb repulsion is relatively weak (Ct 0 < 1), so the reduction of Φ LL has not yet reached the electrostatic-dominant regime.As for the flux between SL pairs shown in Fig. 4(c), the opposite trend is seen when plotting the flux ratio as a function of Ct 0 (Fig. 5(c)).The increase of Φ SL becomes significant when Ct 0 becomes larger than one.
We now propose a model to quantify the influence of the electrostatic force on the mean particle Blues solid lines are model predictions from Eqs. 3.16 and 3.17 with fitted . Brown dash-dotted lines are predictions with the fixed  = 1.
flux.For particle pairs with a separation distance , the mean radial relative velocities between different kinds of particle pairs have been shown in Fig. 3.We then assume that, for each kind of particle pairs, the relative velocity between neutral pairs follows the Gaussian distribution.
Taking SS particle pairs as an example, the probability density function  r,SS of the relative velocity  r,SS () is The standard deviation is determined by the mean relative velocity at  as  SS = √︁ /2 • ⟨| r,SS |⟩().It should be noted that, the relative velocity distribution is actually non-Gaussian (Sundaram & Collins 1997;Wang et al. 2000;Ireland et al. 2016a).However, a Gaussian distribution is sufficient for the first-order assumption, which has already been applied in previous studies (Pan & Padoan 2010;Lu & Shaw 2015).
We then evaluate the RDF of charged pairs.For neutral small-small pairs, all particle pairs with an inward relative velocity ( r,SS < 0) could approach each other.In comparison, when they are identically charged, only those SS pairs whose inward relative velocity exceeds a critical velocity ( r,SS < − C,SS ) are able to get close.By balancing the Coulomb energy barrier and the relative kinetic energy, the critical velocity can be given as (3.8) The critical velocity for LL or SL pairs can be given by replacing  SS =  S /2 with  LL =  L /2 or  SL =  S  L /( S +  L ).Note that the critical velocity is derived from the interaction between a pair of particles, which could reasonably describe the major effect of Coulomb force in the current dilution suspension.If the particle concentration becomes higher, the particle number density needs to be included for a more accurate prediction (Boutsikakis et al. 2022(Boutsikakis et al. , 2023)).By assuming a sharp cut-off at  = − C,SS , Lu & Shaw (2015) related the RDF of charged SS pairs  SS ( |) to that of neutral SS pairs  SS ( | = 0) as where erf is the error function.However, as the relative velocity magnitude  r becomes smaller than  C,SS , the corresponding flux contribution does not drop sharply to zero.Instead, a smooth transition would be expected, so the ratio of the charged RDF to the neutral RDF should be written as Here, Θ SS () is the electrical kernel that gradually transients from one to zero when the magnitude of  r drops below  C,SS .The proportion of charged particles that could overcome the Coulomb repulsion is then obtained from the convolution of Θ SS and  r,SS , which is the numerator of the right-hand term in Eq.3.10.To account for the smooth transition without adding significant complexity, we still adopt the sharp cut-off assumption, but the effective cut-off velocity  SS  C,SS is used.Here,  SS ∼  (1) is the correction factor that ensures an accurate prediction as: (3.11)The difference between the effective cut-off and the actual electrical kernel will be further discussed in Sec.3.3.Consequently, Eq.3.9 becomes To obtain the mean relative velocity between charged SS pairs, we start by computing the mean relative velocity in the velocity interval [−∞, − SS  C,SS ] between neutral SS pairs: Then the mean relative velocity of charged SS pairs can be obtained by subtracting the Coulomb potential energy from the mean relative kinetic energy in the velocity interval Here, the correction factor  SS is also added to the last term on the right-hand side to account for the smooth transition.Taking Eq.3.13 into Eq.3.14 then yields (3.15) Finally, by multiplying Eqs.3.9 and 3.15 and taking into the definition of Ct 0 (Eq.3.6), the flux ratio can be given as The flux ratios for LL pairs and SL pairs can be derived in a similar way and are written as Eqs.3.16 and 3.17 are then fitted as blue lines in Fig. 5, which show good agreement with the simulation results.The fitted correction factors are  SS = 0.193,  LL = 0.739, and  SL = 0.5415, satisfying the assumption that the correction factors are of the order of unity.The predictions with the fixed  = 1 are also shown as brown dash-dotted lines, which correspond to the original model that assumes the sharp cut-off occurs at the critical velocity  C,SS/LL/SL .Although the trends are similar, models with the fixed  = 1 underestimate the critical Ct 0 at which the transition occurs.Moreover, the proposed model underestimates the SS flux when Ct 0 ≈ 10 2 in Fig. 5(a).To find the origin of this discrepancy, the predicted    () using Eq.3.12 and the predicted ⟨| r,SS |⟩() using Eq.3.15 are also plotted as grey dashed lines in Fig. 2 and Fig. 3, respectively.It can been seen that, the predicted RDFs are comparable to the simulation results, so the discrepancy of Φ SS at Ct 0 ≈ 10 2 mainly comes from the underestimation of the mean relative velocity at  ∼  in Fig. 3(a).Since the PDF of  r is assumed Gaussian in our model, the influence of intermittency is not considered.As a result, the proportional of particle pairs with large relative velocity is significantly underestimated.

Particle flux density for charged particles
In Sec.3.2, the mean Coulomb-turbulence parameter Ct 0 is defined using the mean relative velocity ⟨| r,  |⟩, which compares the importance of Coulomb force to the mean relative motion caused by turbulence.However, it has been known that the relative velocity between particle pairs may vary within a wide range, and the mean Coulomb-turbulence parameter Ct 0 may not be sufficient to fully reveal the physics.
In this section it would be of our interest to examine the impacts of Coulomb force on particle pairs with different relative velocities.For different kinds of particle pairs, the particle flux Φ   in the relative velocity space is expanded as (3.18) Here,    is the density of particle flux within each relative velocity interval, which measures the contribution to the total particle flux Φ   by particle pairs whose relative velocity is within the interval  r ± Δ r /2.By comparing Eqs.3.18 and 3.3,    can be given as We start with the effect of Coulomb force on the PDFs of relative velocity    ( r,  ), because  the distribution of    in Eq.3.19 is strongly dependent on    ( r,  ).Fig. 6(a) illustrates PDFs of  r,SS at the separation distance interval  ∈ [1.15, 2.5).The Gaussian distribution defined by Eq. 3.7 is also shown as the black dashed line.By comparing the neutral PDFs and the Gaussian curves, it is clear that the Gaussian assumption serves as a reasonable approximation at | r | < 3  , but significantly underestimates the probability of large | r |.Therefore, the Gaussian assumption is only suitable for modeling low-order effects, not for higher-order statistics.Besides, as will be discussed below, since the Gaussian curve is symmetric, it is not able to capture the asymmetry of the relative velocity between approaching and departing pairs.For SS pairs, compared to the neutral PDF, charged PDFs become wider as  increases, indicating a lower/higher probability of finding particle pairs with low/high relative velocity within the separation interval.This is consistent with the increase of the mean relative velocity with  in Fig. 3(a).If we look at the symmetry of  SS , the neutral PDF is found negatively-skewed.This asymmetry is attributed to: (1) the fluid velocity derivative is negatively-skewed, and the relative motion between particle pairs with  = 1 is partially coupled with local flow and thus exhibits a similar feature (Van Atta & Antonia 1980;Wang et al. 2000); (2) the asymmetry of particle path history gives rise to larger relative velocity between approaching pairs than departing pairs (Bragg & Collins 2014).However, the asymmetry of  SS curves are significantly reduced once particles are charged (see Table 2), which results from the symmetric nature of the Coulomb force.The magnitude of Coulomb force relies only on the interparticle distance , so the approaching or departing SS pairs experience the same amount of repulsion as long as  is the same, which makes the approaching-then-departing process more symmetric.Therefore, introducing Coulomb force could effectively increase the standard deviation but reduce the skewness of  r,SS .
For LL pairs, since Coulomb repulsion is only able to repel pairs with small relative velocity,  LL drops slightly as  r approaches zero, but remains almost unchanged at larger  r .As for  SL , the difference between neutral and charged case is insignificant, which again demonstrates that the velocity difference between particles of different sizes mainly comes from the differential inertia effect, while the influence of Coulomb force is limited.Besides, the movement of large particles with  = 10 has very weak couplings with local flow fields, so the skewness of  LL and  SL in Table 2 are negligible.
We now turn to the distribution of flux density    in the velocity space.Fig. 7(a),(c),(e) plot the flux densities between SS, LL, and SL pairs, respectively.Different from the unimodal PDFs of the relative velocity shown in Fig. 6,    are bimodal because the maximum flux density is reached when the product | r,  | •    ( r,  ) is the largest.Therefore, it is the particle pairs with the intermediate relative velocity that contributes the most to the overall particle flux.
To better describe the impacts of Coulomb force, we define the Coulomb kernel Θ   as the ratios of the charged flux densities to their neutral values, i.e., Θ   () =    ()/   ( = 0), which are displayed in the right panel of Fig. 7.As shown in Fig. 7(b), the value of Θ SS varies by more than one order of magnitude, indicating that the influence of Coulomb force changes significantly as  r changes.When | r | is small, Θ SS decreases drastically as  increases because Coulomb repulsion dominates.In addition, Θ SS is found to be asymmetric.As discussed above, neutral particles with  = 1 tend to separate at a lower relative velocity compared with the approaching pairs.When they are charged, particle pairs originally separating at low speeds will be accelerated and pushed apart by Coulomb repulsion, which effectively shifts the high flux density from the small positive  r to a larger positive  r in the velocity space.Consequently, Θ SS experiences a more significant decrease at  r slightly greater than zero, but quickly exceeds one when  r ≈ 2  .Moreover, Θ SS becomes larger than one at large | r | for both approaching and departing pairs.The explanation is that, with the increase of , small particles are more likely to get attracted towards the large ones and follow their motions.Since the relative velocity between LL pairs is generally much larger, this effect could increase the SS flux density at large | r |.However, the contribution of the increased flux at large positive | r | is negligible compared to the decrease at small | r | (see Fig. 7(a)), so the major effect of Coulomb force is to reduce the SS flux through the small-small repulsive force.
The kernel Θ LL between LL pairs (Fig. 7(d)) also drops quickly at small | r |, but recovers to unity when | r | ⩾ 2  .Besides, because of the limited effects of local fluid velocity gradient mentioned above, the approaching and departing processes are more symmetric for neutral LL pairs, and adding Coulomb force still retains this symmetry.As for Θ SL shown in Fig. 7(f), the change is still most significant at small  r = 0.However, different from the kernels of SS/LL pairs where the impact of the same-sign repulsion is only obvious within a narrow interval (i.e., | r |/  ⩽ 2), the opposite-sign attraction seems to increase Θ SL in a much wider range of  r .For instance, Θ SL is increased even at  r /  = −10 for  = 2 × 10 −14 C in Fig. 7(f).As discussed in Sec.3.2, the main effect of the opposite-sign attraction is to enhance small-large clustering.Then particles of different sizes, though staying close to each other, will develop a large relative velocity as a result of their different responses to local fluctuations, which leads to the increase of Θ SL in a wide range of  r .
The discussion above has shown that, the influence of Coulomb force is different as the particle relative velocity  r,  changes.Therefore, instead of using the mean relative velocity ⟨| r,  |⟩ (Eq.3.6), we adopt the median relative velocity  r,  in each  r bin to estimate the relative kinetic energy.For a given separation distance  and a certain relative velocity bin with the median  r,  , the extended Coulomb-turbulence parameter is given as (3.20) We then examine the dependence of the electrical kernel Θ   on the extended Coulombturbulence parameter Ct.For the data points ( r,  , Θ   ) shown in the right panel of Fig. 7, the corresponding values of Ct are calculated using Eq.3.20, and the results are re-plotted as points (Ct, Θ   ).Note that to ensure meaningful statistics, only data points satisfying certain criterion are considered and the reasons are as follows: (1) In addition to the separation interval  ∈ [1.15, 2.5) shown in Fig. 7, data from two more separation intervals, i.e., [0.85, 1.15) and [2.5, 3.5), are also added.At these separation distances, the effect of electrostatic force is already significant, while the number of samples is large enough for statistics.(2) We only consider pairs with their relative velocity in a certain range to avoid using the more scattered data points at large  r .For SS pairs, the relative velocity range is [−5  , 5  ], while for LL/SL pairs the range considered is [−10  , 10  ].The particle flux within the above ranges contribute to at least 97.1%/95.6%/96.3% of the total SS/LL/SL particle flux in the neutral case, which reflects the major change of Φ   in each case.(3) As shown in Fig. 7(b), Θ SS is not symmetric about  r = 0.When particles are departing, Coulomb repulsion will redistribute  SS in the velocity space, which may distort the Θ SS −Ct relationship.We therefore only use the data from approaching pairs ( r < 0) for later analysis.
Since the outward flux is balanced by the inward flux in the steady state, the total flux could still be evaluated as Fig. 8 plots the dependence of Θ   on Ct for different particle pairs.Despite different particle charge  and separation , the data points for both SS (Fig. 8(a)) and LL (Fig. 8(b)) pairs show a similar trend.When Ct is small, Coulomb force is weak compared with particle-turbulence interaction, so Θ SS and Θ LL stay close to one.When Ct becomes larger than one, the same-sign repulsion starts to significantly reduce the corresponding particle flux density, and a decrease of Θ is observed.In Sec.3.2, it has been assumed that a sharp cut-off occurs at the effective velocity  SS  C,SS or  LL  C,LL for SS or LL pairs.However, the transition of Θ SS and Θ LL in Fig. 8(a) turn out to be smooth.Thus, the sharp cut-off can be understood as the first-order approximation of the influence of Coulomb repulsion, which can be written as By taking into the fitted results ( SS = 0.193 and  LL = 0.739) from Section 3.2, the critical value can be given as Ct crit,SS = 26.85 and Ct crit,LL = 1.83, respectively.However, different from Θ SS and Θ LL , there is no clear dependence of Θ SL on Ct in Fig. 8(c).The reason is as follows.The opposite-sign Coulomb force will attract SL pairs with a low departing velocity together and promote clustering.However, although these SL pairs have a low relative velocity at first, they will develop a much larger relative velocity over time.As a result, the relative kinetic energy  Kin becomes independent of the electrical potential energy  Coul , and no clear relationship can be found between Θ SL and  SL .

Conclusions
In this study, we investigate the effects of charge segregation on the dynamics of tribocharged bidispersed particles in homogeneous isotropic turbulence by means of DNS.Using radial distribution function    (), we show that Coulomb repulsion/attraction significantly inhibits/promotes particle clustering within a short range, while the clustering at a large separation distance is still determined by particle-turbulence interaction.For same-sign particles, Coulomb repulsion repels approaching pair with low relative velocity, giving rise to the increase of mean relative velocity ⟨| r,  |⟩ as the separation distance  decreases.In comparison, the relative velocity between different-size particles is almost unchanged for all , which suggests that the differential inertia effect contributes predominantly to ⟨ r,SL ⟩.
By defining the particle flux Φ   as the number of particles crossing the collision sphere per area per unit time, we are able to quantify the particle collision frequency without a prescribed collision diameter  C .As expected, the Coulomb repulsion/attraction is found to reduce/increase the total particle flux Φ   when particles are close.When plotted as a function of the mean Coulombturbulence parameter Ct 0 that measures the relative importance of electrostatic potential energy to the mean relative kinetic energy, the particle flux ratio Φ   ()/Φ   ( = 0) is shown to follow a similar trend.Specifically, by assuming that the PDF of relative velocity follows a Gaussian distribution, the particle flux can be well modelled by a function of Ct 0 .The total particle flux Φ   is then expanded in the relative velocity space as the flux density    , and the relative change Θ   =    ()/   ( = 0) (also termed the Coulomb kernel) in each relative interval exhibits a significant difference.Finally, the extended Coulomb-turbulence parameter Ct is defined using the median  r,  in each relative velocity bin, which better describes the competition between Coulomb force and the turbulence effect.Then for same-size particle pairs, a clear relationship is found between Θ   and Ct.And the smooth transition of Θ   is observed near the critical value Ct crit,  = 1/ 2   .For SL pairs, however, the relative velocity will grow larger because of the predominant differential inertia effect, so Θ SL shows no clear dependence on  r,SL (and therefore on Ct). per 0 1 2 3 4  FMM 1.12 × 10 −1 2.13 × 10 −3 1.70 × 10 −3 1.65 × 10 −3 1.63 × 10 −3 Here,  is the Ewald parameter, erfc is the complimentary error function, and  ′ = 1 is the relative dielectric constant of surrounding medium.
Since Ewald summation is computationally expensive, a smaller-scale system with  p = 2000 oppositely charged particles is used in validation cases.In each case, ten parallel computations with different particle locations are performed to ensure reliable statistics.Table 3 lists parameters used in Ewald summation.The dimensionless product  C equals  to ensure high accuracy in both real and Fourier spaces.The cut-off radius( C )/wavenumber( C ) in the real/Fourier space is then determined by  C = ( C )/ 1/2  1/6 p and  C = 1.8(C ) 2 / C to balance the computation cost of  (r)   and  (k)  (Fincham 1994).When performing FMM computation, the layer number  per is varied from 0 to 4 to show the influence of adding image domains.The relative error of FMM compared to Ewald summation is given as where the norm of force is defined as the root mean square of the force components following Deserno & Holm (1998).The dependence of  FMM on  per is shown in Table 4.The relative error is significant ( FMM > 10%) if periodicity is no considered.After adding image domains, the relative error drops significantly and almost saturates after  per ⩾ 2. Hence,  per = 2 is chosen to guarantee sufficient accuracy ( FMM = 0.17%) while avoiding expensive computation cost.

Figure 2 .
Figure2.RDFs between (a) small-small, (b) large-large, and (c) small-large particle pairs with different particle charge .The inset in (a) shows the full scale of  SS () in the neutral case.Dashed lines from light to dark correspond to predicted RDFs using Eq.3.12 (or its counterpart for the oppositely-charged case) for  = 5 × 10 −15 C, 1 × 10 −14 C, and 2 × 10 −14 C, respectively.

Figure 3 .
Figure3.Mean radial relative velocity ⟨| r |⟩ normalized by   between (a) small-small, (b) large-large, and (c) small-large particle pairs with different particle charge .The inset in (a) shows the full scale of⟨| r,SS |⟩ in the neutral case.The red dashed line denotes the velocity scaling ∝  in the dissipation range, and the red dash-dotted line denotes the velocity scaling ∝  1/3 in the inertial range.Dashed lines from light to dark correspond to predictions using Eq.3.15 (or its counterpart for the oppositely-charged case) for  = 5 × 10 −15 C, 1 × 10 −14 C, and 2 × 10 −14 C, respectively.
.2) Here,  C =   +   is the collision radius, which equals to the sum of the radii of particles  and .   ( C ) is the RDF at  =  C , and the mean relative velocity is ⟨| r,  |⟩ = ∫ ∞ −∞ | r,  |    ( r,  )d r,  with    ( r,  ) the probability density function (PDF) of  r,  at  =  C .

Figure 4 .
Figure 4. Particle fluxes between (a) small-small, (b) large-large, and (c) small -large particles pairs with different particle charge .

Figure 5 .
Figure5.Ratios of charged particle flux Φ() to neutral particle flux Φ( = 0) as a function of the mean Coulomb-turbulence parameter Ct 0 for (a) small-small, (b) large-large, and (c) small-large particle pairs.Blues solid lines are model predictions from Eqs. 3.16 and 3.17 with fitted . Brown dash-dotted lines are predictions with the fixed  = 1.

Figure 6 .
Figure 6.PDFs of the relative velocity between (a) small-small, (b) large-large, and (c) small-large pairs at  ∈ [1.15, 2.5).Zoom-in views around  r = 0 are shown in the insets.Black dashed lines denote the Gaussian distributions defined by Eq. 3.7 for different kinds of particle pairs.
) with  (•) being the Heaviside step function.The critical value of Ct describing where Θ   steps down can be related to the fitted correction factors for the effective cut-off velocity as Ct crit,  = 1/ 2   .(3.22)

Table 4 .
Relative errors of FMM computation compared to Ewald summation.