Finite volume method for calculation of electrostatic fields in electrostatic precipitators

Journal of Electrostatics 67 (2009) 37–47

Contents lists available at ScienceDirect

Journal of Electrostatics journal homepage: www.elsevier.com/locate/elstat

Finite volume method for calculation of electrostatic fields in electrostatic precipitators N. Neimarlija a, I. Demird zic´ b, S. Muzaferija c, * a

Elektroprivreda BiH dd Sarajevo, Vilsonovo sˇetalisˇte 15, 71000 Sarajevo, Bosnia and Herzegovina Masˇinski falkultet Sarajevo, Vilsonovo sˇetalisˇte 9, 71000 Sarajevo, Bosnia and Herzegovina c ¨ rnberg, Germany CD-adapco, Nordostpark 3-5, 90411 Nu b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 22 February 2007 Received in revised form 14 October 2008 Accepted 14 October 2008 Available online 18 November 2008

This paper presents a numerical method for calculation of coupled electric and space-charge density fields in electrostatic precipitators. It is based on the finite volume discretization of the solution domain by arbitrary polyhedral control volumes and employs an iterative segregated solution procedure of the resulting set of algebraic equations, amounting to a simple, accurate and efficient numerical technique. The method is tested on a number of cases for which analytical solution, numerical and/or experimental results exist. Also, shown are the results of calculation of a 3D model of electrostatic precipitator with spike discharge electrodes. Ó 2008 Elsevier B.V. All rights reserved.

Keywords: Electrostatic precipitator Finite volume method Polyhedral mesh Electrostatic fields

1. Introduction The electrostatic precipitator (ESP) is one of the most common devices used to control fly ash emissions from boilers, incinerators and some other industrial processes. ESP has many advantages compared to mechanical devices such as bag filters or cyclones, because the ESPs can operate over a wide range of gas temperatures and can achieve high particle collection efficiency. Electrostatic precipitation involves several complicated, coupled physical processes: creation of non-uniform electric and space-charge density fields, ionic and electronic charging of particles moving in combined electro- and hydro-dynamic fields, and turbulent transport of charged particles to collection electrodes. Here, we limit ourselves to modeling the steady electric and space-charge density fields and focus on the ability to model real geometries of ESPs with the final aim to extend the methodology to be capable of predicting coupled electro- and hydro-dynamic fields in ESPs. The problem of predicting the distributions of space-charge density and electric field in ESPs has been the subject of much research over the years [1–13]. Leaving out the early classical period of electrostatic field investigation, the use of computers to simulate ESPs processes becomes more and more interesting because of the computer ever increasing capabilities. In the beginning, the finite

* Correspondence to. Tel.: þ49 911 946 433; fax: þ49 911 946 4399. E-mail address: [email protected] (S. Muzaferija). 0304-3886/$ – see front matter Ó 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.elstat.2008.10.007

difference method (FDM) was used. Later, the finite element method (FEM) established itself as a widespread numerical method and found its application in electrostatics. During the past 10–15 years, the finite volume method (FVM), which dominates the field of fluid dynamics, has started to be used for solving electromagnetics problems. Several papers [14–19] deal with the finite volume method using various discretization schemes and algorithms, different from the ones presented in this paper. One of the aims of this paper is to enable numerical simulation of corona discharge and to predict the current–voltage characteristics in 3D models of real, geometrically complex ESPs. Experimental determination of these characteristics is usually expensive in terms of both time and money, and numerical simulations of corona current begin to dominate. Among numerous papers on electrostatic precipitation, only some deal with numerical calculation of electrical characteristics in 3D domains with real geometry of discharge electrodes, e.g.Ref. [20–22]. The experience gained with FVM application to fluid flow [23– 25], solid stress analysis [26–29] and fluid–solid interaction [30,31] is used in this article to develop a method for calculating electric field strength and space-charge density in real ESPs. In Section 2, a mathematical model for corona discharge in an electrostatic field is presented. After that, the finite volume procedure of discretization of the governing equations and the solution algorithm are described. Assessment of the derived model is conducted through a number of cases for which analytical, experimental and/or numerical results exist. Finally, a numerical calculation of a 3D

38

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

model of an ESP is performed and results compared with available experimental data.

2. Mathematical model Several assumptions are made in the model: only the steady state is considered; the effects of charged fly ash particles and the fluctuations of corona are ignored. Also, we ignore space corona discharge and assume zero thickness of the ionization zone.

2.1. Governing equations

rffiffiffi Esc ¼ Ad þ B

d

(6)

a

where Esc is the electric field at the start of corona, A and B are empirical constants, a is the radius of discharge electrode, and d is the relative density of air at standard conditions. Constants A and B for air, and negative corona, usually have values 32.3  105 V/m, and 0.846  105 V/m1/2 [16], respectively, and we assumed d ¼ 1. The Peek’s law (6) was originally developed for cylindrical discharge electrodes, but here a is assumed to be the local radius of curvature of the electrode. Finally, on the symmetry planes, it is assumed that all dependent variables have a zero gradient in the direction normal to the boundary.

