A new computational algorithm for the interaction between electrically charged particles

In industrial applications, the dynamics of particles is frequently controlled through electrical charges, e.g., in electrostatic precipitators or during powder coating. However, the electrification of particulates can also cause the formation of deposits, hazardous sparks, and dust explosions. The objective of our work is to propose a new computational method which captures accurately the dynamics of interacting electrically charged particles during their approximation. This model focuses especially on the precise prediction of the contribution of the electrostatic and collisional forces whose time-scales are typically much smaller than the numerical time-step used to compute forces other than the Coulomb force. To this end, binary particle interaction is calculated through the local adaptive refinement of the time-step. We present results that demonstrate the capability of the method to accurately and efficiently describe binary and multiple particle interaction. Further, through comparison with benchmark solutions we elaborate on the conditions, in terms of particle charges and sizes, for which our model is superior to previously employed approaches.


Introduction
The dynamic interaction of electrically charged particles is an elementary mechanism of many industrial processes. Examples include the control of particle trajectories during electrostatic precipitation [1,7], powder coating [30], pneumatic conveying [29], or printing [25]. As regards granular media, recent studies demonstrated the huge impact of an excess of electric charges on the interactions between grains [20,23]. Not only solid material, also droplets can acquire a net excess of electric charges [3]. The influence of these electric charges on the collision between droplets is also primordial as it affects natural phenomena such as thunderclouds [4]. Nonetheless, electrostatics may also alter particle dynamics in an unwanted way such as the attraction of charges on the particle and its image on the surface of transport pipes leading to the formation of deposits [17,18].
The fundamental physics and mathematical equations underlying this interaction are well-established. Namely, Newton's second law of motion describes inertial acceleration, Gauss' law (also known as Maxwell's first equation) gives the electric field forces, and the Navier-Stokes equations in case aerodynamic interaction is relevant. Albeit being generally accepted, the numerical solution of a coupled system of these equations poses challenges. Let alone the computation of collisions between uncharged particles, two approaches became popular, namely the hardsphere and the soft-sphere approach. In the hard-sphere approach [5], the particle positions are advanced within one numerical time-step until the time instance when next inter-particle collision is expected. The outcome of the collision can be calculated using the laws of classical physics. However, this procedure becomes numerically expensive when considering one of the above mentioned real-scale systems consisting of a large number of particles which collide frequently. In the soft-sphere approach [6] the timestep is kept constant which makes it more suitable for large systems. At each time-step the collisional forces are computed based on the overlapping between particles, e.g., through the elastic theory of Hertz. Nevertheless, only a small overlapping is allowed to obtain physically realistic results. Thus, the restriction in the maximum particle deformation directly leads to a rather low upper limit of the computational time-step.
Although the hard-and soft-sphere approaches function satisfactory for certain flow systems, their focus lies on the correct prediction of pure collisional dynamics. However, if charges are involved, particles interact also strongly while they are not in physical contact. In fact, electrostatic forces increase exponentially with decreasing distance between charges. The effect of this increase is visualized in Fig. 1 which gives the schematics of the interaction between two inertial particles under the influence of an electrostatic field. Therein the exact, i.e., correct trajectory of particle 1 is highly deflected by the repelling inter-particle forces. On the other hand, the numerically computed particle trajectory is qualitatively and quantitatively different since the intermediate change of electrostatic forces during the finite time-step is not accounted for. Therefore, the time-step would be required to be significantly lowered so that not only each collision event but the complete trajectories during particle approximation and departure are well resolved. Since this would result in high computational expenses, most often (e.g., [19,21,28]) inter-particle electrostatic effects are neglected during their close interaction which may be reasonable for low-charged systems. However, in the case of moderately and highly charged particles both these concepts, hard-sphere and soft-sphere approach in their standard formulations, fail.
The authors of the current paper develop the numerical tool [24] which aims to compute the electrostatic charging of powder flows. Previously, we used this tool to conduct a wide range of studies regarding charge build up in circular pipe flows [9][10][11][12], channel flows [14], cubical boxes [15], and even atmospheric flows [13,16] (for an overview, please refer to Ref. [8]). Therein, we applied a simple model to predict the charge exchange between a particle and a surface upon contact. In the present paper we investigate the behavior of particles once they have acquired a charge. More specifically, we propose a new numerical method to compute the interaction between particles under the influence of electrostatic forces. The following section details the new computational method. In Sect. 3.1 the method is tested for the case of only two particles. Afterward, in Sect. 3.2 multiple particles are considered. Finally, in the last section the study is concluded.

Methods
In the present paper, improvements of numerical schemes for the treatment of charged particles are proposed. The first, which we call here standard model since we applied it previously (e.g., [14]) to compute the outcome of inter-particle collisions, is described in Sect. 2.1. This method is further extended through the precise inclusion of electrostatic forces, which is detailed in Sect. 2.2. Both these approaches are introduced separately because their numerical predictions are separately evaluated and compared in Sects. 3.1 and 3.2. It is noted that for the sake of completeness the models are described including aerodynamic forces whereas for simplicity the tests are performed assuming vacuum to fill the free space.

Structure of the standard approach
In free space the nth particle position, x j+1 n , and velocity, u j+1 n , are evolved from time-step j to j + 1 by t using an Eulerforward scheme. The physical time is defined as i.e., t j is the time instant at time-step j. Thus, the positions and velocities of each particle (while dropping the index n) are propagated as where m is the particle's mass. In the above equations, the aerodynamic force acting on a particle is given by F j ad , the electric field force by F j el , and the gravitational force by F g . In detail, the forces read Therein, and gas are the material density of the particle and its surrounding gas. The particle's radius is denoted by r, the particle drag coefficient by C d , the particle velocity relative to the fluid by u rel , the particle's charge by Q, and the gravitational acceleration by g.
The local electric field strength, E j , can be computed through superposition of the electric field of the N particles present in the system, i.e., where = 8.854 × 10 −12 F/m is the electrical permittivity of free space and z(i) j the vector pointing from particle n to i.
The numerical algorithm used to solve these equations is visualized in Fig. 2. As it can be inferred from the flow chart, during the update of the particulate phase it is evaluated for each particle whether it collides with another particle within t or not. To reduce computational expense, collisions are detected via ray casting [26]. In case a collision between particles is predicted, the post-collision positions and velocities are computed analytically. After the update of the position and velocities, the treatment of particle n is finished and it is proceeded with particle n + 1 . Once all N particles are evolved, the computation of time instance t + t is completed and the complete flow system is further propagated by t.

Structure of the adaptive approach
As elaborated in the introduction, in case of charged particles a fixed time-step calculation (as done in the softsphere approach and the method outlined in Sect. 2.1) fails. Reason is that contrary to contact forces, body forces do not only act at the time instance of collision, but affect already substantially the trajectories on which two particle approximate each other. Electrostatic forces, contrary to gravitational forces, change tremendously during the approximation. We tackle this issue by proposing a new computational scheme which is depicted schematically in Fig. 3. To this end, the concept of an outer and inner time-step is introduced. The outer timestep refers to the global time-step of the simulation i.e., the time increment during which the complete flow system is propagated by the time increment t . This includes both the update of the gaseous and the particulate phase. Analogous to the method outlined in the previous section, for each particle it is checked whether it collides with another particle during t or not. In case the particle is not expected to collide, it is propagated as given by Eqs.
(2)- (7). It is noted that the adaptive approach, as well as the standard approach, is limited to binary particle interaction. Thus other particles in the immediate vicinity that may alter the trajectories are not considered. Such a situation occurs, for example, during dilute pneumatic conveying. However, in case the particle collides, for this specific particle an inner time-step, t * , is introduced. As t * is chosen so that t * ≪ t this inner time-step can be interpreted as a local refinement of the time discretization of the particle trajectory whereas the time increment for the propagation of the non-colliding particles is not affected.
The procedure visualized in Fig. 3 involves two parts which require further elaboration in the following text, namely the definition of the size of the inner time step (Sect. 2.3) and the particle propagation (Sect. 2.4) during the inner iteration.

Definition of the inner time-step, t *
In order to capture the interaction of charged particles, the computational time-step is reduced only for the respective particle pair which is in proximity. This so-called inner time-step, t * , is chosen based on several considerations which are governed by the wish to use a value that is sufficiently small to obtain an accurate trajectory and sufficiently large to limit computational expenses. Thus, for each particle pair and each approximation, the ratio between the outer and inner time-steps, s = t∕ t * is determined individually. Let j be the last outer time-step before the possible particle collision which occurs at time instance t j + t c , i.e., t c < t . We introduce two algorithmic bounds to the inner iteration procedure. The first bound to t * limits its lower value by the condition t * < t c . In other words, the inner timestep size needs to be smaller than the time duration between the outer time-step and the collision. Otherwise the approximation would not be resolved and no benefits would be expected from the new approach. The second consideration regards the question which fraction of the outer time-step size is to be resolved through the inner iteration. The limit case of t c → 0 implies that the collision takes place short after t j . In this case also t * → 0 due to its above mentioned lower bound. This would result in a very large number of inner iteration steps to predict the complete outer time-step. Therefore, to avoid excessive computations, the trajectories are resolved symmetrically around the collision point whereas the rest of the outer time-step (i.e., the time span t − 2 t c ) is linearly approximated.
Finally, t * is defined within the given bounds depending on the nature of the interaction. Namely, the time increment for the inner iterations is refined the most for the most relevant case where the particles' kinetic energy ( E kin ) and Coulomb potential ( E pot ) are of similar size, since then they alter each others trajectories the strongest. Equating both quantities yields for the case of particles with the indices 1 and 2 The integral boundaries in Eq. (8a) represent the distance of both particles, |z(1, 2)| , and the minimal distance between monodisperse spheres, 2r. Furthermore u = |u rel,∥ | denotes the particles' absolute relative velocity component parallel to the Coulomb force.
In essence, Eq. (8e) serves as a lower boundary to the inner time-step for the case of u → 0 or q 1 q 2 → ∞ . However, to avoid very small values for t * in the case of u → ∞ or q 1 q 2 → 0 an additional lower boundary of 30 ns is introduced in consistence with [22]. This choice is considered conservative, since [22] modeled even smaller scale interactions then we do in the present study.
To sum up, t * is derived via the relation An overview over the resulting inner time-step is displayed in Fig. 4 depending on the relative velocity and the (equal) charge of two particles. It is noted that according to the figure t * transitions smoothly from one regime to the other.

Particle propagation during inner iteration
For simplification, in the present work only electrostatic forces acting on the particles are considered. Thus, particle acceleration is given by , 30 ns , 1∕t c .
However, the electrostatic field which is required to calculate F j el is only available at the time instances t and t + t but not for the intermediate time-steps during the inner iteration. Thus, the electric field is decomposed into an external field generated by surrounding particle cloud and the local electric field generated by the two collision partners, i.e., The term in the above equation which represents the force between the two particles is fast to solve, thus, it is calculated for each inner iteration. Also, it varies significantly which requires frequent updates to obtain accurate results. The value of E ext , on the other hand, is expensive to compute since it involves the positions of all particles present in the system. Furthermore, it is not expected to vary significantly on a small spatial scale. Therefore, the variation of E ext is approximated linearly between the two known values at the initial positions of both particles.
Having calculated the particle acceleration, its position is evolved using the velocity-Verlet algorithm [27] as follows.
In the above set of equations, (12a) is the first initialization step to get into the new time-step regime whereas the subsequent equations continue the propagation. Compared to the simple Euler-forward, the error of the velocity-Verlet algorithm is smaller at no significant additional computational cost. Further, other non-velocity dependent forces such as Van-der-Waals or polarization forces are easy to include via superposition.

Simulation results and discussion
In the following the behavior of the newly proposed numerical algorithm is evaluated for a variety of conditions. The varied parameters include the particles' size, number density, charge, and initial position. Test cases regarding the collision of two particles are presented in Sect. 3.1 whereas Sect. 3.2 treats the interaction of a large number of particles. In order to keep the set-up simple only electrostatic forces and collisional forces act on the particles. Also, the charge of each particle does not change during collisions. We simulate each condition using three different numerical approaches: A. The standard approach as described in Sect. 2.1. For the computation we apply a constant numerical timestep of t = 0.5 ms since we consider this value to be within the typical order of magnitude for simulations of comparable systems which include the solution of the fluid phase. B. The adaptive approach as described in Sects. 2.2 and 2.3. For the outer iterations, also, a constant numerical time-step of t = 0.5 ms as in A is applied. C. We compute particle trajectories very close to the exact ones by using the standard approach and applying a time-step of t = 0.5 μ s. Through this small time-step the related numerical error diminishes. It is noted that this choice of a time-step is suitable for the herein considered idealized settings but would result in unbearable computational expenses when handling real-scale technical flows.
Thus, approach C delivers a benchmark solution based on which the accuracy of approaches A and B can be established.

Binary particle interaction
First, we investigated the interaction of two particles both of a size of r = 50 μ m, see Fig. 5. Particle 1 is placed with an initial velocity of |u 0 1 | = 1 m/s in x-direction. Particle 2 is initially at rest, i.e., |u 0 2 | = 0 m/s. While the particle charges are constant during each individual simulation, in-between different simulations they were varied systematically from − 0.1 to 0.1 pC. In all considered cases both particles carried an equal absolute charge value, however, their polarity may be opposite, i.e., q 1 = |q 2 | . In other words, our simulation included both attractive and repelling forces.
As it can be seen in Fig. 5, the initial particle distance was varied in x-direction. More precisely, x 0 2 − x 0 1 lies between 1.9 t |u 0 1 | and 3.1 t |u 0 1 | . Due to this choice, the time-span between the discretized time-steps and the actual physical time instance of particle collision changes. In fact, we varied the further up to 4.1 t |u 0 1 | but found that the results are periodic with respect to t |u 0 1 | . Thus, we present the data only for one period. The particle spacing in y and z-direction, on the other hand, was kept constant y 0 1 − y 0 2 = 50 μ m and z 0 1 = z 0 2 = 0 , i.e., a non-central collision was investigated.
To quantify the error of the standard, respectively the adaptive approach, we compare the predicted final particle positions with the exact positions via the relation The error measure in the above equations is normalized by the final absolute velocity, this way it is expressed in units of time-increments. are given on the abscissas. Here, a positive value relates to simulations where both particles are charged with the same polarity and a negative value means that they carry a charge of opposite polarity. The ordinate depicts the initial particle spacing expressed in numerical time-steps it takes for the particles to propagate from their initial location to the point of collision.
The main conclusion to be drawn from Fig. 6 is that over the complete investigated parameter space the adaptive model predicts the final particle position significantly more accurate than the standard model. However, interestingly the error of the new adaptive model is the largest when the initial particle spacing is a whole multiple of the initial velocity times the time-step, i.e., at the ordinate values of 2.0 and 3.0. This is related to the implementation of the new model where the inner iteration is active during one outer iteration (cf. Fig. 2). Thus, if for example the collision happens close to the end of a global time-step, the trajectory of the particle towards the contact point is wellresolved by the inner iteration. The propagation of the particle away from the contact point, on the other hand, may not be equally well resolved if a major part takes place during the following time-step.
To further analyze the model behavior we define an additional error measure based upon the conservation of the system energy which is given by The related error is defined as where E 0 is the sum of the kinetic and potential energy of all N particles present in the system at the beginning of the simulation. An overview over the error field of the adaptive model depending on the particles' charges and their initial spacing (analogous to Fig. 6) is plotted in Fig. 7a. This analysis confirms the general accuracy of the new model, but pronounces more distinctive the error for the cases where the collisions occur close to the beginning or end of a global time-step. As elaborated earlier, this cases result in an incomplete resolution of the steep increase of electrostatic inter-particle forces. Thus the particle is placed unphysically close to the other particle. Due to the high potential energy, the particle is accelerated to a higher kinetic energy prior to the interaction.
Since this setting represents the only visible deficiency of the adaptive model, we further investigated in detail its precise contribution. Therefore, we ran an additional set of simulations where the particle charge was fixed to q 1 = q 2 = 0.1 pC which is according Fig. 7a the worst case in terms of the resulting error size. The initial particle spacing, on the other hand, was varied in very small steps to resolve the resulting error around a critical spacing where particle collision is expected after three time-steps. The obtained values of 2 for 19 simulations resolving the collision time instance from 2.99 to 3.01 t is shown in Fig. 7b. The plot reveals that the error may be up to 6. However, it also ascertains that 2 is only for a very small range of a significant size. Thus, only a very small fraction of particle interactions will be affected. For a multiple-particle system it can be expected that most interactions and, where N j coll denotes the number of collisions between particle up to time step j and N is the number of present particles.
The data plotted in Fig. 9 corresponds to a system of 512 particles each of a diameter of 50 μ m and charge of 0.05 pC. It is noteworthy that the values of 3 of the detailed simulation are very low, more precisely below 0.02%. This fact confirms the assumption taken throughout this paper, namely that the detailed simulation delivers nearly exact results which justifies their usage as benchmark solution. Thus, it can be seen by comparing the curves for N coll of the different models that the adaptive approach delivers results which are approximately four times better than the standard approach for the herein studied conditions. Interestingly, the standard model behaves even worse at approximately t = 16 ms. Here, the kinetic energy explodes instantaneously due to a too tight placement of some particles, which is the same effect that has been already observed for binary particle interaction. However, when comparing the number of collisions it is ascertained that both the standard and adaptive approach deliver satisfactory results. This indicates that the error mainly stems from few particles which do not contribute significantly to the collision frequency.
To get statistical data which is representative for a wide range of condition, we conducted a large number of simulations and analyzed the resulting error distribution. In detail, we varied the size of the 512 particles between 50 and 1000 μ m. Further, we considered between 125 (= 5 3 ) and 1331 (= 11 3 ) particles while keeping their size constant 50 μ m. The resulting pdfs of 4 are depicted in Fig. 10. The shaded areas demonstrate the superiority of the adaptive approach which is clearly substantiated by the distinct peek of the green area at Additionally, in Fig. 10 the normal distributions are plotted which inherit the mean and standard deviations from the shaded areas. According to these approximations, the mean of the standard model is 4 = 0.11 whereas the mean of the adaptive model is located at 4 = 0.017 . Further, the curve of the standard model has a variance of 0.31 compared to a variance of 0.17 of the adaptive approach. These values further demonstrate the improved accuracy of the new model for most investigated conditions. Finally, we evaluate the computational price of the additional inner iteration. To do so, we performed simulations with the standard and adaptive approach considering between 125 and 1728 particles. The resulting computational runtimes when calculating 200 timesteps are plotted in Fig. 11. The observed run-time increase is of the order of O( √ N) which is reasonable considering the improved accuracy.

Conclusions
Collisions between solid particles under the influence of electrostatic forces are of importance in manifold industrial and natural flows. Nevertheless, their precise numerical prediction poses a challenge for development of suitable algorithms. In the present paper, we proposed a new numerical model to compute the interaction between electrically charged particles. Therein, binary particle interaction is calculated through adaptive local refinement of the temporal resolution. The major advantage of this adaptive model is the ability to accurately predict a curved trajectory under strong influence of inter-particle electrostatic forces. The main limitation of the approach is the restriction to binary particle collisions. The new model is of low computational effort and was demonstrated to be superior over the standard model. This model may enhance the accuracy of the prediction of a wide range of charged particle-laden flows subject to an external electric field. In particular, the future implementation of the new approach in our CFD solver [24] will enable us to improve the prediction of the charge accumulation of powder flows significantly. We also anticipate that the new model can be applied to other fields than macroscopic particles, for example to electrolyte solutions [2].