On the accuracy of the numerical computation of the electrostatic forces between charged particles

Powder Technology 322 (2017) 185–194

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

On the accuracy of the numerical computation of the electrostatic forces between charged particles Holger Grosshans* , Miltiadis V. Papalexandris Institute of Mechanics, Materials and Civil Engineering, Université catholique de Louvain, Louvain-la-Neuve 1348, Belgium

A R T I C L E

I N F O

Article history: Received 26 July 2017 Received in revised form 29 August 2017 Accepted 6 September 2017 Available online 11 September 2017 Keywords: Electrostatics Particle dynamics Numerics Coulomb’s law Gauss’ law

A B S T R A C T The accurate calculation of the electrostatic forces between charged particles requires a reliable representation of the emerging electric field. To this end, two principle formulations are available, namely Coulomb’s and Gauss’ laws. While both approaches are in this specific case mathematically equivalent, the numerical solution of Gauss’ law on a grid is contaminated by a spatial discretization error. In this paper, we explore under which conditions this error depreciates the computed results of the particle dynamics. More specifically, we introduce a length-scale of the electrostatic interaction which serves as a base to evaluate the required grid resolution. Also, according to our numerical study, this error reduces in case the particle dynamics is contact dominated. Furthermore, we introduce a computationally efficient hybrid scheme which combines the advantages of both formulations and examine in which situations its accuracy is superior to schemes based on the numerical solution of Gauss’ law. © 2017 Elsevier B.V. All rights reserved.

1. Introduction The electrification of the particulate phase is a frequently occurring phenomenon in flows of solid-fluid mixtures. This may be attributed to various phenomena such as triboelectric charging when dissimilar materials come in contact, induction charging through the presence of an external electric field, and chemical reactions such as combustion leading to ionization. Relevant examples in nature include the electrification of sand storms and of ash erupted by volcanoes [26], as well of aeolian transported grains [18]. Also powders in industrial flows accumulate charge, e.g. when transported pneumatically in pipes [2,9,15]. In many of these cases, electrostatic forces, arising due to the charge of the particles, play an important role on the particle dynamics [24]. Examples include the formation of geological patterns such as razorbacks observed on Mars surface [31] and pollen capture by plants [5]. Furthermore, the electrostatic interaction between particles is utilized in industrial applications, e.g. the separation of different kinds of insulating materials [23] or in electrophotography [29]. For this reason, the accurate computation of electrostatic forces between particles is of great interest. To this end, there are typically

* Corresponding author at: Physikalisch-Technische Bundesanstalt (PTB), 38116 Braunschweig, Germany E-mail address: [email protected] (H. Grosshans).

http://dx.doi.org/10.1016/j.powtec.2017.09.023 0032-5910/© 2017 Elsevier B.V. All rights reserved.

two approaches available. The first one is the estimation of the forces from Coulomb’s law. However, this approach is computationally very expensive. More specifically, it requires the calculation of a force term for each particle pair and, subsequently, their superposition. Thus, the computational effort scales by O(N2 ) with N being the number of particles [22]. The development of Fast Multipole Methods (FMM), initiated by Rokhlin [28], reduced the computational cost to O(NlogN). In fact, FMM is frequently applied to areas such as molecular-dynamics simulations [4,8], but the error constant associated with it is too high. As a result, various computational “tricks” are required to reduce the cost to affordable levels which, however, increase significantly the computational complexity of the method [6]. For this reason, this approach is not well-suited for the computation of charged particle dynamics in real-life situations that consisting of a large number of particles. To tackle this problem, Hogue et al. [17] proposed the so-called screened Coulomb’s force which accounts only for the particles in the proximity of the particle under consideration. They argued that the electrostatic forces between particles can be neglected if their distance is large enough. Thus, they introduced a cut-off distance in their equations while taking care that the cumulative force due to all particles beyond this cut-off were small. A similar approach was subsequently adopted by Korevaar et al. [20] in a LES (Large-Eddy Simulation) framework to simulate the charge build-up of powder in a pipe flow through triboelectric charging. Also, Yang et al. [32] applied a screened Coulomb force term to compute electrostatic effects during particle elutriation in fluidized beds.

186

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