Equations which define the electrostatic field in an ESP are Gauss’s law of electrostatics and the current continuity equation:

div E ¼

3. Numerical method

r eo

(1)

div J ¼ divðrKE  DgradrÞ ¼ 0

(2)

where E ¼  grad f is the electric field, r is the ionic space-charge density, eo ¼ 8.854  1012 F/m is the dielectric permittivity of vacuum, J ¼ rKE  Dgrad r is the ionic current density, K is the ion mobility, and D is the coefficient of diffusivity due to molecular and turbulent transport of ions, while f is the electric potential. These equations can be written in the following integral form valid for an arbitrary part of continuum of the volume V bounded by the surface S, with the surface vector s pointing outwards:

Z

Z

eo gradf$ds ¼ 

Z

rdV

(3)

V

S

rKE$ds ¼

S

Z

Dgradr$ds

S

S

jvj $ds ¼

Z

3.1. Solution domain discretization In order to obtain discrete counterparts of Eq. (5) the solution domain is divided into a number of contiguous, non-overlapping control volumes (CVs) which completely fill the domain. Every internal face is shared by two CVs, and the computational points are placed in the CV’s centroids. The boundary nodes, required for the specification of boundary conditions, reside at the centers of boundary cell faces. An example of control volume of arbitrary polyhedral topology, with characteristic geometric quantities, is shown in Fig. 1.

(4) 3.2. Equations discretization

The above equations make a closed set of two equations with two unknowns, f and r. In order to simplify the presentation of the numerical method, they are written in the form of a so-called generic transport equation describing a general convectiondiffusion transport of a conserved quantity j:

Z

In this section the generic equation (5) is discretized by employing the finite volume method described in detail in Ref. [30], and an algorithm for the solution of discretized equations outlined.

Gj gradj$ds þ

Z

After the space discretization, the generic equation (5) is written for each CV:

P1

qj dV

(5)

s1

V

S

Pj

The meaning of the generic variable j and the coefficients vj, Gj, and qj are given in Table 1.

sj

S1 Sj

dj

s2 S2

2.2. Boundary conditions On the collecting electrode, zero electric potential f and zero gradient of r are applied. On the surface of the discharge electrode, the values of both fo and ro are prescribed. However, the value of ro on the discharge electrode is not known in advance and it must be adjusted during calculation. This is done by employing Peek’s law [32], which says that the electric field intensity cannot exceed the following value:

P0 x3

Sn sn i3

Table 1 The meaning of j, vj, Gj, and qj in generic equation (5). Equation Electric potential Space-charge density

j f r

vj 0 KE

Gj e0 D

qj r 0

x1

i1

rPo i2

x2

Fig. 1. A control volume of general polyhedral topology.

Pn

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47 n Z X j¼1

jvj $ds ¼

n Z X j¼1

Sj

|fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl}

Gj gradj$ds þ

Sj

|fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

Convection

Diffusion

Z

qj dV

(7)

V |fflfflfflffl ffl{zfflfflfflfflffl} Source

where n is the number of CV faces. The gradients in Eq. (7) are evaluated by fitting a linear shape function to a set of sampling points consisting of nearest neighbours of point P0, and using the least square method. In order to calculate the surface integrals using the midpoint rule, the variable values at cell-face centers need to be known, for which the following formula is used:

jSO j ¼

 1h  i   1 j þ jPj þ ðgradjÞP0 rj  rP0 þ ðgradjÞPj rj  rPj 2 P0 2 (8)

where rj is the position vector of the cell-face center j. The use the midpoint rule and formula (8) results in a second-order discretization scheme. 3.2.1. Convection Convective flux of variable j through face j of a CV represents the rate at which j is transported in/out the control volume. This term is nonlinear and requires linearization before the process of solving algebraic equations. For cell-face j the convective flux of variable j can be approximated as:

Z

jvj $dszV_ jj j*j

(9)

Sj

where j*j is the value of the dependent variable j in the center of the cell-face j, and V_ jj is the ‘volume flux’ through the surface Sj:

V_ jj ¼

Z

vj $dszvjj $sj

(10)

Sj

In order to suppress the dispersive nature of the second-order expression (8), the cell-face value jj* is calculated blending the FO second-order jSO j and the first-order jj values



SO FO j*j ¼ jFO j þ gj j j  j j



(11)

where gj is the blending factor with a value between 0 and 1, and jFO j is calculated as:

