1. Introduction
The preferential concentration of solid particles in wall-bounded turbulent flows is of primary importance in a wide range of technological applications. For this reason, it has been the subject of numerous experimental and numerical studies over the years. In general, particle accumulation results from the complex interplay between the forces acting on the particles: aerodynamic, inertial, collisional and gravitational forces. In practice, however, in confined flows particles usually acquire a certain amount of electrostatic charge through triboelectric effects, i.e. via collisions of particles with the walls or with other pre-charged particles. In some applications, such as powder coating (Yang et al. Reference Yang, Ma, Zhu, Chow and Shi2016), electrostatic precipitation (Horender, Schaub & Sommerfeld Reference Horender, Schaub and Sommerfeld2014), triboelectric separation (Labair et al. Reference Labair, Touhami, Tilmatine, Hadjeri, Medles and Dascalescu2017) and spray painting (Böttner & Sommerfeld Reference Böttner and Sommerfeld2002), particle charging takes place deliberately to control the particles’ trajectories. However, in some other powder-handling facilities, such as fluidized beds, it occurs incidentally and leads to operational problems (Fotovat, Bi & Grace Reference Fotovat, Bi and Grace2017, Reference Fotovat, Bi and Grace2018). For example, due to particle–wall adhesion the charged particles can coat vessel walls, which requires frequent cleaning (Sippola et al. Reference Sippola, Kolehmainen, Ozel, Liu, Saarenrinna and Sundaresan2018). Also, the dynamics of fluidized beds is greatly influenced by electrostatic interactions (Jalalinejad, Bi & Grace Reference Jalalinejad, Bi and Grace2012, Reference Jalalinejad, Bi and Grace2015). In pharmaceutical devices, e.g. dry powder inhalers, electrostatic effects deteriorate the effectiveness of the final product (Wong, Kwok & Chan Reference Wong, Kwok and Chan2015). During pneumatic conveying, where the emerging charge levels are particularly high (Klinzing Reference Klinzing2018), high concentrations of charged particles may lead to spark discharges, with devastating consequences. In the past, spark discharges have caused numerous dust explosions in chemical and process industries (Eckhoff Reference Eckhoff2003). Thus far, the effect of triboelectric charging and the ensuing electrostatic forces on particle concentrations has not been investigated in detail. Therefore, many questions on this topic still remain open.
Uncharged solid particles are well known to distribute non-uniformly in wall-bounded turbulent flows. This is related to their turbophoretic drift, i.e. the tendency of particles to migrate towards regions of reduced turbulent kinetic energy (Caporaloni et al. Reference Caporaloni, Tampieri, Trombetti and Vittori1975; Reeks Reference Reeks1983). The design of advanced flow solvers and the ever-increasing power of modern computer has made possible to study the phenomenon of turbophoresis via direct numerical simulations (DNS). For example, one of the earliest DNS on this topic was that of McLaughlin (Reference McLaughlin1989). Turbophoresis was further explored by Marchioli & Soldati (Reference Marchioli and Soldati2002) who performed DNS of particle-laden flow over a smooth flat plate. Therein, they reported that, as the characteristic particle response time scale becomes shorter, the particles congregate faster and closer to the flat wall. In the same study, the authors also identified the coherent structures that are responsible for the entrapment of particles in the near-wall region. More recently, the study of Wang (Reference Wang2010) on turbophoresis in channel flows elucidated the preferential location for particles of different inertia and number density.
On the other hand, the properties of duct flows are quite different to those of flows over a flat plate or even channel flows. This is due to the absence of a homogeneous direction which results in the development of secondary flow structures. More specifically, as shown schematically in figure 1, duct flows are characterized by the formation of Prandtl's secondary flows of the second kind (Prandtl Reference Prandtl1927; Bradshaw Reference Bradshaw1987; Nezu Reference Nezu2005). These secondary flows are formed due to turbulence-driven gradients of the Reynolds stress and are related to the variations of the probability distribution of the coherent flow structures (Uhlmann et al. Reference Uhlmann, Pinelli, Kawahara and Sekimoto2007; Pinelli et al. Reference Pinelli, Uhlmann, Sekimoto and Kawahara2010; Kawahara, Uhlmann & VanVeen Reference Kawahara, Uhlmann and VanVeen2012).
Figure 1. Symmetry planes of the duct and streamlines of induced Prandtl's secondary flow of the second kind. The arrows indicate the direction of the possible particle motion in the cross-section following the gas phase.
Although these secondary flows involve only small-amplitude fluid motion, they have a significant effect on the overall momentum and mass transport. For example, in a duct with square cross-section, the mean secondary flow follows the well-known 8-vortex pattern, as can be observed in figure 1. The intensity of the spanwise motion inside these vortices is between 1 % and 2 % of the bulk fluid velocity (Pirozzoli et al. Reference Pirozzoli, Modesti, Orlandi and Grasso2018). Nonetheless, the role of these vortices is important because they transport high-momentum fluid from the duct core towards the corners.
It is also important to mention that the turbulence modelling of such secondary flows is particularly challenging because the eddy viscosity can no longer be assumed to be isotropic. Standard turbulence models, however, are based on the assumption of isotropy and, therefore, cannot satisfactorily handle secondary flows induced by Reynolds-stress gradients (Speziale Reference Speziale1982; Mani et al. Reference Mani, Babcock, Winkler and Spalart2013).
In the case of particle-laden flows, these secondary motions of the second type have a strong effect on the local concentration of particles. Actually, particle-laden flows in ducts have been the subject of several numerical studies that have appeared in the literature, albeit they are not as much investigated as single-phase duct flows. For example, Wang, Zhao & Yao (Reference Wang, Zhao and Yao2019) studied an open duct flow confined by three walls of the Reynolds number 
$Re=83\,000$. Therein, it was demonstrated that relatively light particles, of Stokes number 
$St$ between 3.2 and 96.6, tend to move diagonally in the direction indicated by the arrow A in figure 1 and to accumulate in the corners. Wang et al. (Reference Wang, Zhao and Yao2019) also showed that, on the other hand, heavy particles (
$St>323$) are much less affected by secondary flows and are more likely to concentrate in the bulk of the duct. Further, Yao & Fairweather (Reference Yao and Fairweather2010) performed large-eddy simulations (LES) to study particle resuspension in a duct flow at 
$Re=250\,000$. In that study, the authors demonstrated that secondary flows enable the displacement of particles from the wall region back into the bulk of the duct, which has a significant impact on the local particle concentration. In the same study, it was also explained that lift forces counteract gravitation and support the transport of particles parallel to the horizontal walls and in the direction pointed by the arrow B in figure 1.
In another study, Sharma & Phares (Reference Sharma and Phares2006) highlighted the importance of secondary flows for particle dispersion in ducts with square cross-sections. They revealed that passive tracers and low-inertia particles are subject to a lateral advective transport that is absent in pipe and channel flows. In this manner, the particles tend to circulate between the core and the boundary of the duct, i.e. following arrows A, B and C that are shown in figure 1. On the other hand, high-inertia particles follow mainly follow the directions of arrows A and B and tend to accumulate close to the wall.
With regard to electrostatics, as mentioned above, particles usually carry a certain amount of electric charge. The emerging electric field significantly modifies the flow patterns of the particulate phase due to electrostatic forces between particles (Dhodapkar Reference Dhodapkar1991). In some cases, these forces have a counterintuitive effect; for example, they may cause particles to move upstream and against the fluid flow (Myler Reference Myler1987). At present, numerical studies on triboelectric charging in particle-laden duct flows are particularly scarce. This may be attributed to the complexity of the system of governing equations which consists of the Navier–Stokes equations for the fluid, the electric field equation, plus the equations of motion of the particulate phase suitably modified so as to take into account the particle–fluid and particle–electric field interactions. As a result, our understanding of this type of flows remains rudimentary.
Recently, Yao, Li & Zhao (Reference Yao, Li and Zhao2020) studied via LES the effect of the electrostatic field on a particle-laden pipe flow at 
$Re=44\,000$, by assuming one-way coupling between particles and carrier fluid. According to these simulations, the electrostatic forces on particles dominate over aerodynamic forces and gravity in the near-wall region. This resulted in higher particle concentrations close to the wall and lower ones in the bulk of the pipe, the latter one been defined as up to 
$95\,\%$ of the pipe's radius. DNS of triboelectric powder charging were performed by Grosshans & Papalexandris (Reference Grosshans and Papalexandris2017a); however, that study treated channel flows and was limited to the initial states of charging. For this reason, the electrostatic charges had a minor impact on the flow patterns.
To our knowledge, the only available studies on the effect of electrostatic forces on particle-laden duct flows are the LES reported by Grosshans (Reference Grosshans2018) and Grosshans et al. (Reference Grosshans, Villafañe, Banko and Papalexandris2019). Those studies, which considered flows at 
$Re=10\,000$ and high-inertia particles, revealed that increasing the electrostatic charge results in (i) higher particle concentration in the corners of the duct, and (ii) more uniform particle concentrations away from them. However, the underlying physical mechanism of this behaviour was not explored in those simulations. It is also worth mentioning that reciprocal investigations, i.e. studies on the influence of different flow parameters (such as 
$Re$ or 
$St$) on the electrostatic field in duct flows are also currently unavailable.
In summary, the two-way interaction between electrostatic field and particle-laden flows in ducts remains largely unexplored. The present work aimed at providing new insight on this interaction by means of DNS. To this end, we employed our newly developed computational tool pafiX (2019); its name stands for particle flow simulation in explosion protection. This tool is available as freeware and is developed specifically for the study of powder and fluid flows under conditions involving a hazard to operational safety. In our study, we considered particle-laden flows at two different density ratios, 
$\rho _{p}/\rho =1000$ and 7500, and two different frictional Reynolds numbers, 
$Re_{\tau }=300$ and 600.
The paper is structured as follows. Section 2 provides an overview of the governing equations and the algorithms that are implemented in pafiX for their numerical treatment. In § 3 we first provide the validation tests of these algorithms and subsequently we present and discuss the results of our DNS study. Finally, § 4 concludes.
2. Mathematical model and numerical method
In this section we outline the mathematical model and the numerical method for its solution. The system of governing equations consists of three strongly coupled parts, namely (i) the Navier–Stokes equations that describe the flow of the carrier gas, (ii) Gauss's law for the electrostatic field and (iii) Newton's law of motion of the solid particles.
For a constant-density, particle-laden flow, the Navier–Stokes equations read
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} {\boldsymbol{u}} =0, \end{gather}
\begin{gather} \frac{\partial {\boldsymbol{u}}}{\partial t} + ({\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla}) {\boldsymbol{u}} = - \frac{1}{\rho} \boldsymbol{\nabla} p + \nu \nabla^2 {\boldsymbol{u}} + {\boldsymbol{F}}_{s} + {\boldsymbol{F}}_{f}. \end{gather}
Using standard notation, 
$\boldsymbol {u}$, 
$\rho$ and 
$p$ stand for the fluid velocity, density and dynamic pressure respectively. Also, 
$\nu$ is the kinematic viscosity of the fluid. Further, 
${\boldsymbol {F}}_{s}$ is the source term accounting for the momentum transfer between the solid particles and the carrier gas. More specifically, its integral over a control volume, e.g. a computational cell, is equal to the opposite of the sum of the aerodynamic forces that act on the particles that are located inside the control volume. Finally, 
${\boldsymbol {F}}_{f}$ represents an externally applied forcing that balances the momentum loss of the fluid due to wall friction.
As is well known, the dynamic pressure enters the constant-density Navier–Stokes equations via its gradient and, therefore, its role is to constrain the fluid motion so that mass continuity (2.1a) is satisfied. In other words, the role of the dynamic pressure is to limit the fluid to isochoric motions. In the proposed algorithm, this constraint is satisfied by coupling the continuity and momentum equations through the distributed Gauss–Seidel scheme originally proposed by Brandt & Dinar (Reference Brandt and Dinar1978). Herein, this algorithm is extended to three-dimensional domains and non-uniform grids. The underlying idea of this scheme is to diminish the error in the continuity equation by iteratively adjusting the velocity field. Afterward, the pressure field is modified accordingly so that the residuals of the momentum equation at all points remain unchanged. Distributive relaxation represents an efficient and intuitive alternative to popular velocity–pressure schemes such as semi-implicit method for pressure-linked equations (SIMPLE) or pressure-implicit with splitting of operators (PISO) (Ferziger & Peric Reference Ferziger and Peric2002).
As mentioned previously, the electric field strength 
${\boldsymbol {E}}$ follows Gauss's law. In the electrostatic approximation, 
${\boldsymbol {E}}$ is defined in terms of the electric potential 
$\varphi _{el}$,
\begin{equation} {\boldsymbol{E}}= -\boldsymbol{\nabla} \varphi_{el}, \end{equation}so that Gauss's law reduces to
\begin{equation} \nabla^2 \varphi_{el} = -\dfrac{\rho_{{el}}}{\varepsilon}, \end{equation}
which is a Poisson equation. In this equation, 
${\rho _{{el}}}$ is the electric charge density. If no external electric field is present, the electric charge density results directly from the positions of the individual particles and their individual charge. The electric permittivity of the solid–gas mixture, 
${\varepsilon }$, can be approximated by the value of the free space either if the permittivity of the particles is low or if the solid volume fraction is very small (Rivas & Iglesias Reference Rivas and Iglesias2007; Rokkam, Fox & Muhle Reference Rokkam, Fox and Muhle2010) which is the case in our simulations. Thus, we apply a value of 
$\varepsilon =8.85 \times 10^{-12}\ \textrm {F}\ \textrm {m}^{-1}$.
The numerical solution of the above equations requires the discretization of the spatial derivatives of 
$\boldsymbol {u}$, 
$p$ and 
$\varphi _{el}$. The convective term of the Navier–Stokes equations is approximated by a fifth-order accurate weighted essentially non-oscillatory scheme (Jang & Shu Reference Jang and Shu1996). The pressure gradient and viscous terms in (2.1b), as well as the velocity derivatives of (2.1a), are discretized via fourth-order central differences. Further, the left-hand side of Gauss's law (2.3) is discretized via second-order central differences. Then the discretization of the Poisson equation results in a linear system that can be solved by standard linear solvers.
In order to keep the high-order discretization schemes simple, the spatial derivatives of the convective terms are transformed via stretching of the spatial coordinates. For a generic flow quantity 
$\phi$, this transformation is performed with the application of the chain rule as follows,
\begin{equation} \dfrac{\textrm{d}\phi}{\textrm{d}x} = \dfrac{1}{x'(\xi)} \dfrac{\textrm{d}\phi}{\textrm{d}\xi}, \end{equation}
where the prime symbol denotes the derivative of the 
$x$ variable with respect to the stretched 
$\xi$ variable. Thus, the non-uniform grid in the physical 
$x$ direction is mapped to the uniform grid in the 
$\xi$ direction on which the derivatives are then discretized.
Further, in order to reduce the required memory and computing time, the deferred-correction method is applied. According to it, the numerical approximation of a quantity 
$\phi$ is expressed by
\begin{equation} \phi = \phi^{{l}} + \left( \phi^{{h}} - \phi^{{l}} \right)^{{old}}, \end{equation}
where the superscript ‘
$l$’ denotes an approximation by a low-order scheme and ‘
$h$’ an approximation by a high-order scheme. The terms indicated by ‘old’ are computed explicitly using values from the previous outer iteration. For the computational cells in the vicinity of a solid boundary there are not a sufficient number of points available to form a five-point stencil for the high-order discretization of the spatial derivatives. In those computational cells, only the low-order schemes are retained.
Time integration of (2.1b) is performed via an implicit second-order scheme using a variable time step, i.e.
\begin{equation} \dfrac{\partial {\boldsymbol{u}}}{\partial t} \approx \dfrac{ \left( 1+\tau^{n+{1}/{2}} \right) {\boldsymbol{u}}^{n+1} - \left( 1+\tau^{n+{1}/{2}}+\tau^{n-{1}/{2}} \right) {\boldsymbol{u}}^{n} + \left( \tau^{n-{1}/{2}} \right) {\boldsymbol{u}}^{n-1} } {\tau^{n+{1}/{2}} \Delta t^{n+1}}, \end{equation}with
\begin{equation} \tau^{n+{1}/{2}} = \dfrac{\Delta t^{n+1}}{\Delta t^{n+1} + \Delta t^{n}} \quad \mathrm{and} \quad \tau^{n-{1}/{2}} = \dfrac{\Delta t^{n}}{\Delta t^{n} + \Delta t^{n-1}} . \end{equation}
In the above equations, the superscript 
$n$ is a label for the current time instance, 
$t^n$, i.e. the most advanced time at which the solution has been computed. Also, 
$\Delta t^{n+1}$ is the time step used to advance the solution from the 
$n$th time level to the next one. The time step is determined via the Courant–Friedrichs–Lewy (CFL) condition. For example, in the simulations presented herein, the maximum Courant number was set equal to 0.2.
The above algorithm has been implemented in staggered grids which are not prone to the well-known odd–even decoupling between the pressure and velocity fields. More specifically, the variables 
$p$, 
${\boldsymbol {E}}$ and 
${\rho _{{el}}}$ are computed at the cell centres whereas the velocity vector 
${\boldsymbol {u}}$ and the source term 
${\boldsymbol {F}}_{s}$ are calculated at the centre of the cell faces. The resulting linear system of equations is solved using the Jacobi method which facilitates the straightforward parallelization of the code. At each time step, (2.1a), (2.1b) and (2.3) are relaxed until the 
$L_2$ norm of the error of each variable drops by three orders of magnitude from its initial value.
The position of the particles is calculated in the Lagrangian frame of reference. More specifically, Newton's second law of motion is solved separately for each particle,
\begin{equation} \dfrac{\mathrm{d} \boldsymbol{u}_{p}}{\mathrm{d} t} = {\boldsymbol{f}}_{{fl}} + {\boldsymbol{f}}_{{coll}} + {\boldsymbol{f}}_{{el}} + {\boldsymbol{f}}_{{g}} . \end{equation}
In this equation, 
${\boldsymbol {f}}_{{fl}}$ represents the acceleration due to the force exerted by the surrounding fluid on the particle (aerodynamic force), 
${\boldsymbol {f}}_{{coll}}$ represents the acceleration due to particle–particle or wall–particle collisions and 
${\boldsymbol {f}}_{{el}}$ is the acceleration due to electrostatic forces. Finally, 
${\boldsymbol {f}}_{{g}}$ stands for the gravitational acceleration.
The acceleration due to the aerodynamic force is computed by the following expression (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2012),
\begin{equation} {\boldsymbol{f}}_{fl} = - \frac{3 \rho}{8 \rho_{p} r} C_{d} \vert {\boldsymbol{u}}_{rel} \vert {\boldsymbol{u}}_{rel}, \end{equation}
where 
$\rho _{p}$ is the particle density, 
$C_{d}$ is the particle drag coefficient and 
${\boldsymbol {u}}_{rel}$ the particle velocity relative to the fluid, 
${\boldsymbol {u}}_{rel} = {\boldsymbol {u}}_{p} - {\boldsymbol {u}}$. The particle drag coefficient, 
$C_{d}$, is computed as a function of the particle Reynolds number, 
$Re_{p} = 2 \vert {\boldsymbol {u}}_{rel} \vert r / \nu$, according to the relation provided by Schiller & Naumann (Reference Schiller and Naumann1933),
\begin{equation} C_{d} = \begin{cases} \dfrac{4}{Re_{p}} \left(6 + Re_{p}^{2/3} \right), & \text{for}\ Re_{p} \leq 1000,\\ 0.424, & \text{for}~Re_{p} > 1000. \end{cases} \end{equation}In addition, the acceleration due to the net effect of gravity on the particle reads
\begin{equation} {\boldsymbol{f}}_{g} = \left( 1- \dfrac{\rho}{\rho_{p}} \right) \boldsymbol{g}, \end{equation}
where 
$\boldsymbol {g}$ is the gravitational acceleration.
As regards the modelling of particle–particle collisions, common methodologies include the hard-sphere approach of Campbell & Brennen (Reference Campbell and Brennen1985) and the soft-sphere approach of Cundall & Strack (Reference Cundall and Strack1979). According to the hard-sphere approach, each particle is propagated to the next collision, so that the time increment is adaptive and independent of the time step used to integrate the Navier–Stokes equations, cf. (2.6). Consequently, if the particulate phase is dense and the collision frequency between particles is high, then the time increment is reduced significantly and the computation of the particle trajectories becomes computationally expensive.
On the other hand, the soft-sphere approach allows the volumes of different particles to occupy at the same time the same space. At each instance, the resulting collisional particle acceleration is computed and deduced from the amount of overlap of their volumes. The advantage of this approach is that it allows for a constant time step. However, to obtain physically realistic results, the particle overlap has to be as limited as possible, which in turn implies that the computational time step has to be kept small.
In our computational tool pafiX, we implemented a variant of the hard-sphere approach, namely, the ray casting method. It was proposed by Roth (Reference Roth1982) in the field of computer graphics to solve a variety of intersection problems and was subsequently extended to collision detection by Schroeder (Reference Schroeder2001). For our purposes, we adapted this approach to detect collisions between spherical particles which allows the use of larger time steps. More specifically, the time step can be the same as the time step 
$\Delta t$ for the Navier–Stokes equations.
Due to the size of the time step, two particles may not be in contact with each other in the next two time instances, say at 
$t^{n+1}$ and 
$t^{n+2}$, even if they collide at some intermediate time instance (between 
$t^{n}$ and 
$t^{n+2}$). According to our implementation of the ray casting method, collisions between two particles (see figure 2) are anticipated when a number of criteria are fulfilled. These criteria are the following, sorted in order of ascending computational cost.
(i) Check whether the particles are in the same or a neighbouring computational cell. This check involves the assumption that particles propagate with a velocity of the order of the fluid phase and therefore do not transverse more than one computational cell per time step.
(ii) Check whether the vector connecting both particle centroids,
${\boldsymbol {z}}_{12}$, and the relative velocity $\boldsymbol {v}^*_1 = \boldsymbol {v}_1 - \boldsymbol {v}_2$ form an acute angle. This test examines if the particles are propagating towards each other, that is, if they are on a collision course.(iii) Determine whether the particles will collide with each other if they propagate with their current velocity, i.e. if they would collide at a later time instance.
(iv) Control if the particle collision anticipated by the previous condition takes place during the following time increment,
$\Delta t^{n+1}$.
Figure 2. Sketch of the parameters used in the ray casting collision detection approach in the rest frame of particle 2.
If any of these conditions is not fulfilled, then the subsequent conditions are not further checked and, instead, the particles are propagated with 
${\boldsymbol {f}}_{coll}=0$.
If all criteria are fulfilled, a collision is anticipated to take place. In general, this will be an oblique collision. The collision parameters are computed in the rest frame of one of the colliding particles. For the collision depicted in figure 2, the particle velocities in the rest frame of particle 2 are 
$\boldsymbol {v}^*_1$ and 
$\boldsymbol {v}^*_2$, with 
$\boldsymbol {v}^*_2 =0$. Then, particle 1 is propagated to the fictitious point 
$\boldsymbol {x}_{{c},1}$ so that it comes into contact with particle 2. The unit vector along the line that connects the centres of the two particles is 
$\boldsymbol {n} = \boldsymbol {x}_{{c},1} - \boldsymbol {x}_2$. Further, the particle velocity can be decomposed into a component in the direction of 
$\boldsymbol {n}$,
\begin{equation} \boldsymbol{v}^*_{{n},1}= \dfrac{\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{v}^*_1}{|\boldsymbol{n}|^2} \boldsymbol{n}, \end{equation}and a component that lies on the tangential plane of contact,
\begin{equation} \boldsymbol{v}_{{t},1}^* = \boldsymbol{v}^*_1 - \boldsymbol{v}^*_{{n},1}. \end{equation}
The latter component, 
$\boldsymbol {v}_{{t},1}^*$, does not change during collision.
The velocity of a particle after collision is given as the sum of the rest-frame velocity and the post-collision velocity component in the direction of 
$\boldsymbol {n}$. Accordingly, the post-collision velocity of particle 1 reads
\begin{equation} \boldsymbol{v}'_1 = \boldsymbol{v}_2 + \boldsymbol{v}_{{t},1}^* + \dfrac{r_1^3 - e r_2^3}{r_1^3 + r_2^3} \boldsymbol{v}^*_{{n},1} , \end{equation}whereas the post-collision velocity of particle 2 reads
\begin{equation} \boldsymbol{v}'_2 = \boldsymbol{v}_2 + (1+e)\dfrac{r_1^3}{r_1^3 + r_2^3} \boldsymbol{v}^*_{{n},1} , \end{equation}
where 
$r_1$ and 
$r_2$ are the radii of the particles and 
$e$ is the restitution ratio. For the derivation of (2.14) and (2.15), it was assumed that the particles are made of the same material and, therefore, the particle densities are equal. Also, for the derivation of (2.15) we took into account the fact that in the rest frame of particle 2 the velocity of this particle is identically zero.
With regard to wall–particle collisions, we remark that when a particle impacts a wall, the wall-normal velocity component changes according to
\begin{equation} \boldsymbol{v}' \boldsymbol{\cdot} {\boldsymbol{n}} = e \boldsymbol{v} \boldsymbol{\cdot} {\boldsymbol{n}} , \end{equation}
with 
$\boldsymbol {n}$ being in this case the unit vector normal to the wall. On the other hand, the tangential velocity component of the particle remains constant.
Finally, the last term in (2.8) describes the electrostatic force acting on a given particle and is given by
\begin{equation} {\boldsymbol{f}}_{{el}} = \dfrac{Q {\boldsymbol{E}}}{m_{p}}, \end{equation}
where 
$Q$ and 
$m_{p}$ are, respectively, the charge and the mass of the particle. This force is calculated by superposing the Coulomb interactions between the particle and its neighbouring particles (say, the particles that at a given time instance reside in the same computational cell) according to the hybrid method proposed by Grosshans & Papalexandris (Reference Grosshans and Papalexandris2017b) which is similar to that of Kolehmainen et al. (Reference Kolehmainen, Ozel, Boyce and Sundaresan2016). This approach is more suitable for wall-bounded flows as compared with Ewald summation and the particle–particle–particle–mesh (P
$^3$M) method (Yao & Capacelatro Reference Yao and Capacelatro2016).
Finally, it is worth mentioning that the entire algorithm is implemented in pafiX in a Fortran 90 computer code parallelized via the message passing interface protocol (MPI). As described in appendix A, the algorithm scales excellently on up to 256 processors.
3. Results and discussion
3.1. Numerical set-up
We consider gravity-driven, particle-laden flows in a duct with square cross-section, as depicted in figure 3(a). The direction of the flow is along the 
$x$-axis and aligned with the gravity vector. In other words, the gravity vector is given by 
${\boldsymbol {g}}=(9.81,0,0)\ \textrm {m}\ \textrm {s}^{-2}$. In order to mimic a very long duct, periodic boundary conditions are applied in the streamwise 
$x$-direction, whereas solid walls are assumed in the 
$y$- and 
$z$-directions. Also, the duct is assumed to be grounded and fully conductive, i.e. zero-Dirichlet boundary conditions are prescribed for the electric potential 
$\varphi _{el}$ at the walls.
Figure 3. (a) Graphic representation and dimensions of the flow domain. The arrow points in the direction of the flow. (b) Computational mesh in the cross-section, consisting of 
$60 \times 60$ cells.
The grid is refined close to the solid walls since the near-wall regions are the ones with the highest velocity gradients. The distribution of grid points in both wall-normal directions (
$y$- and 
$z$-directions) is performed in a fashion similar to that employed by Kim, Moin & Moser (Reference Kim, Moin and Moser1987) for grid refinement in the wall-normal direction in channel flows. Thus, the grid points in the 
$y$-direction are located at 
$y_j\!=\!\cos \theta _j$ with 
$\theta _j \!=\! (\,j-1) {\rm \pi}/ (N\!-\!1)$ and 
$j=1,\ldots ,N$ where 
$N$ is equal to the total number of points in the 
$y$-direction. The number of points, 
$N$, and the grid-point distribution in the 
$z$-direction are the same as in the 
$y$-direction. The resulting grid for 
$N=60$ is depicted in figure 3(b). On the other hand, the mesh is uniform in the streamwise 
$x$-direction. The DNS presented in this paper have been performed on the grid consisting of 
$180 \times 120 \times 120$ cells.
The length of the domain, 
$L$, equals six times the duct height, 
$H$, which ensures the statistical independence of the flow quantities in the streamwise direction, i.e. the flow quantities at 
$L/2$ are not correlated to the flow quantities at the periodic boundaries. Therefore, our results are independent of the size of the chosen computational domain. In the simulations presented herein, the dimensions of the duct were 
$H=4\ \textrm {cm}$ and 
$L=24\ \textrm {cm}$.
The flow parameters in our numerical simulations are presented in table 1. The friction Reynolds number is defined as 
$Re_{\tau } = u_{\tau } H / \nu$, where 
$u_{\tau } = \sqrt {\tau _{w}/\rho }$ is the friction velocity and 
$\tau _{w}$ is the stress at the wall. The Stokes number is defined as the ratio of the particle response time scale to the flow time scale, i.e. 
$St=\tau _{p}/\tau _{f}$. The particle response time scale is given by 
$\tau _{p}=2 \rho _{p} r^2 / (9 \rho \nu )$ with 
$\rho _{p}$ being the particle density. The flow time scale is based on the friction velocity, namely 
$\tau _{f} = \nu / u^2_{\tau }$. Finally, 
$Q$ and 
$\rho _N$ stand, respectively, for the electrostatic charge and number density of the particles. These quantities are expressed in a dimensionless form by the electrostatic Stokes number that we define similar to Boutsikakis et al. (Reference Boutsikakis, Fede, Pedrono and Simonin2020) as the ratio of the particle response time scale to a time scale characterizing the Coulombic particle interaction, i.e. 
$St_{el}=\tau _{p}/\tau _{el}$. Therein, 
$\tau _{el}$ was derived as 
$2 ({\rm \pi} \varepsilon m_{p}/(\rho _N Q^2))^{1/2}$. Alternatively, 
$Q$ and 
$\rho _N$ can be non-dimensionalized as the ratio of the Coulombic to gravitational forces (Kolehmainen et al. Reference Kolehmainen, Ozel, Gu, Shinbrot and Sundaresan2018), 
$F_{el}/F_{g}=Q^2\rho _N^{2/3}/(4{\rm \pi} \varepsilon m_{p}g)$. In the definition of 
$St_{el}$ and 
$F_{el}/F_{g}$ the reference length scale is the distance between the centres of two particles assuming a uniform distribution in the domain.
Table 1. Parameters of the numerical simulations.
According to table 1, we investigate the effect of electrical forces for three different flow conditions: First, a flow of 
$Re_{\tau }=600$ and 
$\rho _{p} / \rho =1000$ is simulated. Further, 
$Re_{\tau }$ is reduced to 300 and, finally, 
$\rho _{p} / \rho$ is increased to 7500. For these three cases, the flow and particle dynamics is studied for uncharged particles. Then, we investigate how these flows are affected if each particle carries an electrostatic charge of 
$Q=0.1\ \textrm {pC}$. More specifically, the particle number density in all cases is 
$10^8\ \textrm {m}^{-3}$. This implies that the particulate phase inside the computational domain consists of 38 400 particles. Also, the particles are monodisperse and spherical, with a radius 
$r=25\ \mathrm {\mu }\textrm {m}$. It is worth mentioning that we assume there is no charge exchange during particle–particle or wall–particle collisions.
The validation of the flow solver is presented in appendix B. In addition to the parameters in table 1, we explored the effect of higher charges on the particles. These results are briefly discussed in appendix C.
3.2. Particle concentrations and dynamics
In this subsection, we present results for the particle concentration and particle dynamics for the six cases mentioned in table 1. For each case, we started the simulation without particles, until the flow turbulence becomes fully developed. Subsequently, we seeded particles at random locations and with velocities equal to those of the fluid in the same locations. In order to save computing time and since particles tend to migrate from high to low turbulence-intensity regions due to turbophoretic drift, the injection was such that the initial particle concentrations were higher close to the walls of the duct. It is noted, however, that the predicted properties of statistically stationary particle-laden flows are totally independent of the initial location of particles. The charge was assigned to the particles instantly when they were embedded in the flow. Accordingly, electrostatic forces between particles modify the dynamics of the particulate phase from the beginning.
First, we examined the temporal evolution of the particle concentrations in the near-wall regions. The results for all six cases are presented in figure 4. In these plots, the time is non-dimensionalized in terms of the frictional velocity, i.e. 
$t^+=t u^2_{\tau }/\nu$. Also, the distance from a wall is given in terms of wall units, i.e. 
$y^+$ and 
$z^+$. For each flow case, the three curves shown in the figure represent the normalized number density 
$C$ of the particles residing within three different distances from the walls of the duct. The number density 
$C$ is normalized with the overall number density of the particles in the duct 
$\rho _N$. For example, for flows at 
$Re_{\tau }= 600$, the black solid lines depict the normalized number density 
$C$ of the particles confined within a very thin layer at the walls. The thickness of this layer is equal to one viscous length scale. Also, the blue dashed curves show the temporal evolution of the particles located within the viscous sublayer, i.e. less than five viscous length scales from the wall; inside this layer viscous stresses dominate over Reynolds stresses. Further, the green dotted lines show the normalized number density of the particles located in the viscous wall region which extends to 50 viscous length scales, where both Reynolds and viscous stresses are significant.
Figure 4. Temporal evolution of the normalized particle number density 
$C$ near the walls: (a) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho = 1000$, 
$F_{el}/F_{g}= 0$; (b) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (c) 
$Re_{\tau }=300$, 
$\rho _{p}/\rho = 1000$, 
$F_{el}/F_{g}= 0$; (d) 
$Re_{\tau }=300$, 
$\rho _{p}/\rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (e) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho = 7500$, 
$F_{el}/F_{g}= 0$; (f) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho = 7500$, 
$F_{el}/F_{g}= 0.004$; (black solid line) 
$y^+$ or 
$z^+$ 
$<1$ for 
$Re_{\tau }= 600$ and 
$y^+$ or 
$z^+$ 
$<0.5$ for 
$Re_{\tau }= 300$: (blue dashed line) 
$y^+$ or 
$z^+$ 
$<5$ for 
$Re_{\tau }= 600$ and 
$y^+$ or 
$z^+$ 
$<2.5$ for 
$Re_{\tau }= 300$; (green dotted line) 
$y^+$ or 
$z^+$ 
$<50$ for 
$Re_{\tau }= 600$ and 
$y^+$ or 
$y^+$ 
$<50$ for 
$Re_{\tau }= 300$.
In figure 4 we readily observe that, in all cases, the concentrations in the viscous wall region reach much faster their final values, i.e. their values when the flow of the mixture has become statistically stationary. The particles are transported from the bulk of the flow towards this region by the large eddies of the flow. Since these eddies have the highest amounts of kinetic energy, the transportation of the particles into this region occurs relatively fast. Once inside the viscous wall region, the particles adopt their preferential locations but are driven by small-scale fluctuations; this is, therefore, a slower process. For the interpretation of the plots shown in figure 4, it is worth recalling that the frictional quantities of the flow of 
$Re_{\tau }= 300$ are different from those of 
$Re_{\tau }= 600$; this has an impact on the normalization of all spatial coordinates and time.
In some cases there is an overshoot of 
$C$ in the early stages of the evolution of the particulate phase, as can be seen in figures 4(a), 4(b) and 4(f). More specifically, the particles are first accelerated from their initial locations towards a location very close to the wall. Subsequently, they move away from the wall and settle at their preferential location at a distance further from it. As mentioned above, the initial particle positions are random and do not necessarily correspond to the actual dynamics of fully developed turbulence. Therefore, this overshoot characterizes the transition of the particulate phase from its initial state to the physically meaningful and statistically stationary state.
Further, the curves in the right column of figure 4 show that the effect of electrostatic forces takes place gradually and not abruptly. Upon comparison with the left column it can be deduced that, when the particles are charged, they assume their stationary-state positions faster than when they do not carry any charge. This is particularly noticeable in the comparison of figure 4(a) with 4(b) as well as in the comparison of figure 4(c) with 4(d). Further, this observation is in accordance with the finding of Grosshans (Reference Grosshans2018) that charges of the same polarity dampen the particle velocity fluctuations.
Also, from figure 4, it is readily inferred that the particle concentrations at the statistically stationary state vary significantly from one case to the other, depending on the underlying turbulence dynamics and the charge they carry. In the following paragraphs, we explore this phenomenon and examine the properties of the particle concentrations when the flow has become statistically stationary. Statistics have been collected for a period of 
$t^+ = 17\,500$. This time period is equivalent to 100 flow-through times, i.e. the ratio of the length of the computational domain to the bulk velocity of the fluid.
Instantaneous snapshots of the particle positions in the cross-section of the duct (
$y$–
$z$ plane) are presented in figure 5 for each case examined herein. These snapshots were taken after the flow became statistically stationary. Also, for visualization purposes, only one in every five particles is shown. These images are complemented with the snapshots of the mean particle concentration 
$\langle C \rangle$ in figure 6. Therein, 
$\langle C \rangle$ is normalized by the maximum local concentration for each case. The bin size for establishing these histograms is 
$30\times 30$ viscous length scales for flows at 
$Re_{\tau }=600$, and 
$15\times 15$ for the flows at 
$Re_{\tau }=300$.
Figure 5. Instantaneous particle positions in the cross-section of the duct (
$y$–
$z$ plane), recorded once the flow has become statistically stationary: (a) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}= 0$; (b) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}= 0.026$; (c) 
$Re_{\tau }=300$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}=0$; (d) 
$Re_{\tau }=300$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}= 0.026$; (e) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =7500$, 
$F_{el}/F_{g}=0$; (f) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =7500$, 
$F_{el}/F_{g}=0.004$. The colours indicate the streamwise velocity of the particles; for each case, the red colour corresponds to the fastest and the blue colour to the slowest particle. For visualization purposes the particle size is scaled up and only one in every five particles is shown.
Figure 6. Mean particle concentration in arbitrary units, 
$\langle C \rangle$, in the cross-section of the duct: (a) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}=0$; (b) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}= 0.026$; (c) 
$Re_{\tau }=300$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}=0$; (d) 
$Re_{\tau }=300$, 
$\rho _{p}/\rho =1000$, 
$F_{el}/F_{g}= 0.026$; (e) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =7500$, 
$F_{el}/F_{g}=0$; (f) 
$Re_{\tau }=600$, 
$\rho _{p}/\rho =7500$, 
$F_{el}/F_{g}= 0.004$. Note that 
$\langle C \rangle$ is normalized by the maximum local concentration of each case, i.e. the vertical axes of the figures are scaled differently.
As a first observation, these figures confirm that particles tend to accumulate close to the walls, as mentioned in the Introduction. Nonetheless, despite this common trend, the resulting particle distributions vary both quantitatively and qualitatively from one case to another. According to the figures, concentration peaks appear at two different locations, namely in the corners and in the bisectors (the symmetry lines of each wall) of the duct. For the case of 
$Re_{\tau }=600$ and 
$\rho _{p}/\rho =1000$, particles accumulate preferably at the bisectors whereas relatively few particles reside near the corners of the duct, as can be seen in figures 5(a) and 6(a). Furthermore, when the particles are charged, figures 5(b) and 6(b), the particle density close to the walls increases significantly, especially the peaks at the bisectors.
In the cases of low frictional Reynolds number, 
$Re_{\tau }=300$, the particle concentration peaks at both locations, i.e. corners and wall bisectors. In particular, the snapshot in figure 5(c) clearly illustrates the impact of secondary flows of the second kind on the particle dynamics. Therein, one can identify a high-concentration region at each wall that is stretching towards the centreline of the duct. These noticeable structures correspond to the arrow C sketched in figure 1 and are induced by the vortical motion of the gas. Thus, at these regions, the particles are ejected from the wall back towards the bulk flow. As can be inferred from figure 5(c), the location of these ejection points is not exactly at the bisectors. Instead, our simulations predicted that these ejection points move with a low temporal frequency back and forth along the wall. When the particles carry electrostatic charge, they are pushed against the wall and no longer follow the turbulent eddies to the bulk of the flow. Consequently, as can be seen in figure 5(d), only a few particles are airborne.
Finally, it is worth noting that, for the case with a high density ratio, figure 6(e), the mean concentration 
$\langle C \rangle$ has peaks only in the corners and not at the walls. Furthermore, the electrostatic charge carried by the particles, figures 5(f) and 6(f), does not appear to alter this general flow pattern.
Further, with regard to the mean particle concentration, comparison of figures 6(a), 6(c) and 6(e) shows that the flow parameters employed in our study result in qualitatively different concentration profiles. Further, the effect of electrostatic charge is distinctly different for each case and depends on the underlying dynamics of the fluid flow. It is conjectured that this difference is linked to the strength of the vortical structures of the carrier gas (i.e. the turbulence intensity) and illustrates the sensitivity of the particles to these structures.
In order to examine in more detail the influence of the electrostatic charge, in figure 7 we present plots of the particle concentrations in four different stations. As indicated in the sketch of figure 1, a duct with square cross-section possesses eight symmetry planes. Owing to this symmetry, we have plotted results for only one quadrant of the duct in figure 7. The horizontal axes on these plots are scaled differently for 
$Re_{\tau }= 300$ and for 
$Re_{\tau }= 600$ and in such a way that the physical coordinates of all cases are consistent with each other. In this figure, the solid lines correspond to flows with uncharged particles, whereas the dashed lines correspond to the equivalent flows with charged particles.
Figure 7. Mean normalized particle number density, 
$\langle C\rangle$, in four different slices: (a) 
$2.5<z^+<5$ for 
$Re_{\tau }= 300$ and 
$5<z^+<10$ for 
$Re_{\tau }= 600$; (b) 
$10<z^+<12.5$ for 
$Re_{\tau }= 300$ and 
$20<z^+<25$ for 
$Re_{\tau }= 600$; (c) 
$27.5<z^+<30$ 
$Re_{\tau }= 300$ and 
$55<z^+<60$ for 
$Re_{\tau }= 600$; (d) 
$147.5<z^+<150$ for 
$Re_{\tau }= 300$ and 
$295<z^+<300$ for 
$Re_{\tau }= 600$. For visibility purposes, in (a) the horizontal axis is linear, whereas in (b)–(d) it is logarithmic. In each figure, the top horizontal axis corresponds to the flows at 
$Re_{\tau }= 300$ and the lower one to flows 
$Re_{\tau }= 600$; (black solid line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}=0$; (black dashed line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (red solid line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
${F_{el}/F_{g}=0}$; (red dashed line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (blue solid line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}=0$; (blue dashed line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}= 0.004$.
Figure 7(a) shows the profile of the mean normalized number density 
$\langle C\rangle$ in a slice very close to the wall of the duct. As noted earlier, for the case with 
$Re_{\tau }= 600$ and 
$\rho _{p}/\rho = 1000$, 
$\langle C\rangle$ does not exhibit a peak in the corners of the duct. Instead, according to figure 7(a), 
$\langle C\rangle$ peaks at the bisectors of the wall. Moreover, the peak value is significantly increased, by a factor of 2.6, when the particles are electrostatically charged. By contrast, the particle concentration in the corners is not significantly affected by the electrostatic charge.
The case of 
$Re_{\tau }= 300$ and 
$\rho _{p}/\rho = 1000$ is characterized by higher concentration peaks in the corners and the walls compared with the other cases. Also, for this specific case, the particle concentration across the other slices, shown in figures 7(b)–7(d), exhibits the highest gradients among all cases examined herein. Thus, in this case, the particle concentration is most influenced by the turbulence dynamics. This implies that the particles are most sensitive to variations of the attacking forces. Accordingly, the presence of electrostatic forces can dramatically influence the particle concentrations. More precisely, the repelling forces push the particles towards the walls of the duct where they accumulate.
By comparison, in the flow at 
$Re_{\tau }= 600$ and 
$\rho _{p}/\rho = 7500$, the particle concentration is only slightly affected, as can be seen in figure 7(a). The importance of inertial forces is also indicated by the relatively high 
$F_{el}/F_{g}$ of 0.026. The electrostatic charge results in an increase in the local concentration of approximately 30 % along the wall. The same can be inferred from figures 7(b) and 7(c) where the concentration in the bulk of the flow is not affected by the electrostatic charges. In fact, in this region the normalized particle number density is close to unity. However, the peak close to the wall, 
$y^+<30$, increases by as much as 30 %. The effect of electrostatic charge in the case 
$Re_{\tau }= 600$ and 
$\rho _{p}/\rho = 1000$ and in the region of the slices examined in figures 7(b) and 7(c) is similar. When the particles are charged, the particle concentration increases close to the walls but is only slightly modified in the bulk of the flow.
The most drastic effect of the electrostatic forces can be observed in the slice which characterizes the centre plane of the duct, as shown in figure 7(d). The highest concentration of particles occurs in the case with 
$Re_{\tau }= 300$ and 
$\rho _{p}/\rho = 1000$ and uncharged particles. The particle concentration decreases with the distance from the wall until it reaches a minimum at 
$y^+\approx 75$. In the corresponding case with electrostatic charging, the particles migrate towards the wall (due to the scaling of the axis, this case is not visible in figures 7d). For the flow at 
$Re_{\tau }= 600$ and 
$\rho _{p}/\rho = 7500$ the effect of the electrostatic charges is still significant albeit less pronounced. The particle concentration at the walls is increased by 54 % and the peak in the bulk of the duct is reduced by 27 %. The influence of electrostatic charges in the case with 
$Re_{\tau }= 600$ and 
$\rho _{p}/\rho = 1000$ is also very significant. More specifically, without electrostatic charge the concentration peak is located at 
$y^+=18$, whereas when the particles are charged, the peak is located right next to the wall. Additionally, this peak is increased by a factor of five. Further, the concentration at the centre is reduced by 31 % and the corresponding peak nearly completely diminishes.
In figure 8 we present plots of the mean particle velocities. These plots clearly substantiate the effect of secondary flows on the particle dynamics. As before, by virtue of the symmetries of the flow domain, we present results for one wall-normal velocity component and in one quadrant of the duct. In general, the flows at low frictional Reynolds number, 
$Re_{\tau }= 300$, exhibit the highest wall-normal velocities (in absolute terms). On the other hand, the flows at high density ratio, 
$\rho _{p} / \rho = 7500$, exhibit the lowest ones. These trends are directly related to the difference in the particles’ Stokes number which measures the sensitivity of the particles to the velocity fluctuations of the surrounding fluid.
Figure 8. Normalized mean particle velocity in a wall-normal direction in four different slices: (a) 
$2.5<z^+<5$ for 
$Re_{\tau }= 300$ and 
$5<z^+<10$ for 
$Re_{\tau }= 600$; (b) 
$22.5<z^+<25$ for 
$Re_{\tau }= 300$ and 
$45<z^+<50$ for 
$Re_{\tau }= 600$; (c) 
$35<z^+<37.5$ for 
$Re_{\tau }= 300$ and 
$70<z^+<75$ for 
$Re_{\tau }= 600$; (d) 
$147.5<z^+<150$ for 
$Re_{\tau }= 300$ and 
$295<z^+<300$ for 
$Re_{\tau }= 600$. In each panel, the top horizontal axis corresponds to the flows at 
$Re_{\tau }= 300$ and the lower one to the flows at 
$Re_{\tau }= 600$: (black solid line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}=0$; (black dashed line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (red solid line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
${F_{el}/F_{g}=0}$; (red dashed line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (blue solid line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}=0$; (blue dashed line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}= 0.004$.
The positive velocities in figure 8(a) at distances larger than 
$y^+\approx 50$ are due to particles moving in the direction of the arrow B of figure 1. In other words, the particles move away from the corner of the duct and towards the bisectors of the walls. As can be inferred from the plots shown in figure 8(a), this particle motion is substantially subdued when the particles carry electrostatic charge. As a matter of fact, for the case at 
$Re_{\tau }= 300$ and 
$\rho _{p} / \rho = 1000$, the dominance of the electrostatic forces over the aerodynamic ones completely inhibits the vortical motion of particles. For this reason, the velocity profiles for this case have not been included in figure 8(a).
The particle velocities in the region 
$y^+<50$ relate to the tail of the structure of particles that are transported from the duct's centreline to the corners in the direction of the arrow A of figure 1. A more detailed perspective of this outward flux is provided in figures 8(b) and 8(c); therein, this outward flux is represented by the peaks of the wall-normal velocity amplitudes at 
$y^+\approx 60$ and 
$y^+\approx 90$, respectively. Interestingly, this part of the vortical motion of the particles is not affected by the electrostatic field; instead, the motion of charged particles is very similar to that of uncharged ones.
The inward particle flux from the wall bisectors towards the centreline of the duct follows the direction of the arrow C in figure 1 and can be observed in figure 8(d). Similarly to the flux along the wall, the particle velocity in this direction is strongly reduced in the case of charged particles. For both types of motions, the extent of the influence of electrostatic forces on the particle velocity is stronger for the cases with 
$\rho _{p} / \rho = 1000$ than for the cases with 
$\rho _{p} / \rho = 7500$.
In summary, both the particle flux along the walls and the inward motion from the wall towards the bulk of the duct are significantly reduced when the particles carry electrostatic charge. On the other hand, the particle transport from the duct's centreline towards the corners is not affected by the electrostatic forces. This finding explains the particle concentration profiles shown in figure 7. In particular, due to electrostatic charges, the particles are driven with roughly the same velocity towards the walls but are slowly ejected back to the centre. As a result, the particle concentration increases at the walls and especially in the corners. Consequently, the notable difference in the particles’ wall-normal velocity between charged and uncharged particles when 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$ substantiates the strong increase of the particle concentration in the vicinity of walls.
To illustrate the direction of the arising electrostatic forces, instantaneous snapshots of the electric potential 
$\varphi _{el}$ are presented in figure 9. We recall that the electric field strength 
${\boldsymbol {E}}$ is equal to 
$-\boldsymbol {\nabla } \varphi _{el}$. Accordingly, one can deduce from these plots that the particle encounters a higher electrostatic force when moving from the bisectors of the walls towards the centreline of the duct than when moving along the diagonal from the centreline to the corner. This implies that, the emerging electric field significantly reduces particle circulation due to secondary flows of a second kind.
Figure 9. Instantaneous snapshots of the electric potential across the cross-section of the duct, induced by charges carried by the particles: (a) 
$Re_{\tau }=600$, 
$\rho _{p} / \rho =1000$; (b) 
$Re_{\tau }=300$, 
$\rho _{p} / \rho =1000$; (c) 
$Re_{\tau }=600$, 
$\rho _{p} / \rho =7500$.
Next, we investigate the influence of electrostatic charges on the streamwise velocity of the particles. Plots of the particle mean streamwise velocity component are provided in figure 10. Since particles follow to a certain extent the gaseous phase, the particle streamwise velocities close to the walls are lower than those in the bulk of the duct. Moreover, as mentioned previously, in ducts with square cross-section, secondary flows play an important role in the transport of particles in the wall-normal directions. In turn, the streamwise momentum of particles is also transported with them along the 
$y$- or 
$z$-axis. More specifically, the vortical flow structures transport fast-moving particles towards the near-wall region and slow-moving ones back to the bulk of the duct. Once the slow-moving particles enter the bulk of the duct, they begin to accelerate. This momentum transport contributes to a more uniform streamwise velocity distribution as can be inferred from figure 10(d). In the slice shown in this figure, which corresponds to particle displacement along the arrow C of figure 1 the particles in the near-wall region (
$y^+<5$ or 
$y^+<20$, depending on the case), are faster than their surrounding fluid owing to the momentum they received when they were in the bulk of the duct. When moving towards the centreline, beyond 
$y^+<5$, respectively 
$y^+<20$, the surrounding gas is faster than the particles which, therefore, experience an acceleration in the 
$x$-direction. Then, once they reach the centreline, their streamwise velocity is maximized. The fact that the maximum particle velocity differs from one case to another and is not equal to the maximum gas velocity is due to the interplay between gravity and the aerodynamic forces acting on the particles.
Figure 10. Normalized mean particle velocity in the streamwise direction in four different slices as a function of the distance to the duct's wall: (a) 
$2.5<z^+<5$ for 
$Re_{\tau }= 300$ and 
$5<z^+<10$ for 
$Re_{\tau }= 600$; (b) 
$22.5<z^+<25$ for 
$Re_{\tau }= 300$ and 
$45<z^+<50$ for 
$Re_{\tau }= 600$; (c) 
$35<z^+<37.5$ for 
$Re_{\tau }= 300$ and 
$70<z^+<75$ for 
$Re_{\tau }= 600$; (d) 
$147.5<z^+<150$ for 
$Re_{\tau }= 300$ and 
$295<z^+<300$ for 
$Re_{\tau }= 600$. For visibility purposes, in (a) the horizontal axis is linear, whereas in (b–d) it is logarithmic. The upper axis in each figure relates to the cases of 
$Re_{\tau }= 300$ and the lower axis to the cases of 
$Re_{\tau }= 600$. In addition, in (d) the mean velocity of the undisturbed gaseous phase is displayed. (black solid line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}=0$; (black dashed line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (red solid line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}=0$; (red dashed line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}= 0.026$; (blue solid line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}=0$; (blue dashed line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}= 0.004$.
Finally, we examine the effect of the electrostatic field on the particles moving close to the wall and in the direction of the arrow B of figure 10(a). Comparison of the computed streamwise velocities corroborates the fact that when the particles carry no charge they move significantly faster. Indeed, as can be observed in figure 8(a), the charged particles are slowed down in the wall-normal direction by the electrostatic forces and remain close to the walls for a long time. Due to longer residence times, more streamwise momentum is exchanged between particles and gas, which results in the slowing down of the particles.
In contrast, the electrostatic forces do not influence the velocity of particles moving outward diagonally, in the direction of the arrow A. With regard to this type of motion, charged and uncharged particles move at the same wall-normal velocity, which in turn implies that they exchange a similar amount of streamwise momentum with their surrounding gas. As a result, their streamwise velocities are also similar, as can be inferred from figures 10(b) and 10(c) beyond 
$y^+\approx 30$. These figures also confirm that the particles in the region 
$y^+ < 30$, i.e. those moving in the direction of the arrow B, are transported slower in the 
$x$ direction.
On the basis of these observations, we conclude that electric charges not only attenuate particle motion induced by secondary flows but they also substantially reduce the streamwise momentum transfer between particles and carrier gas over the cross-section of the duct.
4. Conclusions
We studied, by means of DNS, the dynamics of particle-laden flows through a duct of a square-shaped cross-section. According to our simulations, when the particles are charged, the electrostatic forces dampen significantly the particles’ vortical motions that are induced by the secondary flows of the carrier gas. The charged particles still migrate in the diagonal direction from the centreline of the duct towards its corners. But on their way back from the wall to the centre their velocity is reduced. This modification to the particle dynamics results in significantly different characteristics of the particle number density. The reduction of the vortical motion of particles in the regions of secondary flows leads to an increase of the particle concentration at the walls, especially at the bisectors of the walls and the corners of the duct. Also, the streamwise momentum transfer in the cross-section, which relies on the wall-normal motion of particles, is significantly attenuated. These results demonstrate the fundamental influence of electric forces on the emerging pattern of dispersed two-phase flows. Such an understanding of the underlying mechanisms opens new perspectives for the control of the flow dynamics of powders. It would be interesting to test this influence also for rectangular ducts of different aspect ratios and to elaborate in detail on the dependence of the vortical particle motion on a successive increase of the powder charge.
Acknowledgements
H.G. and C.B. gratefully acknowledge the financial support from the DECHEMA Max Buchner Research Foundation and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 947606 PowFEct). M.C. and M.V.P. gratefully acknowledge the financial support of the National Research Fund of Belgium (FNRS) under the FLOW-CHARGE grant.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Scalability
As a first test, we examined the scalability of pafiX. More specifically, we simulated flows in a duct with square cross-section for which the computational domain, depicted in figure 3(a), was discretized by 
$256\times 120\times 120$ grid cells. In our tests, two different cases were considered. The first one is single-phase flow, i.e. without solid particles, at 
$Re_{\tau }= 600$. In this manner, only the scalability of the Navier–Stokes solved is accessed. The second one is flow laden with 
$10^8\ \textrm {m}^{-3}$ charged particles of 
$Q=0.1\ \textrm {pC}$; this is also the second case given in table 1.
With regard to domain decomposition, the domain was decomposed in the 
$x$-direction. For the geometry considered herein, this is the most efficient domain decomposition in terms of load balance between processors. Evidently, this choice also poses an upper limit to the number of processors that can be used. For the grid employed herein this limit is 256 processors. The speed-up factor of the simulation, 
$S_k$, is defined as the ratio of the run time of the parallel code running on a single processor, 
$T_1$, to the run time of the same code running on 
$k$ processors,
\begin{equation} S_k = \dfrac{T_1}{T_k}. \end{equation}
The resulting speed-up for both cases and for 
$k=2^n$ (
$n=0,1, \ldots , 8$) is plotted in figure 11. According to figure 11, the code scales very well in the range of number of processors considered herein. Our tests confirmed that, as expected, the fraction of the MPI communication time, i.e. the time needed to send and receive data packages from one processor to the other, compared with the total computing time, increases with the degree of domain decomposition. More specifically, for the case of single-phase flow, it increases from 1.3 % for 
$k=2$ to 54.5 % for 
$k=256$. The fact that the speed-up is not significantly reduced can be attributed to the non-optimized memory access procedure of the code, which results in a rather expensive calculation when 
$k$ is small and the variable matrices are large.
Figure 11. Scalability of the fluid solver and the fluid plus the particle and electrostatic solver. The straight dashed lone denotes the ideal linear scaling.
Moreover, a particularly long computing time was recorded for the simulation of the second test case (particle-laden flow with charged particles) on only one processor. Upon detailed inspection, it was revealed that this is directly related to the criterion (i) of the collision algorithm described in § 2. In fact, the computational cost of this criterion scales with 
${O}(N_p^2)$ where 
$N_p$ is the number of points per subdomain and scales with 
$1/k$. Accordingly, the computational expense of the solver for the particulate phase is significantly reduced if 
$k>1$.
It is also worth mentioning that, for each case, the results of these simulations were identical, i.e. independent of 
$k$. In summary, the proposed algorithm scales very well for 
$k>1$. The numerical simulations reported in the following were performed each on 32 processors.
Appendix B. Code validation
For validation purposes of the Navier–Stokes solver, we performed DNS of duct flows without particles and compared them with earlier DNS that appeared in the literature. For flow at 
$Re_{\tau }= 600$, we used five different meshes, ranging from 
$60\times 60\times 60$ to 
$200\times 140\times 140$ cells in the 
$x$-, 
$y$- and 
$z$-directions, respectively. In figure 12(a) we show plots of the mean streamwise velocity profile at the centre plane of the duct, 
$y=H/2$, as functions of the distance from the wall measured in wall units, 
$y^+$. For comparison purposes, the DNS data of Huser & Biringen (Reference Huser and Biringen1993) are also plotted in this figure. It can be readily inferred that the velocity profiles converge as the grid resolution is increased. In particular, the results obtained on the two finer grids, 
$180\times 120\times 120$ and 
$200\times 140\times 140$ cells respectively, are identical and match the DNS data of Huser & Biringen (Reference Huser and Biringen1993).
Figure 12. DNS of flow in a duct without particles: (a) 
$Re_{\tau } = 600$; (b) 
$Re_{\tau } = 300$. Comparison of the mean streamwise velocity at the centre plane of the duct with the DNS data of Huser & Biringen (Reference Huser and Biringen1993) and Sharma & Phares (Reference Sharma and Phares2006). The profiles for 
$Re_{\tau }= 600$ are computed on five different grid resolutions, as shown in the legend of the figure.
Moreover, we performed a simulation of a flow at 
$Re_{\tau }= 300$ on a grid consisting of 
$180\times 120\times 120$ cells. Our numerical predictions for the mean streamwise velocity at the centre plane are shown in figure 12(b). As can be seen, our predictions match those of Sharma & Phares (Reference Sharma and Phares2006).
Further, the root-mean-square (r.m.s.) values of the velocity components are presented in figure 13. From the results for the flow at 
$Re_{\tau }= 600$ (figures 13a–13c) we can infer once again that numerical grid convergence is attained with grid refinement. For the fine grid resolutions, only slight discrepancies are observed between our predictions for the r.m.s. of the wall-normal velocity components, 
$v_{rms}$ and 
$w_{rms}$ and those of Huser & Biringen (Reference Huser and Biringen1993) and Zhu et al. (Reference Zhu, Yang and Chen2009).
Figure 13. DNS of flow in a duct without particles: (a) 
$Re_{\tau } = 600$; (b) 
$Re_{\tau } = 600$; (c) 
$Re_{\tau } = 600$; (d) 
$Re_{\tau } = 300$; (e) 
$Re_{\tau } = 300$; (f) 
$Re_{\tau } = 300$. Comparison of the r.m.s. velocity components with the DNS data of Huser & Biringen (Reference Huser and Biringen1993), Zhu, Yang & Chen (Reference Zhu, Yang and Chen2009) and Sharma & Phares (Reference Sharma and Phares2006). The profiles for 
$Re_{\tau }= 600$ are computed on five different grid resolutions, as shown in the legend of the figure. Note that Sharma & Phares (Reference Sharma and Phares2006) did not provide data for 
$v^+_{rms}$.
The r.m.s. of the velocity components for the flow at 
$Re_{\tau }= 300$ are plotted in figures 13(d)–13(f). In general, our results compare favourably with the DNS data Sharma & Phares (Reference Sharma and Phares2006) (note that Sharma & Phares (Reference Sharma and Phares2006) did not provide the fluctuations of the velocity component 
$v$). The most noticeable difference is in the 
$u_{rms}$ profile far from the wall for which the predictions of Sharma & Phares (Reference Sharma and Phares2006) are 10 %–20 % lower than ours.
On the basis of these tests and the favorable comparison with earlier results, the DNS presented below have been performed on the grid consisting of 
$180\times 120\times 120$ cells.
Appendix C. Higher charge levels
Exemplary results of an exploratory study considering a higher charge level of 0.2 pC on each particle are presented in figure 14. In this figure, we replotted figures 7(b) and 7(d) and added the mean particle concentration profiles of the additional cases. For this new charge level qualitatively the same phenomena occur as for the lower particle charge. The vortical motion of the particles in the cross-section is even more retarded by the arising electrostatic forces. Consequently, even more particles accumulate close to the walls.
Figure 14. Mean normalized particle number density, 
$\langle C\rangle$, replotted from figures 7(b) and 7(d) with additional cases considering a higher charge level: (a) 
$10<z^+<12.5$ for 
$Re_{\tau }= 300$ and 
$20<z^+<25$ for 
$Re_{\tau }= 600$; (b) 
$147.5<z^+<150$ for 
$Re_{\tau }= 300$ and 
$295<z^+<300$ for 
$Re_{\tau }= 600$; (black dotted line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}=0.104$; (red dotted line) 
$Re_{\tau }= 300$, 
$\rho _{p} / \rho = 1000$, 
$F_{el}/F_{g}=0.104$; (blue dotted line) 
$Re_{\tau }= 600$, 
$\rho _{p} / \rho = 7500$, 
$F_{el}/F_{g}=0.017$.
For clarity of the presentation, we focus on the main results section of this paper on one charge level per flow condition.