Even though the introduction of any type of cut-off may result in significant computational savings, there is a steep price to pay in terms of accuracy of the calculations. Whereas the computation of collisional forces is of local character, in the sense that it is reasonable (and common practice) to consider collisions between particles of the same or neighboring cells, electrostatic forces are of non-local nature and, therefore, one has to consider the interactions between distant particles too. Moreover, triboelectrically charged particles usually have the same polarity, in which case the effect of the far field on an individual particle does not equal out. The second approach is based on the computation of the electric field from Gauss’ law. According to it, one computes the charge density instead of the charges of individual particles as is required by the first approach. Since the charge density is a field described in Eulerian framework, this approach does not require the calculation of the forces between all pairs of particles. Thus, it is suitable to simulate real-life systems consisting of a large number of particles. For example, Kolniak and Kuczynski [19] used this approach in conjunction with a RANS (Reynolds-Averaged Navier-Stokes) simulation of the gaseous phase to analyze the behavior of charged particles during pneumatic transport operations. In a similar fashion Dastoori et al. [7], Arif et al. [3] and Labair et al. [21] modeled the particle trajectories in electrostatic precipitators and in freefall separators, respectively. Also, the authors of the present paper used Gauss’ law in conjunction with LES and DNS (Direct Numerical Simulations) extensively to model triboelectric charging of particles in pipe [10-12,14,30] and channel [13] flows. We remark that none of the studies provides a rigorous validation of the numerical implementation of the source term for the electrical forces. However, the numerical solution of Gauss’ law requires the discretization of the electric field on the computational mesh. Since the forces between particles increase with the inverse square of their distance, they are especially important when the particles are close to each other. As an example, let us assume a cloud of particle of a size of 25 microns. For the electric field shall to be accurately resolved while particles are close to each other, the grid resolution should at least be of the order of the particle size. This means that a cubic domain of 1 cm3 would require 54 million cells assuming an equidistant grid. From this example it becomes evident that the required resolution for practical applications, whose dimensions are at the order of meters, is simply not attainable. The criteria for required resolution in certain areas of computational science and engineering years have long been established. A typical example is fluid turbulence, for which the smallest occurring length-scale, the Kolmogorov scale (see Pope [27] for an overview over the topic), provides a firm lower limit for grid resolution. Thus, a simulation produces fully resolved results if the grid size is smaller than this scale. On the other hand, such firm limits are hard to define in electrostatics because the amplitude of the Coulomb forces increases as the particles move closer to each other. This problem was already recognized by Aboud et al. [1] for another field of physics, namely the electrostatic interactions between ions in electrolytic solutions. Thus, they proposed a P3 M algorithm [16] where the long-range components of the forces are computed based on Gauss’ law and the short-range components based on Coulomb’s law. Therein, the double counting of charges due the overlap of both ranges was corrected through an additional force term, the so-called reference force. Despite its similarity to the situation considered in the present paper, the dynamics of ions exhibits some important differences: Due to their atomic size, van der Waals forces were required to be taken into account by Aboud et al. [1] for the short range. Moreover, the size and the charge of ions of the same type is equal. Both is not the case for solid particles which are much larger, polydisperse and possibly of a strongly homogeneous charge distribution. Besides leading to qualitatively different particle interactions, this prevents the precomputation and tabulation the

reference force which was required by Hockney and Eastwood [16] and Aboud et al. [1] to reduce the computational burden. Also, it is worth to mention that P3 M, as well as the above mentioned FMM is cumbersome to implement in non-periodic domains. The main objective of this paper is twofold. First, to analyze the accuracy of the numerical approach for electrostatic particle interaction based on Gauss’ law. Second, to propose an alternative hybrid approach according to which the interaction between neighboring particles is computed via Coulomb’s law and the far-field effects are computed via Gauss’ law. Whereas this hybrid approach was inspired by the work of Aboud et al. [1] it will be shown in the following that the above mentioned differences has some important implications. The paper is organized as follows. In Section 2 we present the governing equations and describe the numerical approaches. In Section 3 we provide an analysis of the length-scales involved in the electrostatic interactions between moving particles. In Section 4 we present and discuss our numerical results and compare the accuracies of the approach based on Gauss’ law and the proposed hybrid approach. Finally, in Section 5 we summarize the conclusions of our study. 2. Governing equations and numerical methods The acceleration of an electrically charged particle in vacuum and in absence of gravitational forces is given by dup = f coll + f el , dt

(1)

where up is the particle velocity, fcoll is the term for collisional accelerations, and fel is the acceleration due to an electric field. The above equation can also be considered a part of an implementation in a computational fluid dynamics solver. Typically, the complete solver requires the solution of the Navier-Stokes equations for the carrier gas. Based on the local characteristics of the gaseous phase further terms representing e.g. drag, virtual mass, Basset, Magnus or Saffman forces are added to Eq. (1) [25]. The term for collisional accelerations accounts for both particleparticle and particle-wall collisions. In the present work we consider fully elastic and binary collisions. For the exact formulation of this term the reader is referred to Grosshans and Papalexandris [13]. The force acting on a particle due to the presence of an electric field, E, is given by f el =

QE , mp

(2)

where Q is the electric charge of the particle and mp its mass. According to Gauss’ law, the electric charge density per unit volume is related to the divergence of the electric field by ∇ •E =

qel . e

(3)

Therein, e = 8.854 F/m denotes the electrical permittivity of the vacuum and qel is the electric charge density. The integral of qel over a control volume V that contains n particles is equal to the sum of the charges of the n particles,  V

qel dV =

n 

Qi .

(4)

i=1

Further, the electric field strength is the gradient of the electric potential v, E = −∇v .

(5)

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

For the case of charged particles, the equivalence of Gauss’ and Coulomb’s laws is established in a straightforward manner as follows. Upon integration of Gauss’ law (Eq. (3)) over a surface S, and by the virtue of the divergence theorem, we obtain  E • dS = S

q . e

(6)

In the above equation, S is a surface enclosing any volume and q is the charge located within this volume. In case of describing the electric field around a charged particle, q becomes the charge carried by the particle, Q. Further, using the spherical symmetry of the electric field surrounding a uniformly charged, spherical particle, allows to solve the integrand. Assuming S to be a spherical surface yields 4 p |r|2

Q r • E(r) = , |r| e

(7)

where r is a vector pointing from the particle center to any point on S. By rearranging the terms of this equation, we obtain Coulomb’s law, E(r) =

Qr . 4 p e |r|3

(8)

