3.1. Clustering and relative velocity of charged particles

Figure 1 compares the spatial distribution of bidispersed particles within a thin slice $|z| \leq 20 \eta$. 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 (figure 1a,b), particle behaviour is solely determined by the particle–turbulence interaction. Small particles ($St=1$) are responsive to fluctuations at the Kolmogorov length scale. As a result, their spatial distribution is highly non-uniform, and small-scale particle clusters can be observed (figure 1a). Meanwhile, large particles ($St=10$) are more inertial, so they are more dispersed in the domain (figure 1b).

Figure 1.

Snapshots of particles in the slice $|z| \le 20 \eta$ with (a) $St=1$ and $q = 0\ \mathrm {C}$, (b) $St=10$ and $q = 0\ \mathrm {C}$, (c) $St=1$ and $q = 2 \times 10^{-14}\ \mathrm {C}$, (d) $St=10$ and $q = 2 \times 10^{-14}\ \mathrm {C}$. Blue/red dots represent small/large particles. The colour bar indicates the vorticity magnitude normalized by its average over the domain.

We employ the radial distribution function (RDF) to quantify how particles from the same (or different) size groups cluster. The RDFs of different kinds of particle pairs, i.e. small–small (SS), large–large (LL) and small–large (SL) pairs, are defined as

(3.1a)$$\begin{gather} g_{{SS}}(r) = \frac{\Delta N_{{SS}}(r) / \Delta V(r)}{N_{{p,S}}(N_{{p,S}}-1)/2/L^3},\quad g_{{LL}}(r) = \frac{\Delta N_{{LL}}(r) / \Delta V(r)}{N_{{p,L}}(N_{{p,L}}-1)/2/L^3}, \end{gather}$$

(3.1b)$$\begin{gather}g_{{SL}}(r) = \frac{\Delta N_{{SL}}(r) / \Delta V(r)}{N_{{p,S}}N_{{p,L}}/L^3}. \end{gather}$$

In the definition of $g_{{SS}}(r)$, $\Delta N_{{SS}}(r)$ is the number of SS pairs with separation distances within the range of $r \pm \Delta r/2$, and $\Delta V(r) = 4 {\rm \pi}[(r+\Delta r/2)^3-(r-\Delta r/2)^3]/3$ is the shell volume in the separation distance bin; $N_{{p,S}}(N_{{p,S}}-1)/2/L^3$ is the average number of SS pairs in the whole domain. Therefore, $g_{{SS}}(r)>1$ means there is a higher density of SS pairs at a given separation distance of $r$ compared with the average pair density. The definitions of $g_{{LL}}(r)$ and $g_{{SL}}(r)$ are similar and thus omitted.