(

jFO ¼ j

jP0 jPj

if if

V_ jj > 0 V_ jj < 0

(12)

the difference between the second-order central difference approximation of the derivative in the direction of vector dj and the value obtained by taking the arithmetic average of the gradients calculated at nodes P0 and Pj. It vanishes if the spatial variation of j is linear or quadratic. Otherwise, its magnitude is proportional to the second-order truncation error of the scheme and as such vanishes with the grid refinement. This correction term detects and smoothes out any unphysical oscillations which could be induced by the numerical procedure. 3.2.3. Source The source term featuring in the generic transport equation (5) is approximated as:

Z

qj dV ¼ qjP0 VP0

Gj gradj$dszGjj ðgradjÞ*j $sj

(15)

V

3.2.4. Boundary conditions Boundary conditions have to be applied on the faces coinciding with the solution domain boundaries. In case of Dirichlet boundary condition, the expressions for diffusion fluxes and sources remain valid, except for replacing jPj by the boundary value jB. On the boundary regions where Neumann boundary conditions are specified, the boundary fluxes are added to the source term, while the variable values at the boundary are obtained using a second-order extrapolation. The algorithm for adjusting the ionic space-charge density ro on the discharge electrode such that the electric field intensity Eo satisfies the Peek’s law (6) is described in detail in Ref. [16] and is used in the present work. 3.2.5. Resulting set of algebraic equations Assembling all the terms featured in Eq. (7) results in one algebraic equation for each CV and each dependent variable. This equation links the value of the variable at the CV center with the values at the centers of the neighboring CVs:

aj0 jP0 

ni X

ajj jPj ¼ bj ;

(16)

j¼1

where ni is the number of internal cell faces surrounding cell P0, and the right hand side bj contains source terms, contributions from boundary faces, and parts of convective and diffusive fluxes which are, for the sake of computational efficiency, treated explicitly:

  s $s ajj ¼ Gjj dj $sj  min V_ jj ; 0 ; Pn j j P aj0 ¼ ajj þ nj¼ 1 V_ jj ; j¼1 " # n X s $s Gjj ðgradjÞj $sj  gradj$dj j j bj ¼ dj $sj j¼1

3.2.2. Diffusion Diffusion flux of variable j through an internal cell-face j can be approximated as:

Z

39



n g X j j¼1



(13)

2

n g X j j¼1

2

V_ jj V_ jj

h  i   rj  rP0 $ðgradjÞP0 þ rj  rPj $ðgradjÞPj h





jPj  jP0 sgn V_ jj

i

þ qjP VP0 þ 0

nB X

ajB jB

j¼1

(17)

Sj

where Gjj stands for the value of diffusion coefficient in the cellface center, and the gradient ðgradjÞ*j is constructed in the following manner:

"

ðgradjÞ*j

¼

ðgradjÞSO j þ

jPj  jP0   dj 

#  gradj$dj dj sj    dj  dj $sj

where nB ¼ n  ni is the number of the boundary faces surrounding cell P0. Note that according to Table 1 and Gauss´s law of electrostatics:

(14)

The second term on the right hand side in [] brackets represents

n X j¼1

V_ jj ¼



0 K

eo rP0 VP0

for electric potential; j ¼ f for space-charge density; j ¼ r

(18)

40

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

3.3. Solution algorithm

4.1. Cylindrical electrostatic precipitator

By integrating Eq. (7) over all N CVs, two mutually coupled sets of N nonlinear, algebraic equations of the form (16) are obtained. Due to the nonlinearity of the underlying equations, the solution of this set of algebraic equations is obtained by employing a segregated iterative algorithm, i.e. these equations are linearized, and sets of equations for each dependent variable temporarily decoupled by assuming that coefficients and sources terms are known (calculated by using dependent variable values from the previous iteration). As a result, a set of linear algebraic equations for each dependent variable is obtained:

A cylindrical ESP consists of a wire as the discharge electrode and an earthed coaxial cylinder (pipe) as the collecting electrode, as shown in Fig. 2. The geometry of ESP is defined by the diameters of the wire and pipe, 2a ¼ 2.77 mm and 2b ¼ 203.2 mm, respectively. For a uniform distribution of electric quantities along the wire, the problem is one-dimensional. However, in order to demonstrate the method’s unique ability to use cells of arbitrary polyhedral topology, a two-dimensional solution domain and the numerical mesh shown in Fig. 2 are used. It is assumed that there is no diffusion transport D ¼ 0 in the ESP.

Aj j ¼ bj

4.1.1. Constant space charge First we consider the case when the corona current is absent but the presence of constant electric space charge r is assumed (implying the ion mobility K ¼ 0). In this case Eq. (2) reduces to an identity and the analytical solution of the Eq. (1) has the form:

(19)

where Aj is an N  N matrix, vector j contains values of dependent variable j at N nodal points and bj is the source vector. The matrices Aj are sparse, with the number of nonzero elements in each row equal to the number of nearest neighbours plus one, ni þ 1, and diagonally dominant, which makes the equation system (19) easily solvable by a number of iterative method which retain the sparsity of the matrix Aj. Here, the conjugate gradient method with an incomplete Cholesky preconditioning is used [33]. There is no need to solve Eq. (19) to a tight tolerance, since the coefficients and sources are only approximate (based on the values of dependent variables from the previous iteration) and reduction of the sum of absolute residuals by an order of magnitude normally suffices. The solution procedure can be described by the following sequence: 1. Provide an initial guess j0 (usually zero values) for all dependent variables. 2. Increase the iteration counter k ¼ k þ 1 3. Assemble and solve Eq. (19) for the electric potential fk. 4. Assemble and solve Eq. (19) for the space-charge density rk. 5. Adjust boundary conditions for rko on the discharge electrode. 6. Return to step 2 and repeat steps 2–5 until a converged solution is obtained, i.e. until the residual norm





krkj k ¼ Akj jk  bkj

1  bj

bj

Dkj

and





  lnr  lnb ro  2 r  b  r 2 þ fo  o b2  a2 4eo 4eo lna  lnb

(22)

where r is radial distance from the wire center. For r ¼ ro ¼ 20 mC/ m3 and applied voltage of fo ¼ 50 kV, the analytical and computational results are shown in Fig. 3. Comparison between the two solutions shows a very good agreement. In order to demonstrate the second-order convergence property of the numerical method, the numerical mesh is systematically refined and the dependence of average numerical error on the characteristic size of numerical mesh is shown in Fig. 3. The average error is found by summing the weighted absolute values of differences between the predicted electric potentials and those given by Eq. (22). The weighting factors are calculated as the ratio of the CV volume and the total volume of the computational domain. The error is normalized with the average value of the corresponding electric potential in the computational domain. A characteristic grid size is obtained by dividing the square root of the total area of the computational domain by the total number of control volumes. It is obvious that the convergence towards the analytical solution is asymptotically of the second order.

(20)

is reduced by a prescribed number (typically three to four) of orders of magnitude. Note that in order to promote stability of the solution method, an under-relaxation is often necessary. This has been introduced by replacing Akj and bkj with:

Akj þ

fðrÞ ¼

bkj þ

1  bj

bj

Dkj jk1

(21)

4.1.2. Stationary corona current flow In this case, it is assumed that corona current occurs and simultaneous solution of Eqs. (3) and (4) is required. The analytical solution of the corresponding Eqs. (1) and (2) for the electric potential is after White [34]:

0 1 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Jo0 b2 2 BC þ 2peo K þ C C Jo0 r 2 Jo0 b2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiC fðrÞ ¼ þ C2  þ C 2 þ ClnB @ A 2peo K 2peo K J0o r 2 þ C2 Cþ 2peo K

r þ Cln b

respectively, where bj is an under-relaxation factor with a value between 0 and 1, Dkj is a diagonal matrix consisting from the diagonal elements of matrix Akj and jk1 is the dependent variable vector from the previous iteration. In addition, this enhances the diagonal dominance of the linearized equations, which improves the rate of convergence of most iterative linear equations solvers.

where C is

4. Test cases

while the distribution of the electric space-charge density can be expressed as:

The method is verified on a simple cylindrical type of ESP and on the more complex problem of a wire-plate electrostatic precipitator through comparison with available analytical, numerical or experimental results.

(23) sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Jo0 2  C ¼ a Esc 2peo K

rðrÞ ¼ 

J00 1 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   2prK 0 Jo a2 þ a2 E2 1  2 pe 2 oK r r2 sc

(24)

(25)

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

41

Smmetry

y Collecting pipe φ=0 grad ρ = 0 Solution domain x

2b

2a

Simmetry

wire φ = φo ρ = ρo

Fig. 2. Model geometry, solution domain and boundary conditions (left) and example of numerical mesh (right).

100

Electric potential, kV

50 Computational Analytical

40

10

30 20 1

10 0

20

40

60

80

100

0.001

Coordinate r, mm

0.01

0.1

Characteristic dimension, m

Fig. 3. Electric potential in a constant space-charge cylindrical ESP (left) and average numerical error (right). Dashed line represents theoretical slope for a second-order method.

where J0 o is obtained from Eq. (23) from the boundary condition fo ¼ f(a), which enables calculation of ro ¼ r(a) from Eq. (25). In this case, a numerical mesh consisting of 72  16 ¼ 1152 CVs with an expansion factor of 1.25 is used. The calculations are

performed for K ¼ 2.2  104 m2/Vs, and fo ¼ 50 kV (J0 o ¼ 1.2 mA/ m). Results of the calculations of the electric potential and spacecharge density are shown in Fig. 4, where a good agreement between numerical and White’s analytical results can be seen.

120

50 Computational Analytical

Density of electric charge, μC/m3

Electric potential, kV

40

30

20

10

0

Computational Analytical

100

80

60

40

20

0

20

40

60

Coordinate r, mm

80

100

0

0

20

40

60

80

Coordinate r, mm

Fig. 4. Electric potential (left) and space-charge density (right) in a cylindrical ESP with corona current flow.

100

42

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

c

Symmetry

A Symmetry O

y

x Symmetry E

Wire φ = φo ρ = ρo B

Symmetry

2a

2b

Solution domain

C Collecting plate φ = 0, gradρ = 0

D

Fig. 5. Geometry, solution domain, and boundary conditions of a wire-plate ESP (left) and details of the numerical mesh near the discharge electrode (right).

Fig. 6. Results of calculation of electric potential f (left) and electric field intensity jEj (right) in wire-plate ESP without space electric charge.

50 Computational Analytical

40

Electric potential, kV

Electric potential, kV

50

30 20 10 0

0

20

40

60

80

100

Distance from discharge electrode (x=0), mm

Computational Analytical

40

30

20

10

0

20

40

60

80

Distance from discharge electrode (y=0), mm

Fig. 7. Variation of electric potential in wire-plate ESP along line OC (left) and along line OE (right).

4.2. Wire-plate electrostatic precipitator 30 FEM Kallio and Stock FDM Kallio and Stock Exp. Penney and Matick FVM (2268 CV)

25

20

c = 495 mm

15

10

Spike discharge electrode Solution domain

5

y 0

0

20

40

60

80

100

120

Distance from discharge electrode, mm Fig. 8. Electric potential in wire-plate ESP with corona discharge at electrode.

O

x

Collecting electrode (plate) Fig. 9. Scheme of ESP with spikes.

2b = 400 mm

Electric potential, kV

The row of wires (discharge electrodes) positioned midway between two plates (collecting electrodes) is a typical configuration of a wire-plate ESP, as shown in Fig. 5, together with the solution domain and boundary conditions.

43

85

15

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

Pike

33.7

Spike

Spike (φ 3.8)

Pike 100

Fig. 10. Discharge electrode shape and its dimensions (in millimeters).

Fig. 11. Solution domain with numerical mesh (left) and mesh detail near the spike tip (right).

In the examples that follow, the analytical results are obtained by neglecting the ion diffusivity (D ¼ 0).

Current density x 10−4, A/m2

4 Unit 5 - 4. field, measured data 3D num. model, Grid 1 3D num. model, Grid 2 3D num. model, Grid 3

3

4.2.1. Case with no space electric charge First, the problem with no corona discharge, i.e. when r ¼ ro ¼ 0, is considered, for which an analytical distribution of the electric potential (for an infinitely thin wire) can be obtained [34]:

h i   p ðx  mcÞ  cos py cos h 2b 2b h i   fðx; yÞ ¼ l ln py p ðx  mcÞ þ cos 2b m ¼ N cos h 2b þN X

2

(26)

where l is the parameter which is determined from the condition: 1

    mc  cos pa cos h p2b 2b     fo ð0; 0Þ ¼ l ln p p mc þ cos 2ba m ¼ N cos h 2b þN X

0 20

30

40

50

60

Electric potential, kV Fig. 12. Comparison V–I characteristics 3D models with measured data [39].

70

(27)

The ESP considered has wires of diameter 2a ¼ 2.77 mm, and distances between discharge electrodes of c ¼ 152.4 mm and between collecting electrodes 2b ¼ 203.2 mm (Fig. 5).

44

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

El.Potential[V]

El.Field[V/m]

4.999e+04 4.515e+04 4.030e+04 3.546e+04 3.062e+04 2.577e+04 2.093e+04 1.608e+04 1.124e+04 6.393e+03 1.548e+03

2.000e+06 1.801e+06 1.603e+06 1.404e+06 1.205e+06 1.007e+06 8.079e+05 6.092e+05 4.105e+05 2.118e+05 1.310e+04

Fig. 13. Electrical potential (left) and intensity of electric field (right, jEjmax ¼ 7.214  106 V/m).

El.Charge[C/m**3]

El.Current[A/m2]

3.000e−04 2.700e−04 2.400e−04 2.100e−04 1.800e−04 1.500e−04 1.200e−04 9.000e−05 6.000e−05 3.000e−05 4.759e−27 y

z

x

5.000e−03 4.500e−03 4.000e−03 3.500e−03 3.000e−03 2.500e−03 2.000e−03 1.500e−03 1.000e−03 5.000e−04 0.000e+00 y

z

x

Fig. 14. Space-charge Density (left, rmax ¼ 1.736  1002 C/m3) and intensity of electrical current density (right, jJjmax ¼ 2.004  1001 A/m2.)

The solution domain is subdivided into 3558 CVs by the numerical mesh whose detail is shown in Fig. 5 (right). Calculations are performed for an applied voltage of fo ¼ 50 kV on the wire surface. Fig. 6 shows contours of electric potential f and electric field intensity jEj obtained by numerical calculation. The results of numerical calculations of electric potential together with the analytical values are shown in Fig. 7. Obviously, numerical and analytical results agree very well. 4.2.2. Case with corona discharge at electrode An analytical solution for this problem does not exist, hence our numerical solution of the coupled equations (3) and (4) is compared with the numerical results of Kallio and Stock [35], who used FEM and FDM, and with the experimental results of Penney and Matick [36]. The geometry of the model used is defined by: a ¼ 0.152 mm, b ¼ 114.3 mm, c ¼ 152.4 mm. The ion diffusivity is set to zero (D ¼ 0), while the ion mobility K ¼ 2  104 m2/Vs. Calculations are performed for the given value of electric current density at the collecting plate Jo ¼ 3.77  104 A/m2 and the electric potential on the wire surface fo is adjusted iteratively by an algorithm similar to one for adjusting ro, as described in Section 3.2.4. The calculation domain is divided into 2268 CVs. Numerical results obtained by the present FVM and FEM and FDM results of Kallio and Stock [35] are shown in Fig. 8. It can be seen that all three numerical methods give nearly the same results which are in good agreement with experimental results of Penney and Matick [36]. Of note is the fact that the FDM predicts a lower

electric potential value at the discharge electrode surface (below 22 kV), while the FEM and FVM methods yield approximately equal values of 25 kV. More details and discussions can be found in Ref. [35] and [37]. Calculations presented in this section show that the implementation of the present method is correct, that discretized equations are consistent with the governing equations of the mathematical model, and that the discrete solution converges towards the exact one in a manner consistent with a second-order scheme.

5. Application to ESP with spike discharge electrodes In this section the three-dimensional corona discharge problem of an ESP with realistic discharge electrode geometry is considered. A parallel plate-wire design with rigid discharge electrodes to which many sharpened spikes are attached has evolved as one of the best discharge electrode designs. Fig. 9 shows a series of wires in a standard ESP configuration. The vertical arrangement and the geometry of spikes are shown in Fig. 10. The hatched area in Fig. 9 of height (normal to the x–y plane) 200 mm is identified as the solution domain. All boundary surfaces except surfaces of the electrodes are assumed to be symmetry planes. The ion mobility was K ¼ 1.6  104 m2/Vs and the diffusivity coefficient D ¼ 2  102 m2/s. The space-charge density ro on the discharge electrode is adjusted by the iterative procedure described in [16].

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

45

El.Potential[V] 4.976e+04 4.495e+04 4.014e+04 3.533e+04 3.052e+04 2.571e+04 2.090e+04 1.609e+04 1.128e+04 6.468e+03 1.658e+03

y

x

z

Fig. 15. Iso-surfaces of electric charge (r ¼ 5, 10, 20, 30, 40, 50 mC/m3) colored by values of electric potential

gradually assume prismatic form with an elliptic base. The main feature of the electric field in this region is that it is not uniform, i.e. there exists a component of vector E parallel to the plate which moves the electric charge r to the border of the ellipse. In order to verify numerical results, the computed V–I characteristics are compared with those measured at the fourth (outlet) field of the Kakanj Power Plant (B&H) [39] (Fig. 12) because the process conditions in this field are approximately equal to those in a particle-free ESP, as the one modeled here. One can see a qualitatively good agreement, especially in the range 20–55 kV of applied potentials. Beyond this range, above 55 kV, the measured V–I characteristics tend to slope upwards as is often characteristic for measured V–I characteristics in flue gas conditions, especially when back corona occurs. Some authors in a similar analysis attributed this discrepancy to the likely ‘spotty’ nature of negative corona on large wires or to the electron current [22,35]. A better quantitative agreement was not expected because the real ESP uses so-called ‘sigma’ collecting electrodes (Fig. 16) whose surface is somewhat bigger than the flat plate in the numerical example, resulting in smaller values of current flux. Also, the presence of even a small concentration of particles in this field can have some influence and decrease the current flux. Fig. 17 shows a good correlation between the computed electric current (flux) intensity on the collecting electrode surface and the photo showing the shape of marks left after removing the collected ash particles. One can see that both calculated and measured contours have elliptic form, with similar ratios of the axes

30 50

According to Peek’s law [32], for the range of applied voltages on the discharge electrodes of fo ¼ 50, 60 and 70 kV, the spike surface (radius 33.7/2 mm; Eo ¼ 38.8 kV/cm, fo ¼ 132.3 kV) does not initiate corona, but the cylindrical part of the spike discharge electrode (radius 3.8/2 mm; Eo ¼ 51.7 kV/cm, fo ¼ 41.3 kV) must be considered as a potential corona emitter. Three systematically refined meshes with 2119 CV (Grid1), 17614 CV (Grid2) and 144460 CV (Grid3) are generated. Grid2, together with a detail showing the locally refined mesh around the spike tip is shown in Fig. 11. In order promote stability of the solution process, it was necessary to use under-relaxation factors (see Eq. (21)) bf ¼ 0.6 and br ¼ 0.83. Fig. 12 shows the V–I characteristics obtained for three numerical meshes considered. It can be seen that the results for Grid2 are close to those of Grid3, indicating that they are close to a gridindependent solution, and the numerical results that follow are obtained on Grid2. Fig. 13 shows the distribution of electric potential and electric field intensity for fo ¼ 50. As can be seen, the electric field intensity is the highest in the vicinity of the spike tip, where the most intensive discharge processes are expected, as shown in Fig. 14. Indeed, if the electric potential on the discharge electrode is successively increased, it is obvious that the discharge processes will start first on the spike tip, so it is logical that the biggest spacecharge density appears on the spike tip. Very high intensity of electric field at the spike tip (which is possible to obtain by numerical calculation without corona process) is suppressed by discharge processes near the surface of the coronating electrode, i.e. the electric field at the electrode surface is prevented from exceeding its critical value by the production of electrons and formation of space charge. As a consequence, a very high density of space electric charge and electric current appears in this zone, thus decreasing the electric field of the spike tip below its local critical value, making the electric charge drop to zero. For a better illustration of the electrostatic fields in 3D space, Fig. 15 shows iso-surfaces for various values of space charge (5, 10, 20, 30, 40 and 50 mC/m3). Further from the spike, the contours

45

35 480

Fig. 16. Sigma collecting electrode of the Kakanj Power Plant (dimensions in mm).

46

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47

Fig. 17. Contours of ionic current density (left, max. ¼ 3.79  04, min. ¼ 2.27  010 A/m2) and particles remained after washing, ESP of TPP Kakanj [39] (right).

Fig. 18. Dust layers on the collecting surfaces [38] (left) and [40] (right).

(computed 1.41, measured 1.44) and approximately the same positions of the centers of the ellipses. Fig. 18 shows the experimental results of Barbarics [38] (who also performed numerical simulations) and of Blanchard et al. [40] who analyzed density and distribution of the dust layers formed on the collecting surfaces. One can see that the computed and measured contours have elliptic forms similar to those in Fig. 17. This is in accord with theoretical and experimental analysis in [38,40] where it is shown that there is a strong correlation between the particle layers and the ionic current density distributions. Of course, this is only an indirect comparison because no measured data for the electric quantities in the real ESP have been found.

6. Conclusions A simple and efficient numerical method is developed for the analysis of electrostatic fields in ESPs. The method is based on the finite volume discretization which uses unstructured meshes made of arbitrary polyhedral control volumes enabling an easy discretization of complex geometries encountered in the ESPs. It employs a segregated, iterative solution algorithm resulting in low memory requirements. The method is amenable to fully scalable

parallel applications, thus allowing very large models to be run efficiently. The test cases presented, chosen for the availability of either analytical solutions or other numerical results, demonstrate a very good accuracy for our method. Calculations of electric potential and electric space-charge density fields in an ESP with spike discharge electrodes show good qualitative agreement with measurements on real and experimental ESPs. To enable a more rigorous comparison, it is necessary to include in the model a two-phase gas-particle flow, as described in Ref. [37]. This change will be the subject of a future publication. References [1] P. Cooperman, A new theory of precipitator efficiency, Atmospheric Environment 5 (1971) 541–551. [2] D.O. Heinrich, Der große gassenabstand im elektrofilterbau, Staub – Reinhalt. Luft 38 (1978) 446–451. [3] G. Cooperman, A unified efficiency theory for electrostatic precipitator, Atmospheric Environment 18 (1984) 277–285. [4] G.A. Kallio, D.E. Stock, Flow visualization inside a wire-plate electrostatic precipitator, IEEE Transactions on Industry Applications 26 (1990) 503–514. [5] B.S. Choi, C.A.J. Fletcher, Computation of particle transport in an electrostatic precipitator, Journal of Electrostatics 40 (1997) 413–418. [6] J.R. McDonald, W.B. Smith, W. Spencer, L.E. Sparks, A mathematical model to predict performance of a plate-wire electrostatic precipitation devices, Journal of Applied Physics 48 (1997) 2231–2243.

N. Neimarlija et al. / Journal of Electrostatics 67 (2009) 37–47 [7] S.H. Kim, K.W. Lee, Experimental study of electrostatic precipitator performance and comparison with existing theoretical prediction models, Journal of Electrostatics 48 (1999) 3–25. [8] P.L. Levin, J.F. Hoburg, Donor cell-finite element descriptions of wire-duct precipitator fields, IEEE Transactions on Industry Applications 26 (1990) 662–670. [9] J.H. Davidson, P.J. McKinny, P. Linnebur, Three-dimensional (3-D) model of electric field and space charge in the barbed plate-to-plate precipitator, IEEE Transactions on Industry Applications 32 (1996) 858–866. [10] W. Egli, U. Kogelschatz, E.A. Gerteisen, R. Gruber, 3D computation of corona, ion induced secondary flows and particle motion in technical ESP configurations, Journal of Electrostatics 40–41 (1997) 425–430. [11] S.H. Kim, H.S. Park, K.W. Lee, Theoretical model of electrostatic precipitator performance for collecting polydisperse particles, Journal of Electrostatics 50 (2001) 177–190. [12] I. Beux, A. Iollo, M.V. Salvetti, A. Soldati, Approximation and reconstruction of the electrostatic field in wire-plate precipitators by a low-order model, Journal of Computational Physics 170 (2001) 893–916. [13] G.R. Offen, R.F. Altman, Issues and trends electrostatic precipitator for U.S. utilities, Journal of Air and Waste Management Association 41 (1991) 222–227. [14] X. Li, I.R. Ciric, M.R. Raghuveer, Highly stable finite volume based relaxation iterative algorithm for solution of DC line ionized fields in the presence of wind, International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 10 (1997) 355–370. [15] H. Lei, L.-Z. Wang, Z.-N. Wu, Applications of upwind and downwind schemes for calculating electrical conditions in a wire-plate electrostatic precipitator, Journal of Computational Physics 193 (2004) 697–707. [16] A.J. Medlin, C.A. J Fletcher, R. Morrow, A pseudotransient approach to steady state solution of electric field-space charge coupled problems, Journal of Electrostatics 43 (1997) 39–60. [17] A.M. Meroth, A.K. Rastogi, A.J. Schwab, Numerical computation of the turbulent particulated flow in an electrostatic precipitator. In: Int. Symp. Filtration and Separation of Fine Dust (1996) 994–1001. [18] A.M. Meroth, T. Gerber, C.D. Munz, P.L. Levin, A.J. Schwab, Numerical solution of nonstationary charge coupled problems, Journal of Electrostatics 45 (1999) 177–198. [19] C.U. Bo¨ttner, The role of the space charge density in particulate processes in the example of the electrostatic precipitator, Powder Technology 135–136 (2003) 285–294. [20] J.A. Houlgreave, K.S. Bromley, J.C. Fothergill, A finite element method for modelling 3D field and current distributions in electrostatic precipitators with electrodes of any shape. In: 6th International Conference on Electrostatic Precipitation (VI ICESP), Budapest (1996). [21] J.A. Houlgreave, K.S. Bromley, J.C. Fothergill, An electroquasistatic finite element model of electric fields and currents generated by coronating electrodes in an electrostatic precipitator, Computation in Electromagnetics 10–12 (April) (1996) 121–126. [22] B.S. Rajanikanth, N. Thirumaran, Modeling the pre-breakdown V–I characteristics of an electrostatic precipitator, Fuel Processing Technology 76 (2002) 159–186.

47

[23] I. Demird zic´, A.D. Gosman, R.I. Issa, A finite volume method for prediction of turbulent flow in arbitrary geometries, Lecture Notes in Physics 141 (1980) 144–150. [24] I. Demird zic´, M. Peric´, Finite volume method for prediction of fluid flow in arbitrary shaped domains with moving boundary, International Journal of Numerical Methods in Fluids 10 (1990) 771–790. [25] S. Muzaferija, A.D. Gosman, Finite volume CFD procedure and adaptive error control strategy for grids of arbitrary topology, Journal of Computational Physics 138 (1997) 766–787. [26] I. Demird zic´, S. Muzaferija, Finite volume method for stress analysis in complex domains, International Journal of Numerical Methods in Engineering 37 (1994) 3751–3766. [27] H. Basˇic´, I. Demird zic´, S. Muzaferija, Finite volume method for simulation of extrusion processes, International Journal of Numerical Methods in Engineering 62 (2005) 475–494. [28] I. Demird zic´, E. D zaferovic´, A. Ivankovic´, Finite-volume approach to thermoviscoelasticity, Numerical Heat Transfer, Part B 47 (2005) 213–237. [29] I. Bijelonja, I. Demird zic´, S. Muzaferija, A finite volume method for incompressible linear elasticity, Computer Methods in Applied Mechanics and Engineering 195 (2006) 6378–6390. [30] I. Demird zic´, S. Muzaferija, Numerical method for coupled fluid flow, heat transfer and stress analysis using unstructured moving meshes with cells of arbitrary topology, Computer Methods in Applied Mechanics and Engineering 125 (1995) 235–255. [31] C.J. Greenshields, H.G. Weller, A. Ivankovic´, The finite volume method for coupled fluid flow and stress analysis, Computer Modeling and Simulation in Engineering 4 (1999) 213–218. [32] F.W. Peek, Dielectric Phenomena in High Voltage Engineering, McGraw-Hill, New York, 1929. [33] J.A. Meijerink, H.A. van der Vorst, Guidelines for the usage of incomplete decompositions in solving sets of linear equations as they occur in practical problems, Journal of Computational Physics 44 (1981) 134–155. [34] H.J. White, Industrial Electrostatic Precipitation, Addison-Wesley, Reading, Massachusetts, 1963. [35] G.A. Kallio, D.E. Stock, Computation of electrical conditions inside wire-duct electrostatic precipitators using a combined finite-element, finite difference technique, Journal of Applied Physics 59 (1986) 1799–1806. [36] G.W. Penney, R.E. Matick, Potentials in dc corrona fields, Transaction AIEE 79 (1960) 91–99. [37] N. Neimarlija, Numerical modeling of separation of ash particles in electrostatic precipitators, PhD Thesis, University of Sarajevo, Bosnian (2006). [38] T. Barbarics, A. Ivanyi, Determination of the charged particle’s movement in electrostatic precipitator, Journal of Electrical Engineering 48 (1997) 39–42. [39] Research Cottrell, Project: Kakanj Unit 5, Archives of Thermo-Power Plant Kakanj, Deutschland GmbH, Kakanj, Bosnia and Herzegovina, 1998. [40] D. Blanchard, P. Atten, L.M. Dumitran, Correlation between current density and layer structure for fine particle deposition in a laboratory electrostatic precipitator, IEEE Transactions on Industry Applications 38 (2002) 832.