In other words, Gauss’ and Coulomb’s law are mathematically equivalent. However, even though Eqs. (3) and (8) are equivalent, the numerical methods for solving them when dealing with a large number of particles may differ considerably in terms of accuracy and computational costs. The computation of the electrostatic forces based on Gauss’ law is computationally more efficient than the direct calculation and superposition of Coulomb’s law. Since Eq. (3) is a field equations, it is formulated in Eulerian framework and needs to be solved once per time-step for the domain under consideration. Thus, the computational costs scales with the physical dimensions of the computational domain and the required grid resolution but, crucially, is independent of the number of particles. On the contrary, Eq. (8) needs to be solved for each binary particle interaction. Thus, for a system of N particles, the number of interactions per particle that need to be computed is N − 1 and the total number of interactions raises to N(N − 1)/2 [22]. If N is small enough this approach is feasible. But since the number of operations scales quadratically with N, the computational cost for real-world problems is prohibitively high. Concerning the accuracy of the numerical solution, we remark that Eq. (8) does not involve any derivatives. Thus, the solution of Coulomb’s law requires only the discretization of the temporal derivative when Eq. (8) is combined with Eq. (1). The contribution of the truncation error resulting from the usage of a finite-order discretization scheme is relatively easy to control through the choice of a sufficiently small time-step. This means that this is a mesh-free approach. Consequently, the computational cost scales only with the number of particles only and does not depend on the size of the domain. On the other hand, Eq. (3) does involve a spatial derivative. Therefore, it needs to be solved on a computational mesh (grid) and its solution is affected by the truncation error of the discretization scheme; for example, in the present study we have chosen a secondorder central difference scheme which is solved using multigrid acceleration. This fact has some important consequences. More specifically, upon inspection of Eq. (8), it is evident that the electric field strength between particles increases as neighboring particles get closer to each other. Thus, one cannot generally assume that mesh refinement will result in numerically converged results. In particular, the electrostatic attractive or repelling forces between particles that temporarily find themselves in the same computational cell are not fully accounted for. This problem is inherent in the

187

coupling of the Eulerian and Lagrangian frameworks and explains why this approach is limited to small numbers of particles or, in the case of particle-laden flows, dilute mixtures. Following the above discussion we propose a new hybrid scheme which combines both approaches. As regards the interaction of an individual particle with the space charge, i.e. with the electric field originating from the particle cloud, Gauss’ law is accurate. Thus, it is employed to compute this part of the influence on the particle. Furthermore, a sub-grid model is derived from the idea to compute the forces between particles present in the same grid cell based on Coulomb’s law. In detail, the algorithm compute the electrostatic forces on a single particle i follows the following steps: 1. Derive the Eulerian charge density field from the Lagrangian particle charges (Eq. (4)). 2. Calculate the electric field at the location of particle i based on Gauss’ law (Eq. (3)). 3. Identify the n particles which reside in the same computational cell as particle i. 4. Solve the Coulomb interaction (Eq. (8)) between particle i and each of the n particles. Sum up the resulting electric field at the location of particle i. 5. Sum up the electric field obtained through steps 2 and 4 and compute the force on the particle (Eq. (2)) based on this sum. To sum up, there are three schemes available to compute electrostatic forces on particles, namely Coulomb’s and Gauss’ laws and the newly proposed hybrid approach. The numerical accuracy of their solution is explored in the following sections. 3. Length-scale analysis In this section we analyze the length-scales involved in the interactions between charged particles. This analysis provides insight as to in which situations the differences between the numerical approaches presented are important and in which situations these differences are negligible. For simplicity purposes, our analysis considers only the case of identical particles with the same polarity. i.e. the electrostatic forces between them are of repelling nature. To this end, we examine the balance between the kinetic energy of a particle and the electrostatic potential due to a neighboring particle. The kinetic energy reads Wkin =

1 2 mp u20 = p qp rp3 u20 . 2 3

(9)

The potential energy of the electric field arising from the charge of a neighboring particle is given by Wpot =

Q1 Q2 . 4pel

(10)

In the above equations, qp , rp , u0 and Q1/2 are the density, radius, velocity and charge of both particles and l is the distance between the mass centers of the two particles. By setting Wkin = Wpot and by using the above relation, we can determine at which particle distance the two energies are equal. This distance, henceforth denoted by k, represents a characteristic length-scale relative to the dynamics of the specific particle. In view of the above relations for the kinetic and potential energies, k is given by the following expression: k=

3 Q1 Q2 . 8 p2 qp u20 e rp3

(11)

The basic meaning of this length-scale is that the kinematics of interaction between particles depends only on k and not on their

188

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

individual properties. In other words, for particle interactions of an equal k the particle trajectories x(t) are similar. This analysis neglects the influence of the initial distance, a, of the particles, which is in general reasonable if a  k, i.e. if the particles in their initial position are little influenced by the electric field surrounding the other particle. In the limiting case of very high values of k, e.g. in the case of light but strongly charged particles carry, it is natural to assume that the repelling electrical forces are sufficiently strong to prevent the particles from getting close to one another. In this case the particle trajectories are strongly influenced by electrostatic forces. As a consequence, inter-particle collisions and subsequent inter-particle charge transfer [13] do not occur. On the other hand, if k is small, the particles may come very close to each other and might even collide if their radii are sufficiently large. To relate these two lengths, we define a new non-dimensional number, namely w = 2 rp /k .

(12)