As shown in figure 2(a), $g_{{SS}}(r)$ for neutral pairs is significantly larger than unity at $r \leq \eta$, indicating strong clustering at the small scale; $g_{{SS}}(r)$ also follows a clear power law $g_{{SS}}(r) \sim (r/ \eta )^{-{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}_{\lambda }$, e.g. ${c}_1 = 0.67$ for ${Re}_{\lambda }=143$ (Saw et al. Reference Saw, Salazar, Collins and Shaw2012) and ${c}_1 = 0.68$ for ${Re}_{\lambda }=140$ (Ireland et al. Reference Ireland, Bragg and Collins2016a). In comparison, $g_{{LL}}(r)$ for neutral pairs is only slightly larger than one, which means the clustering of large particles is relatively weak. However, $g_{{LL}}(r)$ decreases slowly with the increase of $r$ and stays above unity until $r \approx 100 \eta$, while $g_{{SS}}(r)$ drops rapidly and approaches unity around $20\unicode{x2013}30 \eta$. This is because the clustering of inertial particles with the relaxation time $\tau _{{p}}$ is driven by vortices whose time scales are comparable to a particle's relaxation time (i.e. $\tau (r) \sim \tau _{{p}}$). In the inertial range, the time scale satisfies $\tau (r)\sim \epsilon ^{-1/3} r^{2/3}$, so the relationship can be given as $r \sim \epsilon ^{1/3} \tau _{{p}}^{3/2}$ (Yoshimoto & Goto Reference Yoshimoto and Goto2007; Bec et al. Reference Bec, Biferale, Cencini, Lanotte and Toschi2011). With the increase of particle inertia, the length and the time scales of vortices that could affect particle behaviours also become larger. As a result, the size of particle clusters also grow larger, leading to the large correlation length in $g_{{LL}}(r)$ (Yoshimoto & Goto Reference Yoshimoto and Goto2007; Liu et al. Reference Liu, Shen, Zamansky and Coletti2020). Besides, when comparing those particles with a 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 $g_{{SL}}(r)$ is close to unity at all separation distance $r$ when particles are neutral.

Figure 2.

The RDFs between (a) SS, (b) LL and (c) SL particle pairs with different particle charge $q$. The inset in (a) shows the full scale of $g_{SS}(r)$ in the neutral case. Dashed lines from light to dark correspond to predicted RDFs using (3.12) (or its counterpart for the oppositely charged case) for $q=5 \times 10^{-15}\ \mathrm {C}$, $1 \times 10^{-14}\ \mathrm {C}$ and $2 \times 10^{-14}\ \mathrm {C}$, respectively.

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 $g_{{SS}}(r)$ at small $r$ drops by several decades as $q$ increases. The same decreasing trend is observed for $g_{{LL}}(r)$, 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 figure 1. In the charged case ($q=2 \times 10^{-14}\ \mathrm {C}$), the small particles become less concentrated (figure 1c), while no obvious difference is seen in the spatial distribution of large particles (figure 1d). It is worth noting that, since the Coulomb force decays rapidly with $r$, its effect is only significant within a relatively short range ($r \approx 1\unicode{x2013}10 \eta$). Beyond this range, both $g_{{SS}}(r)$ and $g_{{LL}}(r)$ in figure 2(a,b) recover to their neutral values again, so clustering could still be observed at large $r$. As for $g_{{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 ($m_{{L}}/m_{{S}}=(St_{{L}}/St_{{S}})^{1.5}\approx 30$), small particles will always be attracted towards the large ones, giving rise to the rapid growth of $g_{{SL}}(r)$ at small separation ($r \leq \eta$). Again, the opposite-sign attraction decays with increasing $r$, so $g_{{SL}}(r)$ approaches one when $r$ 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. For a pair of particles $i$ and $j$ with the separation $r = |\boldsymbol {x}_i-\boldsymbol {x}_j|$, the radial component of the relative velocity is defined as $w_{{r},ij} = (\boldsymbol {v}_i-\boldsymbol {v}_j) \boldsymbol {\cdot } (\boldsymbol {x}_i-\boldsymbol {x}_j)/|\boldsymbol {x}_i-\boldsymbol {x}_j|$. Taking the ensemble average over SS/LL/SL particle pairs then yields the mean relative velocities between three kinds of particle pairs as $\langle |w_{{r,SS}}| \rangle$, $\langle |w_{r,LL}| \rangle$ and $\langle |w_{{r,SL}}| \rangle$, respectively.

Figure 3(a) compares the mean radial relative velocity between SS pairs, $\langle |w_{{r,SS}}| \rangle$, for various particle charge $q$. For neutral SS pairs, when $r$ is in the inertial range ($\eta \ll r \ll L_{{int}}$), $\tau _{{p,S}}$ is much smaller than the characteristic time scales of turbulent fluctuations at this length scale. The radial relative velocity between SS pairs, $\langle |w_{{r,SS}}| \rangle$, thus follows that of fluid tracers $\delta u_{{r}}$. Here, the fluid relative velocity is obtained from the second-order longitudinal structure function as $\delta u_{{r}} = {C_2}^{1/2} \epsilon ^{1/3} r^{1/3}$ and the constant is $C_2=2.13$ (Yeung & Zhou Reference Yeung and Zhou1997) (shown as dash-dotted lines in figure 3). If $r$ is within the dissipative range, the time scale of local fluctuations becomes comparable to $\tau _{{p,S}}$. Small neutral particles cannot perfectly follow the background flow at this separation, so $\langle |w_{{r,SS}}| \rangle$ deviates from the fluid relative velocity $\delta u_{{r}}= (\epsilon /15 \nu )^{1/2} r$ (dashed lines in figure 3).

Figure 3.

Mean radial relative velocity $\langle |w_{{r}}| \rangle$ normalized by $u_{\eta }$ between (a) SS, (b) LL and (c) SL particle pairs with different particle charge $q$. The inset in (a) shows the full scale of$\langle |w_{{r,SS}}| \rangle$ in the neutral case. The red dashed line denotes the velocity scaling ${\propto }r$ in the dissipation range, and the red dash-dotted line denotes the velocity scaling ${\propto }r^{1/3}$ in the inertial range. Dashed lines from light to dark correspond to predictions using (3.15) (or its counterpart for the oppositely charged case) for $q=5 \times 10^{-15}\ \mathrm {C}$, $1 \times 10^{-14}\ \mathrm {C}$ and $2 \times 10^{-14}\, \mathrm {C}$, respectively.

When particles are charged, since the influence of Coulomb force is negligible at large $r$, curves for different $q$ collapse. In contrast, $\langle |w_{{r,SS}}| \rangle$ is found to rise significantly as particles get close to each other. When $q$ gets larger, the increase of $\langle |w_{{r,SS}}| \rangle$ becomes more pronounced and the deviation from the neutral result also occurs at a larger $r$. 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 $q$ is large, the strong Coulomb repulsion causes biased sampling at small $r$: 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 at a high velocity as the electrostatic force pushes them apart. As a result, pairs with a large $|w_{{r,SS}}|$ take a large proportion of the close pairs as $q$ increases, so the mean radial relative velocity between both approaching and departing pairs, $\langle w_{{r},ij} \rangle$, shows the increasing trend in figure 3(a). This increasing trend with $q$ will be further discussed in § 3.3.

Compared with small particles with $St=1$, large particles with $St=10$ are very inertial and insensitive to fluctuations even at large separations. The curves of $\langle |w_{{r,LL}}| \rangle$ are therefore much flatter. Besides, the constant relative velocity between LL pairs at small $r$ is also significantly larger, because their relative velocity comes from the energetic large-scale motions. When particles are charged, the same increasing trend with $q$ is also observed but is less significant, which can be attributed to the large kinetic energy compared with the electrostatic potential energy.

As for SL pairs, $\langle |w_{{r,SL}}| \rangle$ is larger than both $\langle |w_{{r,SS}}| \rangle$ and $\langle |w_{{r,LL}}| \rangle$ in the neutral cases. This is caused by the different responses of small/large particles to turbulent fluctuations, which is termed the differential inertia effect (Zhou et al. Reference Zhou, Wexler and Wang2001). Interestingly, there is no obvious change when comparing curves between different charges $q$ in figure 3(c). As shown in figure 2(c), the 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 large relative velocities over time because their response to local fluctuations is very different. Therefore, the influence of charge on $\langle |w_{{r,SL}}| \rangle$ is weak.

3.2. Effect of Coulomb force on particle flux

Once the RDF 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 Reference Sundaram and Collins1997; Wang et al. Reference Wang, Wexler and Zhou2000)

(3.2)\begin{equation} \varGamma_{ij}(R_{{C}}) = 2 {\rm \pi}R_{{C}}^2 \cdot g_{ij}(R_{{C}}) \cdot \langle |w_{{r},ij}|\rangle (R_{{C}}). \end{equation}

Here, $R_{{C}}=r_{i}+r_{j}$ is the collision radius, which equals the sum of the radii of particles $i$ and $j$; $g_{ij}(R_{{C}})$ is the RDF at $r=R_{{C}}$, and the mean relative velocity is $\langle |w_{{r},ij}|\rangle = \int _{-\infty }^{\infty } |w_{{r},ij}| p_{ij}(w_{{r},ij})\, \mathrm {d}w_{{r},ij}$, with $p_{ij}(w_{{r},ij})$ the probability density function (PDF) of $w_{{r},ij}$ at $r=R_{{C}}$.

Even though other particle parameters (e.g. the Stokes number) are set the same, different choices of $R_{{C}}$ will affect the final outcome of $\varGamma _{ij}$ by changing the collision geometry. Therefore, instead of directly using $\varGamma _{ij}$, we define the particle flux $\varPhi _{ij}$ as the ratio of the collision kernel $\varGamma _{ij}(r)$ to the area of the collision sphere $4 {\rm \pi}r^2$ at $r$ as

(3.3)\begin{equation} \varPhi_{ij}(r) = \frac{\varGamma_{ij}(r)}{4 {\rm \pi}r^2} = \frac{1}{2} g_{ij}(r) \cdot \langle |w_{{r},ij}|\rangle (r). \end{equation}

Here, $\varPhi _{ij}$ can be understood as the number of particles crossing the collision sphere per area per unit time, which is independent of the prescribed $R_{{C}}$.

Apart from (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

(3.4)\begin{equation} \varPhi_{ij}^{{in}}(r) = g_{ij}(r) F(w_{{r},ij}<0|r) \cdot \langle w_{{r},ij}^{-}\rangle (r). \end{equation}

Here, $F(w_{{r},ij}<0|r)$ is the fraction of particles moving inwards at $r$, and the mean inward relative velocity is $\langle w_{{r},ij}^{-}\rangle = - \int _{-\infty }^{0} w_{{r},ij} p_{ij}(w_{{r},ij})\, \mathrm {d}w_{{r},ij}$. After the system reaches equilibrium, RDFs between particle pairs no longer change, suggesting that the inward flux should balance the outward flux at all $r$. Here, the outward particle flux is

(3.5)\begin{equation} \varPhi_{ij}^{{out}}(r) = g_{ij}(r) F(w_{{r},ij}>0|r) \cdot \langle w_{{r},ij}^{+}\rangle (r), \end{equation}

with the mean outward relative velocity given by $\langle w_{{r},ij}^{+}\rangle = \int _{0}^{\infty } w_{{r},ij} p_{ij}(w_{{r},ij})\, \mathrm {d}w_{{r},ij}$. The definition of $\varPhi _{ij}$ is based on a pair of particles $i$ and $j$ without specifying their sizes. If the average is taken over all SS pairs, the result is the SS particle flux denoted by $\varPhi _{{SS}}$. The flux between LL pairs ($\varPhi _{{LL}}$) and SL pairs ($\varPhi _{{SL}}$) can also be obtained by taking the average over corresponding particle pairs.

The SS fluxes $\varPhi _{{SS}}$ defined by (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 (3.3) uses the information of both approaching and departing particle pairs, it is adopted in the following sections for better statistics.

Figure 4(a) compares the SS fluxes with various $q$. The flux for neutral pairs remain constant when $r$ is comparable to $\eta$. This indicates that, although both $g_{{SS}}$ (figure 2a) and $\langle |w_{{r,SS}}| \rangle$ (figure 3a) keep changing as $r$ decreases, their product is almost constant for small $r$. For inertial particle pairs with small $r$, it has been shown that $g(r) \propto r^{D_2-3}$, with $D_2$ the correlation dimension (Bec et al. Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007), while the relative velocity follows $\langle w_{{r,SS}} \rangle \propto r^{3-D_2}$ (Gustavsson & Mehlig Reference Gustavsson and Mehlig2014). Therefore, the dependence of $\varPhi _{{SS}}$ on $r$ is cancelled out. Besides, in practical situations the collision diameter $R_{{C}}$ between micron particles/droplets is smaller than $\eta$, so the constant flux at such separation distance leads to a quadratic dependence of collision kernel on $R_{{C}}$ as $\varGamma _{SS}=4 {\rm \pi}R_{{C,SS}}^2 \varPhi _{{SS}}(R_{{C,SS}})$ (Sundaram & Collins Reference Sundaram and Collins1997).

Figure 4.

Particle fluxes between (a) SS, (b) LL and (c) SL particle pairs with different particle charge $q$.

When particles are charged, the SS flux is found to decrease rapidly as $r$ drops. To characterize the effect of Coulomb force on the reduction of $\varPhi _{{SS}}$, we need to quantify the competition between the driving force and the resistance of particle collisions. For a pair of same-sign particles $i$ and $j$ separated by $r$, the driving force can be evaluated by the relative kinetic energy $E_{{Kin}}(r)=m_{ij} \langle |w_{{r},ij}| \rangle ^2/2$. Here, $m_{ij}=(m_i^{-1}+m_j^{-1})^{-1}$ is the effective mass, $\langle |w_{{r},ij}| \rangle (r)$ is the mean radial relative velocity between neutral pairs with a separation of $r$. The resistance is the electrical energy barrier $E_{{Coul}} = q_i q_j /4 {\rm \pi}\varepsilon _0 r$. The energy ratio is then given as

(3.6)\begin{equation} {Ct}_0 \equiv \frac{E_{{Coul}}}{E_{{Kin}}} = \frac{q_i q_j}{2 {\rm \pi}\varepsilon_0 r m_{ij} \langle |w_{{r},ij}| \rangle^2}. \end{equation}

Note that, in previous studies regarding the clustering of charged particles with weak inertia (Lu et al. Reference Lu, Nordsiek, Saw and Shaw2010a,Reference Lu, Nordsiek and Shawb), the approaching velocity between particle pairs can be directly derived through perturbation expansion in the Stokes number (Chun et al. Reference Chun, Koch, Rani, Ahluwalia and Collins2005). In this work, however, particles are very inertial, so the mean relative velocity $\langle |w_{{r},ij}| \rangle$ between neutral pairs is used instead. By using the corresponding effective mass (e.g. $m_{{SS}}=m_{{S}}/2$, $m_{{LL}}=m_{{L}}/2$, and $m_{{SL}}=m_{{S}}m_{{L}}/(m_{{S}}+m_{{L}})$) and the mean relative velocity ($\langle |w_{{r,SS}}|\rangle$, $\langle |w_{{r,LL}}|\rangle$, and $\langle |w_{{r,SL}}|\rangle$), (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 (3.6), the value of ${Ct}_0$ can be varied as the particle charges or the separation distance changes. Therefore, for each data point in figure 4(a), we plot the ratio of $\varPhi _{{SS}}(q)$ to its corresponding value in the neutral case $\varPhi _{{SS}}(q=0)$ as a function of ${Ct}_0$ in figure 5(a). The flux ratios for different $q$ 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 $\varPhi _{{LL}}$ in figure 4(b) are analysed in a similar way, and the results are plotted in figure 5(b). Because of the large kinematic energy between LL pairs, the Coulomb repulsion is relatively weak (${Ct}_0<1$), so the reduction of $\varPhi _{{LL}}$ has not yet reached the electrostatic-dominant regime. As for the flux between SL pairs shown in figure 4(c), the opposite trend is seen when plotting the flux ratio as a function of ${Ct}_0$ (figure 5c). The increase of $\varPhi _{{SL}}$ becomes significant when ${Ct}_0$ becomes larger than one.

Figure 5.

Ratios of charged particle flux $\varPhi (q)$ to neutral particle flux $\varPhi (q=0)$ as a function of the mean Coulomb-turbulence parameter ${Ct}_0$ for (a) SS, (b) LL and (c) SL particle pairs. Blue solid lines are model predictions from (3.16) and (3.17) with fitted $\beta$. Brown dash-dotted lines are predictions with the fixed $\beta =1$.

We now propose a model to quantify the influence of the electrostatic force on the mean particle flux. For particle pairs with a separation distance $r$, the mean radial relative velocities between different kinds of particle pairs have been shown in figure 3. We then assume that, for each kind of particle pair, the relative velocity between neutral pairs follows the Gaussian distribution. Taking SS particle pairs as an example, the PDF $p_{{r,SS}}$ of the relative velocity $w_{{r,SS}}(r)$ is

(3.7)\begin{equation} p_{{r,SS}}(w_{{r,SS}}) = \frac{1}{\sqrt{2 {\rm \pi}} \sigma_{{SS}}} \exp{\left(-\frac{w_{{r,SS}}^2}{2 \sigma_{{SS}}^2}\right)}. \end{equation}

The standard deviation is determined by the mean relative velocity at $r$ as $\sigma _{{SS}} = \sqrt {{\rm \pi} /2} \cdot \langle |w_{{r,SS}}| \rangle (r)$. It should be noted that the relative velocity distribution is actually non-Gaussian (Sundaram & Collins Reference Sundaram and Collins1997; Wang et al. Reference Wang, Wexler and Zhou2000; Ireland et al. Reference Ireland, Bragg and Collins2016a). However, a Gaussian distribution is sufficient for the first-order assumption, and has already been applied in previous studies (Pan & Padoan Reference Pan and Padoan2010; Lu & Shaw Reference Lu and Shaw2015).

We then evaluate the RDF of charged pairs. For neutral SS pairs, all particle pairs with an inward relative velocity ($w_{{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 ($w_{{r,SS}}<-w_{{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)\begin{equation} w_{{C,SS}}=\left( \frac{q^2}{2 {\rm \pi}\varepsilon_0 r m_{{SS}}}\right)^{1/2}. \end{equation}

The critical velocity for LL or SL pairs can be given by replacing $m_{{SS}}=m_{{S}}/2$ with $m_{{LL}}=m_{{L}}/2$ or $m_{{SL}}=m_{{S}}m_{{L}}/(m_{{S}}+m_{{L}})$. Note that the critical velocity is derived from the interaction between a pair of particles, which could reasonably describe the major effect of the 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. Reference Boutsikakis, Fede and Simonin2022, Reference Boutsikakis, Fede and Simonin2023). By assuming a sharp cutoff at $w=-w_{{C,SS}}$, Lu & Shaw (Reference Lu and Shaw2015) related the RDF of charged SS pairs $g_{{SS}}(r|q)$ to that of neutral SS pairs $g_{{SS}}(r|q=0)$ as

(3.9)\begin{align} g_{{SS}}(r|q) &= g_{{SS}}(r|q=0) \cdot \frac{\displaystyle\int_{-\infty}^{{-}w_{{C,SS}}} p_{{r,SS}}(w)\,\mathrm{d}w}{\displaystyle \int_{-\infty}^{0} p_{{r,SS}}(w)\,\mathrm{d}w} \nonumber\\ &=g_{{r,SS}}(r|q=0) \cdot \left[1-\mathrm{erf}\left(\frac{w_{{C,SS}}}{\sqrt{2} \sigma_{SS}}\right)\right], \end{align}

where $\mathrm {erf}$ is the error function. However, as the relative velocity magnitude $w_{{r}}$ becomes smaller than $w_{{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

(3.10)\begin{equation} \frac{g_{{SS}}(r|q)}{g_{{SS}}(r|q=0)} = \frac{\displaystyle\int_{-\infty}^{0} \varTheta_{{SS}}(w) p_{{r,SS}}(w)\,\mathrm{d}w}{ \displaystyle\int_{-\infty}^{0} p_{{r,SS}}(w)\,\mathrm{d}w}. \end{equation}

Here, $\varTheta _{{SS}}(w)$ is the electrical kernel that gradually transitions from one to zero when the magnitude of $w_{{r}}$ drops below $w_{{C,SS}}$. The proportion of charged particles that could overcome the Coulomb repulsion is then obtained from the convolution of $\varTheta _{{SS}}$ and $p_{{r,SS}}$, which is the numerator of the right-hand term in (3.10). To account for the smooth transition without adding significant complexity, we still adopt the sharp cutoff assumption, but the effective cutoff velocity $\beta _{{SS}}w_{{C,SS}}$ is used. Here, $\beta _{{SS}} \sim O(1)$ is the correction factor that ensures an accurate prediction as

(3.11)\begin{equation} \int_{-\infty}^{0} \varTheta_{{SS}}(w) p_{{r,SS}}(w)\,\mathrm{d}w \approx \int_{-\infty}^{-\beta_{{SS}}w_{{C,SS}}} p_{{r,SS}}(w)\,\mathrm{d}w. \end{equation}

The difference between the effective cutoff and the actual electrical kernel will be further discussed in § 3.3. Consequently, (3.9) becomes

(3.12)\begin{equation} g_{{SS}}(r|q) = g_{{SS}}(r|q=0) \cdot \left[1-\mathrm{erf}\left(\frac{\beta_{{SS}} \cdot w_{{C,SS}}}{\sqrt{2} \sigma_{SS}}\right)\right]. \end{equation}

To obtain the mean relative velocity between charged SS pairs, we start by computing the mean relative velocity in the velocity interval $[-\infty,-\beta _{{SS}} w_{{C,SS}}]$ between neutral SS pairs

(3.13)\begin{align} \langle |w_{{r,SS}}^{\prime}| \rangle (r|q=0) &= \frac{\displaystyle-\int_{-\infty}^{-\beta_{{SS}} w_{{C,SS}}} w \cdot p_{{r,SS}}(w)\,\mathrm{d}w}{ \displaystyle\int_{-\infty}^{-\beta_{{SS}} w_{{C,SS}}} p_{{r,SS}}(w)\,\mathrm{d}w} \nonumber\\ &= \langle |w_{{r,SS}}| \rangle (r|q=0) \cdot \exp{\left(\frac{-\beta_{{SS}}^2 w_{{C,SS}}^2}{2 \sigma_{SS}^2}\right)\Bigg/ \left[1-\mathrm{erf}\left(\frac{\beta_{{SS}} \cdot w_{{C,SS}}}{\sqrt{2} \sigma_{SS}}\right)\right]}. \end{align}

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 $[-\infty,-\beta _{{SS}} w_{{C,SS}}]$

(3.14)\begin{equation} \frac{1}{2}m_{{SS}} {\langle |w_{{r,SS}}| \rangle^2 (r|q)} = \frac{1}{2}m_{{SS}} {\langle |w_{{r,SS}}^{\prime}| \rangle^2 (r|q=0)} - \frac{\beta_{{SS}}^2 q^2}{4 {\rm \pi}\varepsilon_0 r}. \end{equation}

Here, the correction factor $\beta _{{SS}}$ is also added to the last term on the right-hand side to account for the smooth transition. Taking (3.13) into (3.14) then yields

(3.15)\begin{align} \langle |w_{{r,SS}}| \rangle (r|q) &= \langle |w_{{r,SS}}| \rangle (r|q=0) \nonumber\\ &\quad \cdot \left[\frac{\displaystyle\exp{\left(-\dfrac{\beta_{{SS}}^2 w_{{C,SS}}^2 }{\sigma_{SS}^2}\right)}}{ \displaystyle\left[1-\mathrm{erf} \left(\dfrac{\beta_{{SS}}w_{{C,SS}}}{\sqrt{2}\sigma_{SS}}\right)\right]^2} - \frac{\beta_{{SS}}^2 q^2}{2 {\rm \pi}\varepsilon_0 m_{{SS}} \cdot \langle |w_{{r,SS}}| \rangle^2 (r|q=0) \cdot r}\right]^{1/2} . \end{align}

Finally, by multiplying (3.9) and (3.15) and taking into account the definition of ${Ct}_0$ (3.6), the flux ratio can be given as

(3.16)\begin{align} \frac{\varPhi_{{SS}}(q)}{\varPhi_{{SS}}(q=0)} &= \frac{g_{{SS}}(r|q) \cdot \langle |w_{{r,SS}}| \rangle (r|q)}{g_{{SS}}(r|q=0) \cdot \langle |w_{{r,SS}}| \rangle (r|q=0)} \nonumber\\ &= \left\{ \exp{\left(-\frac{2 \beta_{{SS}}^2 {Ct}_0}{\rm \pi}\right)} - \beta_{{SS}}^2 {Ct}_0 [1-\mathrm{erf}(\beta_{{SS}} \sqrt{{Ct}_0 / {\rm \pi}})]^2 \right\}^{1/2}. \end{align}

The flux ratios for LL pairs and SL pairs can be derived in a similar way and are written as

(3.17a)$$\begin{gather} \frac{\varPhi_{{LL}}(q)}{\varPhi_{{LL}}(q=0)} = \left\{\exp{\left(-\frac{2 \beta_{{LL}}^2 {Ct}_0}{\rm \pi}\right)} - \beta_{{LL}}^2 {Ct}_0 [1-\mathrm{erf}(\beta_{{LL}} \sqrt{{Ct}_0 / {\rm \pi}})]^2\right\}^{1/2}, \end{gather}$$

(3.17b)$$\begin{gather}\frac{\varPhi_{{SL}}(q)}{\varPhi_{{SL}}(q=0)} = \left\{ \exp{\left(-\frac{2 \beta_{{SL}}^2 {Ct}_0}{\rm \pi}\right)} + \beta_{{SL}}^2 {Ct}_0 [1+\mathrm{erf}(\beta_{{SL}} \sqrt{{Ct}_0 / {\rm \pi}})]^2\right\}^{1/2}. \end{gather}$$

Equations (3.16) and (3.17) are then fitted as blue lines in figure 5, which show good agreement with the simulation results. The fitted correction factors are $\beta _{{SS}}=0.193$, $\beta _{{LL}}=0.739$ and $\beta _{{SL}}=0.5415$, satisfying the assumption that the correction factors are of the order of unity. The predictions with the fixed $\beta =1$ are also shown as brown dash-dotted lines, which correspond to the original model that assumes the sharp cutoff occurs at the critical velocity $w_{{C,SS/LL/SL}}$. Although the trends are similar, models with the fixed $\beta =1$ underestimate the critical ${Ct}_0$ at which the transition occurs. Moreover, the proposed model underestimates the SS flux when ${Ct}_0 \approx 10^2$ in figure 5(a). To find the origin of this discrepancy, the predicted $g_{ij}(r)$ using (3.12) and the predicted $\langle |w_{{r,SS}}| \rangle (r)$ using (3.15) are also plotted as grey dashed lines in figures 2 and 3, respectively. It can be seen that the predicted RDFs are comparable to the simulation results, so the discrepancy of $\varPhi _{{SS}}$ at ${Ct}_0 \approx 10^2$ mainly comes from the underestimation of the mean relative velocity at $r \sim \eta$ in figure 3(a). Since the PDF of $w_{{r}}$ is assumed Gaussian in our model, the influence of intermittency is not considered. As a result, the proportion of particle pairs with large relative velocity is significantly underestimated.

3.3. Particle flux density for charged particles

In § 3.2, the mean Coulomb-turbulence parameter ${Ct}_0$ is defined using the mean relative velocity $\langle |w_{{r},ij}| \rangle$, which compares the importance of the Coulomb force with 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 interest to examine the impacts of the Coulomb force on particle pairs with different relative velocities. For different kinds of particle pairs, the particle flux $\varPhi _{ij}$ in the relative velocity space is expanded as

(3.18)\begin{equation} \varPhi_{ij} = \int_{-\infty}^{\infty} \phi_{ij}(w_{{r}})\,\mathrm{d} w_{{r}}. \end{equation}

Here, $\phi _{ij}$ is the density of particle flux within each relative velocity interval, which measures the contribution to the total particle flux $\varPhi _{ij}$ by particle pairs whose relative velocity is within the interval $w_{{r}} \pm \Delta w_{{r}}/2$. By comparing (3.18) and (3.3), $\phi _{ij}$ can be given as

(3.19)\begin{equation} \phi_{ij}=\tfrac{1}{2} g_{ij} \cdot |w_{{r},ij}|p_{ij}(w_{{r},ij}). \end{equation}

We start with the effect of the Coulomb force on the PDFs of relative velocity $p_{ij}(w_{{r},ij})$, because the distribution of $\phi _{ij}$ in (3.19) is strongly dependent on $p_{ij}(w_{{r},ij})$. Figure 6(a) illustrates PDFs of $w_{{r,SS}}$ at the separation distance interval $r \in [1.15\eta, 2.5\eta )$. The Gaussian distribution defined by (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 $|w_{{r}}|<3u_\eta$, but significantly underestimates the probability of large $|w_{{r}}|$. Therefore, the Gaussian assumption is only suitable for modelling 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.

Figure 6.

The PDFs of the relative velocity between (a) SS, (b) LL and (c) SL pairs at $r \in [ 1.15\eta, 2.5\eta )$. Zoom-in views around $w_{{r}}=0$ are shown in the insets. Black dashed lines denote the Gaussian distributions defined by (3.7) for different kinds of particle pairs.

For SS pairs, compared with the neutral PDF, charged PDFs become wider as $q$ 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 $q$ in figure 3(a). If we look at the symmetry of $p_{{SS}}$, the neutral PDF is found negatively skewed. This asymmetry is attributed to: (i) the fluid velocity derivative is negatively skewed, and the relative motion between particle pairs with $St=1$ is partially coupled to the local flow and thus exhibits a similar feature (Van Atta & Antonia Reference Van Atta and Antonia1980; Wang et al. Reference Wang, Wexler and Zhou2000); (ii) the asymmetry of particle path history gives rise to a larger relative velocity between approaching pairs than departing pairs (Bragg & Collins Reference Bragg and Collins2014). However, the asymmetry of $p_{{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 the Coulomb force relies only on the interparticle distance $r$, so the approaching or departing SS pairs experience the same amount of repulsion as long as $r$ 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 $w_{{r,SS}}$.

Table 2.

Normalized skewness of the radial relative velocity $w_{{r}}$ at $r \in [ 1.15\eta, 2.5\eta )$.

For LL pairs, since Coulomb repulsion is only able to repel pairs with small relative velocity, $p_{{LL}}$ drops slightly as $w_{{r}}$ approaches zero, but remains almost unchanged at larger $w_{{r}}$. As for $p_{{SL}}$, the difference between the 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 the Coulomb force is limited. Besides, the movement of large particles with $St=10$ has very weak couplings to the local flow fields, so the skewness of $p_{{LL}}$ and $p_{{SL}}$ in table 2 is negligible.

We now turn to the distribution of flux density $\phi _{ij}$ in the velocity space. Figure 7(a,c,e) plots the flux densities between SS, LL and SL pairs, respectively. Different from the unimodal PDFs of the relative velocity shown in figure 6, $\phi _{ij}$ are bimodal because the maximum flux density is reached when the product $|w_{{r},ij}| \cdot p_{ij}(w_{{r},ij})$ is the largest. Therefore, it is the particle pairs with the intermediate relative velocity that contribute the most to the overall particle flux.

Figure 7.

Particle flux density in the velocity space, $\phi (w_{{r}})$, between (a) SS, (c) LL and (e) SL pairs at $r \in [1.15\eta, 2.5\eta )$. Zoom-in views around $w_{{r}}=0$ are shown in the insets. (b,d,f) Dependence of Coulomb kernel $\varTheta _{ij}=\phi _{ij}(q)/\phi _{ij}(q=0)$ on $w_{{r},ij}$ corresponding to (a,c,e), respectively.

To better describe the impacts of the Coulomb force, we define the Coulomb kernel $\varTheta _{ij}$ as the ratios of the charged flux densities to their neutral values, i.e. $\varTheta _{ij}(q) = \phi _{ij}(q)/\phi _{ij}(q=0)$, which are displayed in the right panel (i.e. (b), (d) and (f)) of figure 7. As shown in figure 7(b), the value of $\varTheta _{{SS}}$ varies by more than one order of magnitude, indicating that the influence of the Coulomb force changes significantly as $w_{{r}}$ changes. When $|w_{{r}}|$ is small, $\varTheta _{{SS}}$ decreases drastically as $q$ increases because Coulomb repulsion dominates. In addition, $\varTheta _{{SS}}$ is found to be asymmetric. As discussed above, neutral particles with $St=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 $w_{{r}}$ to a larger positive $w_{{r}}$ in the velocity space. Consequently, $\varTheta _{{SS}}$ experiences a more significant decrease at $w_{{r}}$ slightly greater than zero, but quickly exceeds one when $w_{{r}} \approx 2 u_\eta$. Moreover, $\varTheta _{{SS}}$ becomes larger than one at large $|w_{{r}}|$ for both approaching and departing pairs. The explanation is that, with the increase of $q$, 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 $|w_{{r}}|$. However, the contribution of the increased flux at large positive $|w_{{r}}|$ is negligible compared with the decrease at small $|w_{{r}}|$ (see figure 7a), so the major effect of the Coulomb force is to reduce the SS flux through the SS repulsive force.

The kernel $\varTheta _{{LL}}$ between LL pairs (figure 7d) also drops quickly at small $|w_{{r}}|$, but recovers to unity when $|w_{{r}}| \geq 2u_\eta$. Besides, because of the limited effects of the local fluid velocity gradient mentioned above, the approaching and departing processes are more symmetric for neutral LL pairs, and adding the Coulomb force retains this symmetry. As for $\varTheta _{{SL}}$ shown in figure 7(f), the change is still most significant near $w_{{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. $|w_{{r}}|/u_\eta \leq 2$), the opposite-sign attraction seems to increase $\varTheta _{{SL}}$ in a much wider range of $w_{{r}}$. For instance, $\varTheta _{{SL}}$ is increased even at $w_{{r}}/u_\eta =-10$ for $q=2\times 10^{-14}\ \mathrm {C}$ in figure 7(f). As discussed in § 3.2, the main effect of the opposite-sign attraction is to enhance SL clustering. Then, particles of different sizes, although 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 $\varTheta _{{SL}}$ in a wide range of $w_{{r}}$.

The discussion above has shown that the influence of the Coulomb force is different as the particle relative velocity $w_{{r},ij}$ changes. Therefore, instead of using the mean relative velocity $\langle |w_{{r},ij}| \rangle$ (3.6), we adopt the median relative velocity $\bar {w}_{{r},ij}$ in each $w_{{r}}$ bin to estimate the relative kinetic energy. For a given separation distance $r$ and a certain relative velocity bin with the median $\bar {w}_{{r},ij}$, the extended Coulomb-turbulence parameter is given as

(3.20)\begin{equation} {Ct} \equiv \frac{q_i q_j}{2 {\rm \pi}\varepsilon_0 r m_{ij} \bar{w}_{{r},ij}^2}. \end{equation}

We then examine the dependence of the electrical kernel $\varTheta _{ij}$ on the extended Coulomb-turbulence parameter ${Ct}$. For the data points $(w_{{r},ij},\varTheta _{ij})$ shown in the right panel (i.e. (b), (d) and (f)) of figure 7, the corresponding values of ${Ct}$ are calculated using (3.20), and the results are re-plotted as points $({Ct}, \varTheta _{ij})$. Note that to ensure meaningful statistics, only data points satisfying a certain criterion are considered and the reasons are as follows:

1. (i) In addition to the separation interval $r \in [1.15\eta,2.5\eta )$ shown in figure 7, data from two more separation intervals, i.e. $[0.85\eta,1.15\eta )$ and $[2.5\eta,3.5\eta )$, are also added. At these separation distances, the effect of the electrostatic force is already significant, while the number of samples is large enough for statistics.

2. (ii) We only consider pairs with their relative velocity in a certain range to avoid using the more scattered data points at large $w_{{r}}$. For SS pairs, the relative velocity range is $[-5 u_\eta, 5 u_\eta ]$, while for LL/SL pairs the range considered is $[-10 u_\eta, 10 u_\eta ]$. The particle flux within the above ranges contributes 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 $\varPhi _{ij}$ in each case.

3. (iii) As shown in figure 7(b), $\varTheta _{{SS}}$ is not symmetric about $w_{{r}}=0$. When particles are departing, the Coulomb repulsion will redistribute $\phi _{{SS}}$ in the velocity space, which may distort the $\varTheta _{{SS}}-{Ct}$ relationship. We therefore only use the data from approaching pairs ($w_{{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 $\varPhi _{{SS}}(q)=2\int _{-\infty }^{0} \varTheta _{{SS}}({Ct}) \phi _{{SS}}(w_{{r,SS}})\, \mathrm {d} w_{{r,SS}}$.

Figure 8 plots the dependence of $\varTheta _{ij}$ on ${Ct}$ for different particle pairs. Despite different particle charge $q$ and separation $r$, the data points for both SS (figure 8a) and LL (figure 8b) pairs show a similar trend. When ${Ct}$ is small, the Coulomb force is weak compared with particle–turbulence interaction, so $\varTheta _{{SS}}$ and $\varTheta _{{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 $\varTheta$ is observed. In § 3.2, it has been assumed that a sharp cutoff occurs at the effective velocity $\beta _{{SS}} w_{{C,SS}}$ or $\beta _{{LL}} w_{{C,LL}}$ for SS or LL pairs. However, the transitions of $\varTheta _{{SS}}$ and $\varTheta _{{LL}}$ in figure 8(a) turn out to be smooth. Thus, the sharp cutoff can be understood as the first-order approximation of the influence of Coulomb repulsion, which can be written as

(3.21)\begin{equation} \varTheta_{ij}({Ct}) = 1-H({Ct}_{{crit},ij}), \end{equation}

with $H(\,\cdot\,)$ being the Heaviside step function. The critical value of ${Ct}$ describing where $\varTheta _{ij}$ steps down can be related to the fitted correction factors for the effective cutoff velocity as

(3.22)\begin{equation} {Ct}_{{crit},ij}=1/\beta_{ij}^2. \end{equation}

By taking into the fitted results ($\beta _{{SS}}=0.193$ and $\beta _{{LL}}=0.739$) from § 3.2, the critical values can be given as ${Ct}_{{crit,SS}} = 26.85$ and ${Ct}_{{crit,LL}}=1.83$, respectively.

Figure 8.

Dependence of the electrical kernel $\varTheta _{ij}$ on the extended Coulomb-turbulence parameter ${Ct}$ for (a) SS, (b) LL and (c) SL particle pairs. Different shapes represent data points from $r \in [0.85\eta,1.15\eta )$ (circles: $\circ$), $r \in [1.15\eta,2.5\eta )$ (diamonds: $\diamond$) and $r \in [2.5\eta,3.5\eta )$ (triangles: $\triangle$), respectively. Open/filled symbols denote data points for approaching/departing pairs.

However, different from $\varTheta _{{SS}}$ and $\varTheta _{{LL}}$, there is no clear dependence of $\varTheta _{{SL}}$ on ${Ct}$ in figure 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 $E_{{Kin}}$ becomes independent of the electrical potential energy $E_{{Coul}}$, and no clear relationship can be found between $\varTheta _{{SL}}$ and $\phi _{{SL}}$.