In order to demonstrate the physical relevance of this quantity, let us consider the limiting case of heavy particles carrying a low charge. This results in a small value for k which means that the particle dynamics is dominated by its own inertia. The electric forces, on the other hand, play a minor role even during close proximity with other charged particles. If in this case the particles are of a large size, resulting in a high value for w, the particle interactions are dominated by collisions between them. On the other hand, if the particles are small, i.e. w is small, the particles come close to each other and interact strongly on a small scale but do not collide with each other. To sum up, we expect the electric forces to be important for particle interaction dynamics if both k and w are sufficiently small. On the contrary, if w is large inertia and collisions are dominating whereas for large k the electric forces prevent particle approaches. In the present study, we are most interested in particle interactions where both k and w are relatively small. Since for these parameters the most important interactions take place on a small scale, the solution of Eq. (3) requires a very fine numerical grid. To evaluate the quality of the grid resolution, h, we define a second non-dimensional parameter, namely f = h/k .

(13)

According to this definition, small values for f represent a well resolved computation while large values relate to a coarse grid. We remark that whereas w is a physical non-dimensional parameter related to charged particle dynamics, f is a numerical nondimensional parameter related to the appropriateness of the grid resolution. We analyzed the performance of the various numerical approaches mentioned above for different values for f through numerical studies the results of which are discussed below.

4. Numerical results and discussion We performed a series of computations concerning charged particles in order to test the accuracy of the numerical methods based on Coulomb’s and Gauss’ laws, as well as the accuracy of the proposed hybrid scheme. Further, we investigated the suitability of the scaling parameters k, w and f, proposed in the previous section, to characterize the particle interaction. In the first part, we studied binary particle interactions whereas in the second part a larger number of particles is computed. An overview over the detailed parameters of the computations are given in Table 1. 4.1. Binary particle interaction In the binary particle computations we considered two particles that are initially placed 10 mm apart. Their initial velocity is u0 and their velocity vectors are pointing toward one another. Thus, if no electrostatic forces were considered all cases would result in head-on collisions. The hypothetical contact point is in the center of computational domain, which is an orthogonal parallelepipedon with dimensions 4cm×2cm×2cm in which the axis of the longest dimension is collinear to the velocity vector. This axis is denoted by the coordinate x. Since the initial velocity vectors are parallel and no additional forces (such as gravity) are taken into account, the particle kinematics are essentially one-dimensional. Eq. (3) was computed on a Cartesian, orthogonal, equidistant mesh. As discussed above, solving the equation of motion of a particle via Coulomb’s law requires no computational mesh, i.e. spatial discretization of the computational domain. Evidently, the numerical solution contains a truncation error related to the temporal discretization. However, this error can be controlled without excessive computational cost. In order to determine the appropriate temporal discretization, we computed case 1 with four different computation steps. The computed histories of the distance between the particles and the center of the domain are visualized in Fig. 1. Since the trajectories of both particles are symmetric with respect to the center of the domain, located at the origin x = 0, in all figures related to binary particle interaction we provide the results for one particle only. In Fig. 1, and all figures of particle trajectories, the abscissa gives the time multiplied by the initial particle velocity, whereas the ordinate gives the respective particle position along the x-axis. In this figure it can be observed that the particles are transported due to their initial kinetic energy toward the symmetry line. Since it resides within the electric field of the other particle, its kinetic energy is transformed into potential electrostatic energy as both particles approach each other. Consequently, the particle is slowing down until its velocity becomes zero at about u0 t = 5.4 mm. Afterward, it accelerates in the opposite direction and the potential energy is transformed back to kinetic energy. Fig. 1 shows that the computation using a step of u0 Dt = 100lm, or smaller, gives accurate results. Based on the results of this convergence study, a

Table 1 Parameters of all simulated cases. Note that the cell size h and its normalized value f are only relevant for the computations of Gauss’ law and the proposed hybrid approach whereas solving Coulomb’s law is mesh-free. The parameters k, w and f are derived from the other given quantities. Case

k/mm

w

f

Q/pC

qp /(kg/m3 )

rp /lm

u0 /(m/s)

h/mm

1 2 3 4 5 6 7 8 9 10 11

2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75 2.75

0.18 0.18 0.18 1.46 1.46 0.18 0.36 0.73 0.18 0.18 0.18

0.73 0.73 0.73 0.73 0.73 0.36 0.73 0.73 0.18 1.46 2.91

100 50 200 1600 800 100 200 400 100 100 100

4000 1000 4000 2000 500 4000 2000 1000 4000 4000 4000

250 250 250 2000 2000 250 500 1000 250 250 250

0.5 0.5 1.0 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5

2 2 2 2 2 1 2 2 0.5 4 8

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

15

u0 t = 200 m u0 t = 100 m u0 t = 50 m u0 t = 25 m

12 x / mm

189

9 6 3 0 0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 1. Particle trajectory in binary particle interaction for case 1, q = 100 pC, qp = 4000 kg/m3 , rp = 250lm, u0 = 0.5 m/s, computed using Coulomb’s law. The results indicate that using a computation step of u0 Dt = 100lm or smaller gives converged results.

computation step of u0 Dt = 50 lm has been used in the numerical tests that we present below. Consequently, the computations based on Coulomb’s law with this stepping are accurate enough to be used as the reference (“exact”) results for the purposes of testing the accuracy of the other schemes. First, we investigate the role of the non-dimensional parameter w based on the direct solution of Coulomb’s law with the aforementioned computation step. The resulting particle trajectories of the cases 1–5 are visualized in Fig. 2. It can be observed that the trajectories of cases 1, 2 and 3 are identical. The value of w for these two cases is the same, even though the particle charge, radius and initial velocity are different (cf. Table 1). The same conclusion can be drawn when comparing the trajectories of cases 4 and 5. These results suggest that w is in fact the appropriate parameter to characterize the nature of interaction between the charged particles. For cases 4 and 5 the parameter w was increased from 0.18 to 1.46 by a corresponding increase of the particle radius. As anticipated in the previous section, this leads to a collision between the particles. This collision is manifested by the sharp change of the trajectories at u0 t = 3.2 mm., where the particle distance from the origin equals its radius. Having verified the appropriateness of the non-dimensional parameter w for comparing the length-scales involved in binary particle interactions, we proceed to examine the accuracy of the different numerical approaches. The particle trajectories for case 1 are visualized in Fig. 3. According to the trajectories calculated from Coulomb’s law, which are the reference results, show that the particles approach each other due to their initial momentum.

16

However, when they come sufficiently close, then the repelling electrostatic forces become strong and force the particles to slow down. It can be observed that the trajectory predicted by the numerical solution of Gauss’ law is qualitatively different from the reference (exact) result based on Coulomb’s law. The explanation for this important error is the following When the particles come close to one another, the electric field between them is discretized by only a few cells. This means that as the distance between the particles decreases, the electric field between them becomes more and more under-resolved. Consequently, the errors due to spatial discretization become more important and contaminate the numerical results. Eventually, the two particles will be located in the same computational cell, which does not allow to account for electrostatic forces at all. This results in the prediction of an unphysical collision at u0 t = 6.0 mm when using Gauss’ law. On the other hand, the results of the hybrid approach that we propose are very similar to those of Coulomb’s law. Therein, the electrostatic forces are computed via Coulomb’s law when the particles are located in the same cell (cf. Section 2). Nonetheless, the hybrid approach is not completely free of errors and, therefore, the trajectory predicted by it does not exactly match the reference trajectory that is computed by Coulomb’s law. The small discrepancy observable in Fig. 3 relates to the implementation of the hybrid approach. More specifically, if both particles are close to each other but located in an adjacent cell they reside in an area of a bad resolved electric field. Nevertheless, if the grid is refined, see the curve for f = 0.36 in Fig. 3, also the solution of Gauss’ law gives accurate results. Based

cases 1, 2 and 3, χ = 0.18 cases 4 and 5, χ = 1.46

x / mm

12

8

4

0 0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 2. Particle trajectories in binary particle interaction computed via Coulomb’s law. Cases 1, 2 and 3 are of the same scaling parameters and result in the same particle trajectories. For cases 4 and 5 the increase in w leads to a collision between both particles at u0 t = 3.2 mm.

190

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

12

Reference solution, case 1 Approach based on Gauss’ law, case 1, = 0.73 Approach based on Gauss’ law, case 6, = 0.36 Hybrid approach, case 1, = 0.73

x / mm

9

6

3

0 0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 3. For case 1 the hybrid approach gives accurate results, whereas the solution of Gauss’ law results even in an unphysical collision. The Gauss’ law solution can be recovered by refining the grid to f = 0.36.

on this resolution no collision between the particles is predicted any longer. It is noteworthy, that for these conditions the hybrid approach gives a similar accuracy than Gauss’ law, even if the applied grid is only half as fine, resulting in a ratio of the number of cells of 1:8. Thus, the hybrid approach provides a computationally cheap alternative to approaches based on Coulomb’s or Gauss’ law which is an important outcome of the present study. The previous example elucidates the differences in the three approaches and indicates that the approach based on Gauss’ law can produce spurious results unless one employs very refined computational meshes. For this reason, it is deemed useful to explore further how the numerical errors of this approach and of the hybrid method scale with respect to the non-dimensional parameters w and f. We start with the scaling of the errors with respect to w. As a measure of the numerical errors of the two approaches, we employ the root-mean-square (rms) of the error of the particle locations at each time instance, i.e. u

0t

erms =

0

|xC − x|2 dt u0 t

1/2

.

(14)

Herein, xC is the “exact” solution, i.e. the one computed from Coulomb’s law and x is the solution computed from another approach. The overall computation length u0 t, was 16 mm. As it can be inferred from Figs. 2 and 3, this computation length covers the

complete interaction between the two particles, i.e. their approach, proximity and separation. In Fig. 4, this error is plotted for a fixed f of 0.73 and varying w. The data points corresponding to w = 1.46 represent cases 4 and 5. The fact that these cases yield the same error even though they correspond to different dimensional parameters demonstrates that err depends directly on the non-dimensional groups w and f. Further, a filled symbol in Fig. 4 denotes that a collision was predicted by the respective numerical method whereas an open symbol means that no collision was predicted. By comparison with the results from the solution of Coulomb’s law, it was identified whether a predicted collision was physical or spurious. If the collision was predicted correctly, the symbol is filled black, if not, it is filled red. It is inferred from Fig. 4 that non-physical collisions occurred when using Gauss’ law for w = 0.18 and 0.36. On the contrary, the hybrid approach predicts accurately that no collisions take place for these values of w. Consequently, the error of the computation is much smaller compared to Gauss’ law. Moreover, since the scaling parameters are identical and the kinematics of the particle interaction is not disturbed by collisions, the hybrid approach gives the same errors for the cases with w = 0.18 and 0.36. However, in both approaches the error drops fast when the predicted collision is physical. A large value of w means that the particle dynamics is collision dominated and, thus, the resolution requirements for an accurate representation of the electric field are not too strict. It should be clarified that the term “collision dominated”

1

Approach based on Gauss’ law Hybrid approach

err / mm

0.8 0.6 0.4 0.2 0 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

Fig. 4. Numerical errors of the hybrid approach and the approach based on Gauss’ law for cases with fixed f = 0.73 and variable w. The data points correspond to the cases 1, 7, 8, 4 and 5, out of which the latter two give an identical result. The filled symbols represent cases where a collision occurs, the empty symbols where no collision occurs. Symbols filled with a red color mark cases where a collision was erroneous predicted by the respective method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

1

191

Approach based on Gauss’ law Hybrid approach

err / mm

0.8 0.6 0.4 0.2 0 0

0.5

1

1.5

2

2.5

3

Fig. 5. Numerical errors of the hybrid approach and the approach based on Gauss’ law for cases with w = 0.18 and variable f. The data points correspond to the cases 9, 6, 1, 10 and 11, respectively. The filled symbols represent cases where a collision occurs, the empty symbols where no collision occurs. Symbols filled with a red color mark cases where a collision was erroneous predicted by the respective method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

12

Reference solution Approach based on Gauss’ law Hybrid approach

x / mm

9

6

3

0

0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 6. Particle trajectories of case 6 (w = 0.18, f = 0.36). The solution of Gauss’ law gives slightly better results than the hybrid approach.

does not imply that collisions take place at high frequencies, since this would require high particle concentration and Stokes number. Instead, in this context the term implies that particle trajectories are influenced mostly by particle collisions than by electrostatic interactions between particles. The fact that collisions also happen for w < 1 can be explained by the small initial distance between particles so that the initial forcing due to electric field is considerable. To sum up, the numerical errors are significantly reduced when the hybrid approach is used, especially in the case of a low w. In Fig. 5, we have plotted the rms error for cases with fixed w = 0.18 but variable f. In other words, the physical parameters are kept the same but the grid resolution varies. The smallest value of f = 0.18 corresponds to the highest resolution. In this case, the mesh spacing was smaller than the minimum distance between particles, 2rp . This means that the electric field is very well resolved and that the approach based on Gauss’ law gives accurate results. Moreover, due to the small cell size, the two particles can never be located in the same cell. Also, at such high resolutions, the hybrid approach reduces to the approach based on Gauss’ law. For this reason both schemes produce the same numerical error. The next coarser grid, f = 0.36, still provides an adequate resolution of the electric field. For this reason, the errors of both approaches are still quite small, albeit higher than in the previous case with f = 0.18. Moreover, for this case the error of the

hybrid approach is slightly higher than the one produced by the approach based on Gauss’ law. This can be explained by the fact that in this particular case and according to the hybrid approach, the computations of Coulomb’s law are superimposed to the predictions from a well-resolved electric field. This essentially results in an over-prediction of the electrostatic forces which is an inherent source of error in the hybrid approach due to the double counting of forces. The effect of this superposition can be seen in the particle trajectories visualized in Fig. 6. However, the error is obviously small and reduces considerably further at larger values of f for which the electric field becomes under-resolved. We may therefore conclude that for very fine grids the numerical solution of Gauss’ law and the hybrid approach give the same results, whereas the hybrid approach produces more accurate results at coarser grids. Nevertheless, there is a small range of in-between grid sizes where the solution of Gauss’ law gives slightly more accurate results than the hybrid approach. The increase of f to 0.73 leads to an under-resolved electric field. As a result, the approach based on Gauss’ law predicts a spurious collision between particles. On the other hand, our hybrid approach gives clearly more accurate results. For a further increase of f to 1.46 the numerical errors in the computation of the electric field become higher. However, according to Fig. 5 the error of the numerical solution of Gauss’ law is small and

192

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

12

Reference solution Approach based on Gauss’ law Hybrid approach

x / mm

9

6

3

0 0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 7. Particle trajectories of case 10 (k = 2.75 mm, w = 0.18, f = 1.46). The non-physical collision predicted based on Gauss’ law corrects the particles’ momentum balance and leads to an accurate trajectory.

only slightly higher than the hybrid approach. The underlying reason can be observed in Fig. 7. More specifically, due to the under-resolved electric field Gauss’ law predicts once again a spurious (non-physical) collision. This collision gives an additional momentum to the particle and forces it to correct its trajectory. Since the momentum exchange due to the collision is similar to the momentum loss due to the low resolution, the resulting error is low. However, this is considered a ‘lucky shot’ and not a good performance of Gauss’ law. This is confirmed by the results for f = 2.93 where the even worse resolved electric field also leads to a non-physical collision and consequently to a high error. Interestingly, the larger f is the smaller is the error of the hybrid approach. In fact, one can assume that for further increases of f the results of the hybrid approach converge asymptotically to Coulomb’s law.

N = 36 particles were distributed with an equidistant spacing of 1 cm. All particles were assigned the same velocity amplitude but their initial direction of motion was random. Identical initial conditions were used and solved numerically via the three approaches that are examined in this paper. In the right image of Fig. 8 visualizes the particle positions predicted by all three schemes after u0 t = 10 mm. As regards the positions of the particles, it can be observed that the predictions of the three numerical approaches are quite close. However, the predictions of the hybrid approach (green) are more closer to the ones by Coulomb’s law (white) than the ones by Gauss’ law (red). In order to evaluate the error quantitatively, compared to Eq. (14) a slightly different definition of the error was used, namely

4.2. Multiple particle interaction

Erms =



N 2 1   xC,i − xi  N

1/2

.

(15)

i=1

We now proceed to examine the accuracy of the numerical methods to systems with many particles. In such systems, each particle does not interact only with its closest neighbor and it is under the influence of the electric field generated by the entire particle cloud. The numerical set-up and the initial conditions of case 1 is visualized in Fig. 8 (left). The computational domain is a cubic box of the size 8cm × 8cm × 8cm. However, all particles were initially placed in the x-z-plane, so the problem reduces to two dimensions. In total,

Thus, Erms represents the sum of the rms errors of the trajectories of all particles at a certain time instance. Since this is an additive quantity, Erms grows with computation length. In fact, a N-body problem as considered here constitutes a chaotic system. Therefore, due to the contribution of any error the locations of the particles predicted by the three approaches will not correlate to each other anymore after a certain amount of time. At this stage, Erms is expected

Fig. 8. Left: Initial conditions of case 1; right: case 11 after u0 t = 10 mm. The white spheres correspond to the solution of Coulomb’s law, the red to Gauss’ law and the green to the hybrid approach. The particles are all of the same size but have been enlarged in these figures to different sizes for visualization purposes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

6

Approach based on Gauss’ law, case 1, Hybrid approach, case 1, Approach based on Gauss’ law, case 11, Hybrid approach, case 11,

5

err / mm

4

193

= 0.73 = 0.73 = 2.91 = 2.91

3 2 1 0 0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 9. Evolution of the numerical error of the cases 1 and 11 (k = 2.75 mm, w = 0.18). A finer grid leads to more accurate results, whereas especially for f = 2.91 the hybrid approach is superior to Gauss’ law.

to obtain a value of approximately half the domain size which we confirmed by a long term computation of case 1. Thus, in order to evaluate the accuracy of a scheme, we compare the transient of the error at the initial stage of the computation. The resulting errors of cases 1 and 11 is plotted in Fig. 9. The curves confirm the conclusion of the binary particle studies that for these conditions the hybrid approach is more accurate than Gauss’ law. Moreover, it is demonstrated how the usage of a finer grid improves the accuracy of the numerical results. When evaluating the errors, one should keep in mind that case 1 requires 64 times more cells than case 11 when computed in a three dimensional domain. Nevertheless, the hybrid approach shows a relatively high error for case 11 compared to the binary particle study, cf. Fig. 5. This is related to the limited resolution of the electric field generated by the surrounding particles which does not exist if only two particles are present. However, also for this conditions the hybrid approach is an obvious improvement to Gauss’ law. Further, the errors related to the computation of case 5 with multiple particles is plotted in Fig. 10. As discussed above, due to its high w-value this case is collision dominated. Since the size of the individual particles is larger than a cell, the hybrid approach and Gauss’ law give identical results. The small errors given in Fig. 10 confirm the conclusion from the binary particle study that if w is high it is sufficient to compute Gauss’ law on a relatively coarse grid.

particles, namely Coulomb’s and Gauss’ laws. An important conclusion is that while being mathematically equivalent their numerical solution is inherently different since the latter is disturbed by spatial discretization errors. We found that if the particle size is large compared to the scale of electrostatic interaction, collisions dominate particle dynamics and Gauss’ law gives accurate results. Moreover, if the affordable cell size is small compared this scale the electric field is well-resolved and the error is small. On the contrary, if the particle size is small and the grid coarse, the imprecise solution of Gauss’ law significantly affects the predicted particle trajectories. For these cases we proposed, based on the P3 M algorithm, a new hybrid scheme which combines the advantages of both formulations. We demonstrated, that compared to the approach based on Gauss’ law the same accuracy can be obtained whereas the computational costs are reduced by a factor of eight. Generally speaking, the new scheme is computationally more efficient than solving Coulomb’s law for a large number of particles and for certain conditions its accuracy is superior to Gauss’ law. The presented study is limited to particles charged with an equal polarity, thus, to repelling forces. The future consideration of attractive forces will extend the applicability of the study also to the analogous fields of physics such as magnetic and gravitational interaction.

Acknowledgments 5. Conclusions In this paper we studied the two commonly applied approaches when computing the dynamic interaction of electrically charged

1.4

The first author gratefully acknowledges the financial support of the National Research Fund of Belgium (FNRS) under the GRANMIX Projet de Recherche grant.

Approach based on Gauss’ law Hybrid approach

1.2

err / mm

1 0.8 0.6 0.4 0.2 0 0

2

4

6

8 u0 t / mm

10

12

14

16

Fig. 10. Evolution of the numerical error of the case 5 (k = 2.75 mm, w = 1.46, f = 0.73). The particle dynamics is collision-dominated and the effect of the electrostatic forces is small. Therefore, both approaches give accurate results.

194

H. Grosshans, M. Papalexandris / Powder Technology 322 (2017) 185–194

References [1] S. Aboud, D. Marreiro, M. Saraniti, A poisson P3 M force field scheme for particle-based simulations of ionic liquids, J. Comput. Electron. 3 (2004) 117–133. [2] I. Adhiwidjaja, S. Matsusaka, S. Yabe, H. Masuda, Simultaneous phenomenon of particle deposition and reentrainment in charged aerosol flow - effects of particle charge and external electric field on the deposition layer, Adv. Powder Technol. 11 (2) (2000) 221–233. [3] S. Arif, D.J. Branken, R.C. Everson, H.W.J.P. Neomagus, L.A. le Grange, A. Arif, CFD modeling of particle charging and collection in electrostatic precipitators, J. Electrostat. 84 (2016) 10–22. [4] J.A. Board, J.W. Causey, J.F. Leathrum, A. Windemuth, K. Schulten, Accelerated molecular dynamics simulation with the parallel fast multipole method, Chem. Phys. Let. 198 (1992) 89–94. [5] G.E. Bowker, H.C. Crenshaw, Electrostatic forces in wind-pollination-part 1: measurement of the electrostatic charge on pollen, Atmos. Environ. 41 (2007) 1587–1595. [6] E. Darve, The fast multipole method: numerical implementation, J. Comput. Phys. 160 (2000) 195–240. [7] K. Dastoori, B. Makin, M. Kolhe, M. Des-Roseaux, M. Conneely, CFD modelling of flue gas particulates in a biomass fired stove with electrostatic precipitation, J. Electrostat. 71 (2013) 351–356. [8] H.Q. Ding, N. Karasawa, W.A. Goddard, III, Atomic level simulations on a million particles: the cell multipole method for Coulomb and London nonbond interactions, J. Chem. Phys. 97 (1992) 4309–4315. [9] W. Fath, C. Blum, M. Glor, C.D. Walther, Electrostatic ignition hazards due to pneumatic transport of flammable powders through insulating or dissipative tubes and hoses - new experiments and calculations, J. Electrostat. 71 (3) (2013) 377–382. [10] H. Grosshans, M.V. Papalexandris, Evaluation of the parameters influencing electrostatic charging of powder in a pipe flow, J. Loss Prev. Process Ind. 43 (2016) 83–91. [11] H. Grosshans, M.V. Papalexandris, Large eddy simulation of triboelectric charging in pneumatic powder transport, Powder Technol. (2016) 1008–1015. [12] H. Grosshans, M.V. Papalexandris, A model for the non-uniform contact charging of particles, Powder Technol. 305 (2016) 518–527. [13] H. Grosshans, M.V. Papalexandris, Direct numerical simulation of triboelectric charging in a particle-laden turbulent channel flow, J. Fluid Mech. 818 (2017) 465–491. [14] H. Grosshans, M.V. Papalexandris, Numerical study of the influence of the powder and pipe properties on electrical charging during pneumatic conveying, Powder Technol. 315 (2017) 227–235.

[15] J. Guardiola, V. Rojo, G. Ramos, Influence of particle size, fluidization velocity and relative humidity on fluidized bed electrostatics, J. Electrostat. 37 (1–2) (1996) 1–20. [16] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles, Adam Hilger. 1988. [17] M. Hogue, C.I. Calle, P.S. Weitzman, D.R. Curry, Calculating the trajectories of triboelectrically charged particles using discrete element modeling (DEM), J. Electrostat. 66 (2008) 32–38. [18] A.K. Kamra, Physical sciences: visual observation of electric sparks on gypsum dunes, Nature 240 (1972) 143–144. [19] P.Z. Kolniak, R. Kuczynski, Numerical modeling of powder electrification in pneumatic transport, J. Electrostat. 23 (1989) 421–430. [20] M.W. Korevaar, J.T. Padding, M.A. van der Hoef, J.A.M. Kuipers, Integrated DEM-CFD modeling of the contact charging of pneumatically conveyed powders, Powder Technol. 258 (2014) 144–156. [21] H. Labair, S. Touhami, A. Tilmatine, S. Hadjeri, K. Medles, L. Dascalescu, Study of charged particles trajectories in free-fall electrostatic separators, J. Electrostat. 88 (2017) 10–14. [22] G. Lapenta, Particle in Cell Methods With Application to Simulations in Space Weather, Citeseer, Astrofysica, Centrum Voor Plasma. 2013. [23] J.C. Laurentie, P. Traoré, L. Dascalescu, Discrete element modeling of triboelectric charging of insulating materials in vibrated granular beds, J. Electrostat. 71 (6) (2013) 951–957. [24] S. Matsusaka, H. Maruyama, T. Matsuyama, M. Ghadiri, Triboelectric charging of powders: a review, Chem. Eng. Sci. 65 (2010) 5781–5807. [25] M.R. Maxey, J.J. Riley, Equation of motion for a small rigid sphere in a nonuniform flow, Phys. Fluids 26 (1983) 883–889. [26] T. Miura, T. Koyaguchi, Y. Tanaka, Measurements of electric charge distribution in volcanic plumes at Sakurajima volcano, Japan, Bull. Vulcanol. 64 (2002) 75–93. [27] S.B. Pope, Turbulent Flows, Cambridge University Press. 2000. [28] V. Rokhlin, Rapid solution of integral equations of scattering theory in two dimensions, J. Comput. Phys. 86 (1990) 414–439. [29] L.B. Schein, Recent advances in our understanding of toner charging, J. Electrostat. 46 (1999) 29–36. [30] N. Schwindt, U. von Pidoll, D. Markus, U. Klausmeyer, M.V. Papalexandris, H. Grosshans, Measurement of electrostatic charging during pneumatic conveying of powders, J. Loss Prev. Process Ind. (2017)In press. [31] T. Shinbrot, K. LaMarche, B.J. Glasser, Triboelectrification and razorbacks: geophysical patterns produced in dry grains, Phys. Rev. Let. 96 (2006)178002. [32] Y. Yang, C. Zi, Z. Huang, J. Wang, M. Lungu, Z. Liao, Y. Yang, H. Su, CFD– DEM investigation of particle elutriation with electrostatic effects in gas-solid fluidized beds, Powder Technol. 308 (2017) 422–433.