- Email: [email protected]
Electrostatic field calculation using R-functions and the method of characteristics in electrostatic precipitator
Journal of
ELECTROSTATICS ELSEVIER
Journal of Electrostatics 38 (1996) 269 282
Electrostatic field calculation using R-functions and the method of characteristics in electrostatic precipitator T. B a r b a r i c s a'*, H. Igarashi b, A. Iv~inyi a, T. H o n m a b "Department of Electromagnetic Theo~, Technical Universi O, of Budapest, H-1521, Budapest, Hungat T b Division Systems and Information Engineering, Facul~' c~[Engineering, Hokkaido Univet:~'i(v, Sapporo, Japan Received 15 December 1995: accepted after revision 12 June 1996
Abstract
This paper deals with calculation of the electric field in an electrostatic precipitator using the R-functions and the method of characteristics for the space charge field analysis. The fields are described by the Poisson equation and the equation of current continuity. The current from the wire is controlled by the charge density on the high-voltage electrode, as determined from Peek's empirical law, in order to obtain the current-voltage characteristics of the corona device. Keywords: Electrostatic precipitator; M e t h o d of characteristics: Finite element calculation: Electrostatic field; Wire-duct precipitator: Peek's formula
1. Introduction The prevention of environmental pollution requires techniques for cleaning up coal-boiler flue gases. In general, electrical precipitators are used to control the emission of particles. These devices help to filter the flying ash, avoiding their emission to atmosphere. The electrically charged ash is collected by grounded plates under the action of the electric field. The ash is charged by the ions formed through the corona effect on the high-voltage emitting electrodes. Since the electric field is strongly affected by space charge, and vice versa, in electric precipitators, these quantities must be self-consistently determined by an appropriate numerical technique.
* Corresponding author. Tel.: + (36-1) 463-28-12; fax: + t36-1) 463-31-98: e-mail: barbarics(w,evtsz.bme.hu. 0304-3886/96/$15.00 Copyright ~ 1996 Elsevier Science B.V. All rights reserved. PI1 S 0 3 0 4 - 3 8 8 6 ( 9 6 ) 0 0 0 3 1 -9
270
71 Barbarics et al./Journal ¢#Ele~trostatics 38 (1996) 269 282
The main relations in the field calculation are based upon Maxwell equations. Introducing the scalar or the vector potentials, the solution of the problem can be reduced to solving differential equations. The results have to fulfill the prescribed boundary conditions. The construction of modern electrostatic precipitators needs more and more accurate models to determine the electric field. To solve the problem of the electric field calculation, several numerical methods are used around the world. Presently, the methods of integral equations and procedures based on the variational principle are in extensive use. Many global and discrete versions, such as global variational methods, boundary element methods, finite element methods are applied in practice. Among the numerical methods the most popular ones are the finite element, the boundary element and the finite difference methods. Cristina et al. [1], Davis and Hoburg [2], and Butler et al. [3] used the method of characteristics (MOC) in their work combined with finite elements; Igarashi et al. [4] found a solution with the help of the boundary element method combined with MOC, while finite difference analysis was used by McDonald et al. [5]. Gregory et al. [-6] present a combined finite element and finite difference technique in their paper. The easy formulation and the flexibility are the advantages of the above-mentioned methods. In adopting these methods, the results are determined at grid points. The number of these points in some cases can affect the accuracy of the solution. In case of the finite element method, the value of the potential on the electrode is matched only at the grid points, and not on the whole surface. This effect can influence the accuracy of the electric field calculation and the space charge distribution. For a better solution, more and more unknowns have to be taken into account by the refined mesh. The increased number of the triangular elements can be seen in the finite element mesh of Davis et al. in [2, Fig. 2], Arthur et al. in [3, Fig. 6] and Gregory et al. in [5, Fig. 2]. The finite element grid around the electrode contains many small triangular elements and by this the number of unknowns increases abruptly. The value of the potential on the high-voltage electrode has an influence on the emitted charges, so it has to be known exactly. In an attempt to avoid this problem, the R-functions method is applied. This method ensures that the solution proves the correct value of the potential at every point of the Dirichlet-type boundaries for any selected approximation function. The possibility of modelling the electrodes, not only in a cylindrical shape but in any artificial forms, is the main advantage of the R-functions method. The Ritz method, applied in the global variational method, provides a solution in semi-analytical form at any point in the examined region. Combining the global variational method with R-functions, the solution of the differential equations exactly fulfilling the prescribed boundary conditions can be determined. In the application of the finite element and finite difference methods simulation of the homogeneous Neumann-type boundary conditions can be fulfilled either as forced or as natural conditions. In finite element or finite difference analysis, it is not easy to take into account the forced Neumann-type boundary conditions. To solve this problem, the number of unknowns must be increased by extending the examined region. With the help of the R-functions method, this problem can be considered
7~ Barbarics et al./Journal of Electrostatics 38 (1996) 269 282
271
easily with the construction of the solution including the prescribed Neumann-type boundary conditions as well. Applying this method, the solution can be determined without an increased number of unknowns. A numerical method based on the global variational method, combined with R-functions and the method of characteristics, is applied in this paper to determine the electrical field and the charge density of an electrostatic precipitator.
2. Definition of the problem and the governing equations The electric field and the charge density can be derived from Maxwell's equations by introducing the scalar potential P
V.E =--,
(1)
V. J -- 0,
(2}
J = itpE,
(3)
E =-
(4)
Co
v~0,
where E is the electric field vector, ~o the electric scalar potential, eo the permittivity of the free space, p the ionic charge density, J the current density vector, and I~ the mobility of ions. Eqs. (1) and (4) yield the Laplace-Poisson equation - A q o = - -P.
(5}
Co
The substitution of (1), (2) and (4) into (3) yields p2
(6)
v. Vp + IL-- = 0, C0
where v is the velocity of ions forced byE. Along the lines of electric force, a solution to (6) can be described by the following form: p =
+ It
,
{7)
where pw is the charge density emitted by the high-voltage electrode and t is the time required to move along the characteristic lines from the field point to the surface of the high-voltage electrode [4]. The space charge density p in the region f2 (Figs. 1 and 2) can be obtained from (7), provided that Pw and the electric field are known. The cross-section of the two-dimensional electrostatic precipitator considered in the numerical simulation can be seen in Figs. 1 and 2. Two types of models are examined to show the benefits of the method used here. The Dirichlet-type boundary
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269--282
272
.
.
.
.
.
.
.
.
\ Fig. 1. The cross-sectionof the first model.
\ . . . .
/
\
/
Fig. 2. The cross-sectionof the second model.
conditions are the same for both models: q~lro, = 0,
(8)
rplro~ = qo2,
(9)
where F m is the grounded electrode, and FD2 models the high-voltage electrode with the potential of ~02. For the first model on the boundaries Fyl and /'N2 the homogeneous Neumanntype boundary conditions are fulfilled as natural boundary conditions of the functional, while in the second case, the symmetries are taken into account as forced Neumann-type boundary conditions o n F N 1 and FN2: ~c~n ~
= 0,
(10)
~nn r~2 = 0 .
(11)
F-,i
The models for the field calculation are selected as in [1] for the simulation of the particle transport. To solve the Laplace-Poisson equation (5), the global variational method is applied, where the solution of the differential equation is the potential function, ensuring the energy-related functional to be a minimum.
E Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
273
In a bounded region, the energy-related functional of the static electric field can be formulated with the scalar potential ~0 as W(+p) = ~l f ~ (e.V~oVq~-2p~p)dQ.
(12)
To fulfill the prescribed Dirichlet-type boundary conditions, the scalar potential is decomposed into two parts, the known ~oa and the unknown q~ terms: ~P = ~Pa + ~P~.
(13)
The known component of the potential function is selected to be continuous in the region, fulfilling the prescribed Dirichlet-type boundary conditions along the surface C D = FD1LdFD2 :
q)aIrD, = 0,
(14)
(/961/'D+, = (P2"
{15)
The unknown component fulfills homogeneous Dirichlet-type boundary conditions along these surfaces:
q~lr,,,l~ = 0.
(16)
In the solution concerning the first model, the homogeneous Neumann-type boundary condition is satisfied during the minimisation process as natural boundary condition of the functional. Having selected the known component of the potential function fulfilling the prescribed Dirichlet-type boundary conditions, the unknown component in accordance with the Ritz method can be approximated by the first n elements of an entire function set ¢Pz = E f k a k ,
k = 1,2 . . . . . n,
(17)
where a is a column vector with aa (k = 1,2, ..., n) elements of the unknown coefficients. The kth element of the approximating function Fk can be constructed to fulfil the prescribed homogeneous Dirichlet-type boundary conditions as
Fk(r,z)=WDfk(r,z),
k = 1,2 . . . . . n,
(18)
wherefk(r, z) is the basis functions of the approximation, which was selected to be the Chebyshev polynomial and WDis the weighting function ensuring that at any selection of the basis functions, the prescribed Dirichlet-type boundary conditions are fulfilled for the unknown component of the scalar potential. The elements of a can be determined with the help of the first derivation of (12), which results in a linear set of equations Aa = b,
(19)
where A is an nth-order quadratic matrix and b is a column vector with n elements formed from the known component of the potential function and from the space
274
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269 282
charges. The elements of the matrices A and b can be determined as A [k,1] = ~ c VFk VF~ d ~ , J~2
(20)
b[k] = - ~ (c Vq)~ VFk + pFk)d(2. d~2
(21)
In the second solution, the satisfaction of the homogeneous Neumann-type boundary conditions is achieved by adding an additional term to the potential function (13): q) = ~Pa + q)~ + q~NWDN,
(22)
where WDNis the weighting function having zero value on the surfaces with prescribed Dirichlet- and Neumann-type boundary conditions, obeying the following set of relations: VWDN ]I"D = 0 ,
(23)
VWDNIr, = 1 ,
(24)
d~p ,oN = V(~06+ ~0,). n - d ~ '
(25)
where n is the normal vector of the boundary, which is selected to be n = - VWDN.
(26)
For the above selection, the components of Eqs. (20) and (21) are modified to be (D~ = (f)6 -~- WDN V(D6" gl,
(27)
F~~ = Fk + WDN VFk" n.
(28)
For the modified potential, the result can be determined with the help of (19)-(21), replacing the functions q~6 and Fk with (pp and F~'. 2.1. The R-functions In this case, the field problem is to find the potential function, satisfying the scalar Laplace Poisson equation on the examined region Q, described by the R-functions. The R-functions can be constructed by the reference to geometry. Each prescribed boundary has its own function prescribed by geometry. The functions have to fulfill the following conditions. They must have a zero value at the point P belonging to its boundary, a positive value if the point P belongs to the examined region, and a negative value if it is outside this region [7]. Usually, the boundaries are not simple; they are formed by the union, the intersection, or by the superposition of more subregions. In such cases, the R-function which describes the region can be constructed through R-conjunction, R-disjunction and R-negation.
T. Barbaric~ et al./Journal o f Electrostatics 38 (1996) 269 282
275
Fig. 3. T h e r e g i o n d e s c r i b e d b y R - c o n j u n c t i o n deft} a n d R - d i s j u n c t i o n (right).
R-conjunction is defined on the basis of triangular inequalities. It describes thc union of two subregions resulting in an R-function for the new region: :
2
W I V W 2 = W 1 + W 2 -}- N / W 1
-~- W ~ .
(29)
The realisation of the R-disjunction is W T A W 2 = W 1 -~- W 2 - - ~
-/- tV 2 .
(30)
which describes the intersection of the two subregions described by the R-functions w 1 and w2. The regions described by the R-conjunction and the R-disjunction are shown in Fig. 3. In the solution of the electrostatic problem, the weighting function wD in (18) has to ensure that in the approximation the unknown components of the potenlial function fulfil the homogeneous Dirichlet-type boundary condition along the surface FD of the examined region. In this case there are two Dirichlet-type boundary conditions, so the R-conjunction has to be applied several times, because the surface FDI can be determined with the help of the three functions, using the geometrical parameters a , h, d, e as indicated in Fig. 4: wl = (a + b) - y,
(31)
w.
(32)
= a -
y,
h b w m = ~ x - y + (a + b) + ~ e,
(33)
Wl = (w~ v w.) v w.i.
(34)
The R-function w~ describes the grounded electrode, so it has zero value on the electrode, positive value in the examined region and negative value outside this region.
T. Barbarics et al./Journal o f Electrostatics 38 (1996) 269-282
276
d
e
b ¸
I'D,
C
a
I"N 1
I~o2 r w
B
A
G
Fig. 4. Geometrical parameters.
The R-function of the high-voltage electrode can be formed as w2 = X ~
+ y2 _ rw,
(35)
where rw is the radius of the high-voltage electrode. The weighting function of the Dirichlet-type boundary conditions can be formed with help of the R-disjunction: WD =
W 1 A W2 .
(36)
WD satisfies the conditions of R-functions, its value is equal to zero on the boundaries having Dirichlet conditions. To fulfil the prescribed Dirichlet-type boundary conditions the known component of the potential function can be constructed as ~oa -
-
(p2WI
-
W 1 +
.
(37) W2
By this selection the known component of the potential satisfies the prescribed Dirichlet-type boundary conditions given in (14) and (15), on the high-voltage electrode w2 = 0, and ~oa = ~oa and on the grounded electrode Wl = 0 and q~a = 0. For the second model, the weighting function of the Neumann-type boundary conditions has to be determined as W~V= (Cz + d + e) + x,
(38)
Wv = f z -- x,
(39)
WN = WlV V WV.
(40)
WN has zero value on the Neumann-type boundaries. The potential on the boundaries with Dirichlet-type conditions is fulfilled with the help of the known component of the potential. Thus, the weighting function WDN in (20) for the forced Neumann-type
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
277
boundaries has to have zero value on the Dirichlet-type boundaries, so it can be constructed as WDN
=
WN V Wo .
{41)
3. The solution of the problem Taking the first step to realise the solution, the charge density p over the region Q is set to zero (Po = 0). The electrical field intensity E1 is determined with the aid of Eqs. (4), (5) and (12) for the potential function (13) or (22). To model the effects of the space charges, a virtual ion is moved from a selected point to the surface of the high-voltage electrode with a velocity of v = - l~El and the moving time t is obtained from the following equation: t = ~ ti = 1 / ( H ~ ( f , / I E , I I)),
(42)
w h e r e / i is the length of the ith step on the characteristic line. The substitution of (42) into (7) determines the first approximation for the ion charge density Pl at the point for a selected value of P~I. Having determined the charge density P l at each point, the electric field intensity E2 is calculated again taking into account the new space charge density pl. Two iteration loops were applied for the solution. In the first one, the charge density pj has to give a convergent solution. In the second loop, to find the value of the charge density p and the electric field the charge distribution of the high-voltage electrode must be modified under the condition that the electric field on the surface of the high-voltage electrode Ew is equal to the corona onset field Ec given by the Peek's law in the following form [8]: E~ = 3.1 x 1066(1 + 3.08 x 10-2/(6rw)°5),
(43)
where 5 is the relative air density and rw is the radius of the high-voltage electrode. The quantities are in SI units. The charge density on the electrode is determined by the numerical solution of the non-linear equation [Ew(Pw)[ = Ec. The flow chart of the procedure is plotted in Fig. 5. As can be seen from the flow diagram, the procedure is repeated until a self-consistent solution is achieved.
4. Numerical data The parameters for the examined models rw, a, b, cl, c2, d, e,J'~ , j ) plotted in Fig. 4 are assumed, respectively, to be equal to 2.5, 125, 25, re, 30, 50, 50, ~ , 120 ram. In the first case, a Dirichlet boundary condition was prescribed at point A (on the surface of the high-voltage electrode, ~P2 = 52.5 kV) and on the line C D - E F (~ol = 0). The homogeneous Neumann-type boundary condition between point B and G along the symmetry axis is satisfied by the selection of the basis functions. In the first model, on the lines B-C and F - G , the homogeneous Neumann-type boundary
278
T. Barbarics et al./Journal o ( Electrostatics 38 (1996) 269-282
_
I Field calculation using the R-functions
]
-Aq~ n= p nq/e o
I
I
Tracing along the characteristic lines ealeulation of [ti ] in each nodes
MOC
I
Calculation o f pn in each point o n = [ 1/pw + g / e o Z t i ]'1
Calculation o f
Pwm
E J Pwm ) = Eerit
I
I
Fig. 5. The flow diagram.
conditions are fulfilled as natural boundary condition of the functional. In the second case, the same boundary conditions are applied, but the Neumann-type boundary conditions between the points B and C, and the points F and G are prescribed to be forced ones. The equipotential lines in the cross-section computed by the present method for (]92 52.5 kV can be seen in Figs. 6 and 7. In these cases there was no applied space charge density. The results agree well with other numerical solutions [1]. The electrostatic potential decreases inversely with the distance from a cylindrical electrode, so the equipotential lines around the electrode have to be more compact than along the grounded electrode. The potential difference between adjacent plotted equipotential lines is 5 kV. The difference between the natural and forced Neumann=
279
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
-150 -i00 -50 0
50 I00
----------~
l
I
I
t
0
50
i00
__
150
-I00
-50
Fig. 6. The equipotential lines of the first model.
-150 -I00
-
-50
5o I00
150
-----'---'--"
~ -i00
t -50
l,
l
1
0
50
i00
Fig. 7. The equipotential lines of the second model.
type boundary conditions can be seen in these figures; as in the second case the equipotential lines are perpendicular to the surfaces FN1, FN2. In Fig. 8, the equipotential lines are shown for the second model, and in this case the charge density is taken into account as well. In this model, the radius of the high-voltage electrode seems to be increasing under the effect of the space charges and for this reason, the potential around the high-voltage electrode in Fig. 8 does not
280
T. Barbarics et al./Journal o[ Electrostatics 38 (1996) 269-282
-150 -i00
J
J
-50
D 50
150
i ~ ,
-
-i00
-50
i
I
I
0
50
I00
Fig. 8. The field lines with free space charge.
J I
I
I %
\
% %
% • •
%
I
I
/
I
• / p
%
I
Fig. 9. The lines of electric force.
change so intensively as it does without the space charges. The effect of the free charges can be followed in this figure well. The lines of the electric force can be seen in Fig. 9. An error estimation was made to ensure the validity of the model. In Figs. 10-12 the potential with respect to the distance from the wire can be seen. In Fig. 10 the number of unknowns are determined by a four-point Gauss quadrature, and in Fig. 11 by a six-point Gauss quadrature. In Fig. 12, the effects of changing
1~ Barbarics et al./Journal o]"Eleetrosmtics 38 (1996) 269-282
281
60,00 kV 60,00 kV •~ - - - N g = 6
40,00 kV
A
I
Ng=61
30,00 kV 20,00 kV 10,00 kV 0,00 kV 20
0
40
60
80
100
120
140
160
ITwn Fig. 10. The v a r i a t i o n of the G a u s s q u a d r a t u r e , Nunknow. = 16.
60,00 kV ~Nun=16 I ---'~Nun=25] Nun:36[
I
50,00 kV 40,00 kV 30,00 kV 20,00 kV 10,00 kV 0,00 kV 0
20
40
60
80
100
120
140
160
mnl Fig. 1 I. The v a r i a t i o n of the n u m b e r of the u n k n o w n s , N o .... = 4.
60,00 kV
- - - O - ~ Nun=16 I
---B--- Nun=2Sl
50,00 kV
=-.
Nun=36]
40,00 kV 30,00 kV 20,00 kV 10,00 kV
o,oo kV 0
20
40
60
60 mm
100
120
140
Fig. 12. The v a r i a t i o n of the n u m b e r of the u n k n o w n s , NG ..... = 6.
160
282
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
the n u m b e r s in the G a u s s q u a d r a t u r e can be examined. As the figures show, there is no significant difference between these solutions.
5. Summary T h e n u m e r i c a l analysis of a t w o - d i m e n s i o n a l electrostatic p r e c i p i t a t o r that applies the m e t h o d of characteristics a n d the R-functions m e t h o d has been described. Ado p t i n g R-functions c o m b i n e d with the Ritz m e t h o d for electrical field calculation, the result can be given by the s o l u t i o n of differential equations. T h e result gives the value of the p o t e n t i a l exactly on the p r e s c r i b e d b o u n d a r i e s w i t h o u t having to increase the n u m b e r of u n k n o w n s as in the case of the finite element or the b o u n d a r y element m e t h o d , where the n u m b e r of u n k n o w n s has to be a u g m e n t e d , which i n t r o d u c e s a m o r e c o m p a c t grid. T h e D i r i c h l e t - t y p e b o u n d a r y c o n d i t i o n s are fulfilled a l o n g the b o u n d a r i e s , a n d in the s e c o n d case, the satisfaction of the h o m o g e n e o u s N e u m a n n type b o u n d a r y c o n d i t i o n s are forced by the solution. D u r i n g this field c a l c u l a t i o n only a t w o - d i m e n s i o n a l m e t h o d has been e x a m i n e d , b u t e x t e n d i n g the R-functions m e t h o d , a t h r e e - d i m e n s i o n a l p r o b l e m can be e v a l u a t e d w i t h o u t difficulty.
Acknowledgements T. B a r b a r i c s w o u l d like to express his t h a n k s to the D i v i s i o n Systems a n d I n f o r m a tion Engineering, F a c u l t y of Engineering, H o k k a i d o U n i v e r s i t y for the i n v i t a t i o n to c a r r y o u t further researches in this topic, a n d to the I n t e r n a t i o n a l E x c h a n g e F u n d of H o k k a i d o U n i v e r s i t y for the financial s u p p o r t .
References [1] S. Cristina, G. Dinelli and M. Feliziani, Numerical computation of corona space charge and V - I characteristic in DC electrostatic precipitator, IEEE Trans. Ind. Appl., 27 (1991) 147 153. [2] J.L. Davis and J.F. Hoburg, Wire-duct precipitator field and charge computation using the finite element and characteristics method, J. Electrostatics, 14 (1983) 187-199. [3] A.J. Butler, Z.J. Cendes and J.F. Hoburg, Interfacing the finite-element method with the method of characteristics in self-consistent electrostatic field models, IEEE Trans. Ind. Appl., 25(3) (1989) 533-537. [4] H. Igarashi, S. Sakai, T. Nakamura, T. Morinaga and T. Honma, A boundary element analysis of space charge fields in a corona device, IEEE Trans. Mag., 29 (1993) 1508 1511. [5] J.R. McDonald, W.B. Smith, H.W. Spencer III and L.E. Sparks, A mathematical model for calculating electrical conditions in wire-duct electrostatic precipitator devices, J. Appl. Phys., 48 (1977) 2231-2243. [6] G.A. Kallio and D.E. Stock, Computation of electrical conditions inside wire-duut electrostatic precipitators using a combined finite-element, finite-difference technique, J. Appl. Phys., 59 (1986) 1799 1806. [7] A. Ivanyi, R-functions in electromagnetism, Technical Report No. TUB-TR-93-EE08, Budapest, 1993. [8] G. Hartmann, Theoretical evaluation of Peek's law, IEEE Trans. Ind. Appl., IA-20 (1984) 1647 1651.
ELECTROSTATICS ELSEVIER
Journal of Electrostatics 38 (1996) 269 282
Electrostatic field calculation using R-functions and the method of characteristics in electrostatic precipitator T. B a r b a r i c s a'*, H. Igarashi b, A. Iv~inyi a, T. H o n m a b "Department of Electromagnetic Theo~, Technical Universi O, of Budapest, H-1521, Budapest, Hungat T b Division Systems and Information Engineering, Facul~' c~[Engineering, Hokkaido Univet:~'i(v, Sapporo, Japan Received 15 December 1995: accepted after revision 12 June 1996
Abstract
This paper deals with calculation of the electric field in an electrostatic precipitator using the R-functions and the method of characteristics for the space charge field analysis. The fields are described by the Poisson equation and the equation of current continuity. The current from the wire is controlled by the charge density on the high-voltage electrode, as determined from Peek's empirical law, in order to obtain the current-voltage characteristics of the corona device. Keywords: Electrostatic precipitator; M e t h o d of characteristics: Finite element calculation: Electrostatic field; Wire-duct precipitator: Peek's formula
1. Introduction The prevention of environmental pollution requires techniques for cleaning up coal-boiler flue gases. In general, electrical precipitators are used to control the emission of particles. These devices help to filter the flying ash, avoiding their emission to atmosphere. The electrically charged ash is collected by grounded plates under the action of the electric field. The ash is charged by the ions formed through the corona effect on the high-voltage emitting electrodes. Since the electric field is strongly affected by space charge, and vice versa, in electric precipitators, these quantities must be self-consistently determined by an appropriate numerical technique.
* Corresponding author. Tel.: + (36-1) 463-28-12; fax: + t36-1) 463-31-98: e-mail: barbarics(w,evtsz.bme.hu. 0304-3886/96/$15.00 Copyright ~ 1996 Elsevier Science B.V. All rights reserved. PI1 S 0 3 0 4 - 3 8 8 6 ( 9 6 ) 0 0 0 3 1 -9
270
71 Barbarics et al./Journal ¢#Ele~trostatics 38 (1996) 269 282
The main relations in the field calculation are based upon Maxwell equations. Introducing the scalar or the vector potentials, the solution of the problem can be reduced to solving differential equations. The results have to fulfill the prescribed boundary conditions. The construction of modern electrostatic precipitators needs more and more accurate models to determine the electric field. To solve the problem of the electric field calculation, several numerical methods are used around the world. Presently, the methods of integral equations and procedures based on the variational principle are in extensive use. Many global and discrete versions, such as global variational methods, boundary element methods, finite element methods are applied in practice. Among the numerical methods the most popular ones are the finite element, the boundary element and the finite difference methods. Cristina et al. [1], Davis and Hoburg [2], and Butler et al. [3] used the method of characteristics (MOC) in their work combined with finite elements; Igarashi et al. [4] found a solution with the help of the boundary element method combined with MOC, while finite difference analysis was used by McDonald et al. [5]. Gregory et al. [-6] present a combined finite element and finite difference technique in their paper. The easy formulation and the flexibility are the advantages of the above-mentioned methods. In adopting these methods, the results are determined at grid points. The number of these points in some cases can affect the accuracy of the solution. In case of the finite element method, the value of the potential on the electrode is matched only at the grid points, and not on the whole surface. This effect can influence the accuracy of the electric field calculation and the space charge distribution. For a better solution, more and more unknowns have to be taken into account by the refined mesh. The increased number of the triangular elements can be seen in the finite element mesh of Davis et al. in [2, Fig. 2], Arthur et al. in [3, Fig. 6] and Gregory et al. in [5, Fig. 2]. The finite element grid around the electrode contains many small triangular elements and by this the number of unknowns increases abruptly. The value of the potential on the high-voltage electrode has an influence on the emitted charges, so it has to be known exactly. In an attempt to avoid this problem, the R-functions method is applied. This method ensures that the solution proves the correct value of the potential at every point of the Dirichlet-type boundaries for any selected approximation function. The possibility of modelling the electrodes, not only in a cylindrical shape but in any artificial forms, is the main advantage of the R-functions method. The Ritz method, applied in the global variational method, provides a solution in semi-analytical form at any point in the examined region. Combining the global variational method with R-functions, the solution of the differential equations exactly fulfilling the prescribed boundary conditions can be determined. In the application of the finite element and finite difference methods simulation of the homogeneous Neumann-type boundary conditions can be fulfilled either as forced or as natural conditions. In finite element or finite difference analysis, it is not easy to take into account the forced Neumann-type boundary conditions. To solve this problem, the number of unknowns must be increased by extending the examined region. With the help of the R-functions method, this problem can be considered
7~ Barbarics et al./Journal of Electrostatics 38 (1996) 269 282
271
easily with the construction of the solution including the prescribed Neumann-type boundary conditions as well. Applying this method, the solution can be determined without an increased number of unknowns. A numerical method based on the global variational method, combined with R-functions and the method of characteristics, is applied in this paper to determine the electrical field and the charge density of an electrostatic precipitator.
2. Definition of the problem and the governing equations The electric field and the charge density can be derived from Maxwell's equations by introducing the scalar potential P
V.E =--,
(1)
V. J -- 0,
(2}
J = itpE,
(3)
E =-
(4)
Co
v~0,
where E is the electric field vector, ~o the electric scalar potential, eo the permittivity of the free space, p the ionic charge density, J the current density vector, and I~ the mobility of ions. Eqs. (1) and (4) yield the Laplace-Poisson equation - A q o = - -P.
(5}
Co
The substitution of (1), (2) and (4) into (3) yields p2
(6)
v. Vp + IL-- = 0, C0
where v is the velocity of ions forced byE. Along the lines of electric force, a solution to (6) can be described by the following form: p =
+ It
,
{7)
where pw is the charge density emitted by the high-voltage electrode and t is the time required to move along the characteristic lines from the field point to the surface of the high-voltage electrode [4]. The space charge density p in the region f2 (Figs. 1 and 2) can be obtained from (7), provided that Pw and the electric field are known. The cross-section of the two-dimensional electrostatic precipitator considered in the numerical simulation can be seen in Figs. 1 and 2. Two types of models are examined to show the benefits of the method used here. The Dirichlet-type boundary
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269--282
272
.
.
.
.
.
.
.
.
\ Fig. 1. The cross-sectionof the first model.
\ . . . .
/
\
/
Fig. 2. The cross-sectionof the second model.
conditions are the same for both models: q~lro, = 0,
(8)
rplro~ = qo2,
(9)
where F m is the grounded electrode, and FD2 models the high-voltage electrode with the potential of ~02. For the first model on the boundaries Fyl and /'N2 the homogeneous Neumanntype boundary conditions are fulfilled as natural boundary conditions of the functional, while in the second case, the symmetries are taken into account as forced Neumann-type boundary conditions o n F N 1 and FN2: ~c~n ~
= 0,
(10)
~nn r~2 = 0 .
(11)
F-,i
The models for the field calculation are selected as in [1] for the simulation of the particle transport. To solve the Laplace-Poisson equation (5), the global variational method is applied, where the solution of the differential equation is the potential function, ensuring the energy-related functional to be a minimum.
E Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
273
In a bounded region, the energy-related functional of the static electric field can be formulated with the scalar potential ~0 as W(+p) = ~l f ~ (e.V~oVq~-2p~p)dQ.
(12)
To fulfill the prescribed Dirichlet-type boundary conditions, the scalar potential is decomposed into two parts, the known ~oa and the unknown q~ terms: ~P = ~Pa + ~P~.
(13)
The known component of the potential function is selected to be continuous in the region, fulfilling the prescribed Dirichlet-type boundary conditions along the surface C D = FD1LdFD2 :
q)aIrD, = 0,
(14)
(/961/'D+, = (P2"
{15)
The unknown component fulfills homogeneous Dirichlet-type boundary conditions along these surfaces:
q~lr,,,l~ = 0.
(16)
In the solution concerning the first model, the homogeneous Neumann-type boundary condition is satisfied during the minimisation process as natural boundary condition of the functional. Having selected the known component of the potential function fulfilling the prescribed Dirichlet-type boundary conditions, the unknown component in accordance with the Ritz method can be approximated by the first n elements of an entire function set ¢Pz = E f k a k ,
k = 1,2 . . . . . n,
(17)
where a is a column vector with aa (k = 1,2, ..., n) elements of the unknown coefficients. The kth element of the approximating function Fk can be constructed to fulfil the prescribed homogeneous Dirichlet-type boundary conditions as
Fk(r,z)=WDfk(r,z),
k = 1,2 . . . . . n,
(18)
wherefk(r, z) is the basis functions of the approximation, which was selected to be the Chebyshev polynomial and WDis the weighting function ensuring that at any selection of the basis functions, the prescribed Dirichlet-type boundary conditions are fulfilled for the unknown component of the scalar potential. The elements of a can be determined with the help of the first derivation of (12), which results in a linear set of equations Aa = b,
(19)
where A is an nth-order quadratic matrix and b is a column vector with n elements formed from the known component of the potential function and from the space
274
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269 282
charges. The elements of the matrices A and b can be determined as A [k,1] = ~ c VFk VF~ d ~ , J~2
(20)
b[k] = - ~ (c Vq)~ VFk + pFk)d(2. d~2
(21)
In the second solution, the satisfaction of the homogeneous Neumann-type boundary conditions is achieved by adding an additional term to the potential function (13): q) = ~Pa + q)~ + q~NWDN,
(22)
where WDNis the weighting function having zero value on the surfaces with prescribed Dirichlet- and Neumann-type boundary conditions, obeying the following set of relations: VWDN ]I"D = 0 ,
(23)
VWDNIr, = 1 ,
(24)
d~p ,oN = V(~06+ ~0,). n - d ~ '
(25)
where n is the normal vector of the boundary, which is selected to be n = - VWDN.
(26)
For the above selection, the components of Eqs. (20) and (21) are modified to be (D~ = (f)6 -~- WDN V(D6" gl,
(27)
F~~ = Fk + WDN VFk" n.
(28)
For the modified potential, the result can be determined with the help of (19)-(21), replacing the functions q~6 and Fk with (pp and F~'. 2.1. The R-functions In this case, the field problem is to find the potential function, satisfying the scalar Laplace Poisson equation on the examined region Q, described by the R-functions. The R-functions can be constructed by the reference to geometry. Each prescribed boundary has its own function prescribed by geometry. The functions have to fulfill the following conditions. They must have a zero value at the point P belonging to its boundary, a positive value if the point P belongs to the examined region, and a negative value if it is outside this region [7]. Usually, the boundaries are not simple; they are formed by the union, the intersection, or by the superposition of more subregions. In such cases, the R-function which describes the region can be constructed through R-conjunction, R-disjunction and R-negation.
T. Barbaric~ et al./Journal o f Electrostatics 38 (1996) 269 282
275
Fig. 3. T h e r e g i o n d e s c r i b e d b y R - c o n j u n c t i o n deft} a n d R - d i s j u n c t i o n (right).
R-conjunction is defined on the basis of triangular inequalities. It describes thc union of two subregions resulting in an R-function for the new region: :
2
W I V W 2 = W 1 + W 2 -}- N / W 1
-~- W ~ .
(29)
The realisation of the R-disjunction is W T A W 2 = W 1 -~- W 2 - - ~
-/- tV 2 .
(30)
which describes the intersection of the two subregions described by the R-functions w 1 and w2. The regions described by the R-conjunction and the R-disjunction are shown in Fig. 3. In the solution of the electrostatic problem, the weighting function wD in (18) has to ensure that in the approximation the unknown components of the potenlial function fulfil the homogeneous Dirichlet-type boundary condition along the surface FD of the examined region. In this case there are two Dirichlet-type boundary conditions, so the R-conjunction has to be applied several times, because the surface FDI can be determined with the help of the three functions, using the geometrical parameters a , h, d, e as indicated in Fig. 4: wl = (a + b) - y,
(31)
w.
(32)
= a -
y,
h b w m = ~ x - y + (a + b) + ~ e,
(33)
Wl = (w~ v w.) v w.i.
(34)
The R-function w~ describes the grounded electrode, so it has zero value on the electrode, positive value in the examined region and negative value outside this region.
T. Barbarics et al./Journal o f Electrostatics 38 (1996) 269-282
276
d
e
b ¸
I'D,
C
a
I"N 1
I~o2 r w
B
A
G
Fig. 4. Geometrical parameters.
The R-function of the high-voltage electrode can be formed as w2 = X ~
+ y2 _ rw,
(35)
where rw is the radius of the high-voltage electrode. The weighting function of the Dirichlet-type boundary conditions can be formed with help of the R-disjunction: WD =
W 1 A W2 .
(36)
WD satisfies the conditions of R-functions, its value is equal to zero on the boundaries having Dirichlet conditions. To fulfil the prescribed Dirichlet-type boundary conditions the known component of the potential function can be constructed as ~oa -
-
(p2WI
-
W 1 +
.
(37) W2
By this selection the known component of the potential satisfies the prescribed Dirichlet-type boundary conditions given in (14) and (15), on the high-voltage electrode w2 = 0, and ~oa = ~oa and on the grounded electrode Wl = 0 and q~a = 0. For the second model, the weighting function of the Neumann-type boundary conditions has to be determined as W~V= (Cz + d + e) + x,
(38)
Wv = f z -- x,
(39)
WN = WlV V WV.
(40)
WN has zero value on the Neumann-type boundaries. The potential on the boundaries with Dirichlet-type conditions is fulfilled with the help of the known component of the potential. Thus, the weighting function WDN in (20) for the forced Neumann-type
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
277
boundaries has to have zero value on the Dirichlet-type boundaries, so it can be constructed as WDN
=
WN V Wo .
{41)
3. The solution of the problem Taking the first step to realise the solution, the charge density p over the region Q is set to zero (Po = 0). The electrical field intensity E1 is determined with the aid of Eqs. (4), (5) and (12) for the potential function (13) or (22). To model the effects of the space charges, a virtual ion is moved from a selected point to the surface of the high-voltage electrode with a velocity of v = - l~El and the moving time t is obtained from the following equation: t = ~ ti = 1 / ( H ~ ( f , / I E , I I)),
(42)
w h e r e / i is the length of the ith step on the characteristic line. The substitution of (42) into (7) determines the first approximation for the ion charge density Pl at the point for a selected value of P~I. Having determined the charge density P l at each point, the electric field intensity E2 is calculated again taking into account the new space charge density pl. Two iteration loops were applied for the solution. In the first one, the charge density pj has to give a convergent solution. In the second loop, to find the value of the charge density p and the electric field the charge distribution of the high-voltage electrode must be modified under the condition that the electric field on the surface of the high-voltage electrode Ew is equal to the corona onset field Ec given by the Peek's law in the following form [8]: E~ = 3.1 x 1066(1 + 3.08 x 10-2/(6rw)°5),
(43)
where 5 is the relative air density and rw is the radius of the high-voltage electrode. The quantities are in SI units. The charge density on the electrode is determined by the numerical solution of the non-linear equation [Ew(Pw)[ = Ec. The flow chart of the procedure is plotted in Fig. 5. As can be seen from the flow diagram, the procedure is repeated until a self-consistent solution is achieved.
4. Numerical data The parameters for the examined models rw, a, b, cl, c2, d, e,J'~ , j ) plotted in Fig. 4 are assumed, respectively, to be equal to 2.5, 125, 25, re, 30, 50, 50, ~ , 120 ram. In the first case, a Dirichlet boundary condition was prescribed at point A (on the surface of the high-voltage electrode, ~P2 = 52.5 kV) and on the line C D - E F (~ol = 0). The homogeneous Neumann-type boundary condition between point B and G along the symmetry axis is satisfied by the selection of the basis functions. In the first model, on the lines B-C and F - G , the homogeneous Neumann-type boundary
278
T. Barbarics et al./Journal o ( Electrostatics 38 (1996) 269-282
_
I Field calculation using the R-functions
]
-Aq~ n= p nq/e o
I
I
Tracing along the characteristic lines ealeulation of [ti ] in each nodes
MOC
I
Calculation o f pn in each point o n = [ 1/pw + g / e o Z t i ]'1
Calculation o f
Pwm
E J Pwm ) = Eerit
I
I
Fig. 5. The flow diagram.
conditions are fulfilled as natural boundary condition of the functional. In the second case, the same boundary conditions are applied, but the Neumann-type boundary conditions between the points B and C, and the points F and G are prescribed to be forced ones. The equipotential lines in the cross-section computed by the present method for (]92 52.5 kV can be seen in Figs. 6 and 7. In these cases there was no applied space charge density. The results agree well with other numerical solutions [1]. The electrostatic potential decreases inversely with the distance from a cylindrical electrode, so the equipotential lines around the electrode have to be more compact than along the grounded electrode. The potential difference between adjacent plotted equipotential lines is 5 kV. The difference between the natural and forced Neumann=
279
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
-150 -i00 -50 0
50 I00
----------~
l
I
I
t
0
50
i00
__
150
-I00
-50
Fig. 6. The equipotential lines of the first model.
-150 -I00
-
-50
5o I00
150
-----'---'--"
~ -i00
t -50
l,
l
1
0
50
i00
Fig. 7. The equipotential lines of the second model.
type boundary conditions can be seen in these figures; as in the second case the equipotential lines are perpendicular to the surfaces FN1, FN2. In Fig. 8, the equipotential lines are shown for the second model, and in this case the charge density is taken into account as well. In this model, the radius of the high-voltage electrode seems to be increasing under the effect of the space charges and for this reason, the potential around the high-voltage electrode in Fig. 8 does not
280
T. Barbarics et al./Journal o[ Electrostatics 38 (1996) 269-282
-150 -i00
J
J
-50
D 50
150
i ~ ,
-
-i00
-50
i
I
I
0
50
I00
Fig. 8. The field lines with free space charge.
J I
I
I %
\
% %
% • •
%
I
I
/
I
• / p
%
I
Fig. 9. The lines of electric force.
change so intensively as it does without the space charges. The effect of the free charges can be followed in this figure well. The lines of the electric force can be seen in Fig. 9. An error estimation was made to ensure the validity of the model. In Figs. 10-12 the potential with respect to the distance from the wire can be seen. In Fig. 10 the number of unknowns are determined by a four-point Gauss quadrature, and in Fig. 11 by a six-point Gauss quadrature. In Fig. 12, the effects of changing
1~ Barbarics et al./Journal o]"Eleetrosmtics 38 (1996) 269-282
281
60,00 kV 60,00 kV •~ - - - N g = 6
40,00 kV
A
I
Ng=61
30,00 kV 20,00 kV 10,00 kV 0,00 kV 20
0
40
60
80
100
120
140
160
ITwn Fig. 10. The v a r i a t i o n of the G a u s s q u a d r a t u r e , Nunknow. = 16.
60,00 kV ~Nun=16 I ---'~Nun=25] Nun:36[
I
50,00 kV 40,00 kV 30,00 kV 20,00 kV 10,00 kV 0,00 kV 0
20
40
60
80
100
120
140
160
mnl Fig. 1 I. The v a r i a t i o n of the n u m b e r of the u n k n o w n s , N o .... = 4.
60,00 kV
- - - O - ~ Nun=16 I
---B--- Nun=2Sl
50,00 kV
=-.
Nun=36]
40,00 kV 30,00 kV 20,00 kV 10,00 kV
o,oo kV 0
20
40
60
60 mm
100
120
140
Fig. 12. The v a r i a t i o n of the n u m b e r of the u n k n o w n s , NG ..... = 6.
160
282
T. Barbarics et al./Journal of Electrostatics 38 (1996) 269-282
the n u m b e r s in the G a u s s q u a d r a t u r e can be examined. As the figures show, there is no significant difference between these solutions.
5. Summary T h e n u m e r i c a l analysis of a t w o - d i m e n s i o n a l electrostatic p r e c i p i t a t o r that applies the m e t h o d of characteristics a n d the R-functions m e t h o d has been described. Ado p t i n g R-functions c o m b i n e d with the Ritz m e t h o d for electrical field calculation, the result can be given by the s o l u t i o n of differential equations. T h e result gives the value of the p o t e n t i a l exactly on the p r e s c r i b e d b o u n d a r i e s w i t h o u t having to increase the n u m b e r of u n k n o w n s as in the case of the finite element or the b o u n d a r y element m e t h o d , where the n u m b e r of u n k n o w n s has to be a u g m e n t e d , which i n t r o d u c e s a m o r e c o m p a c t grid. T h e D i r i c h l e t - t y p e b o u n d a r y c o n d i t i o n s are fulfilled a l o n g the b o u n d a r i e s , a n d in the s e c o n d case, the satisfaction of the h o m o g e n e o u s N e u m a n n type b o u n d a r y c o n d i t i o n s are forced by the solution. D u r i n g this field c a l c u l a t i o n only a t w o - d i m e n s i o n a l m e t h o d has been e x a m i n e d , b u t e x t e n d i n g the R-functions m e t h o d , a t h r e e - d i m e n s i o n a l p r o b l e m can be e v a l u a t e d w i t h o u t difficulty.
Acknowledgements T. B a r b a r i c s w o u l d like to express his t h a n k s to the D i v i s i o n Systems a n d I n f o r m a tion Engineering, F a c u l t y of Engineering, H o k k a i d o U n i v e r s i t y for the i n v i t a t i o n to c a r r y o u t further researches in this topic, a n d to the I n t e r n a t i o n a l E x c h a n g e F u n d of H o k k a i d o U n i v e r s i t y for the financial s u p p o r t .
References [1] S. Cristina, G. Dinelli and M. Feliziani, Numerical computation of corona space charge and V - I characteristic in DC electrostatic precipitator, IEEE Trans. Ind. Appl., 27 (1991) 147 153. [2] J.L. Davis and J.F. Hoburg, Wire-duct precipitator field and charge computation using the finite element and characteristics method, J. Electrostatics, 14 (1983) 187-199. [3] A.J. Butler, Z.J. Cendes and J.F. Hoburg, Interfacing the finite-element method with the method of characteristics in self-consistent electrostatic field models, IEEE Trans. Ind. Appl., 25(3) (1989) 533-537. [4] H. Igarashi, S. Sakai, T. Nakamura, T. Morinaga and T. Honma, A boundary element analysis of space charge fields in a corona device, IEEE Trans. Mag., 29 (1993) 1508 1511. [5] J.R. McDonald, W.B. Smith, H.W. Spencer III and L.E. Sparks, A mathematical model for calculating electrical conditions in wire-duct electrostatic precipitator devices, J. Appl. Phys., 48 (1977) 2231-2243. [6] G.A. Kallio and D.E. Stock, Computation of electrical conditions inside wire-duut electrostatic precipitators using a combined finite-element, finite-difference technique, J. Appl. Phys., 59 (1986) 1799 1806. [7] A. Ivanyi, R-functions in electromagnetism, Technical Report No. TUB-TR-93-EE08, Budapest, 1993. [8] G. Hartmann, Theoretical evaluation of Peek's law, IEEE Trans. Ind. Appl., IA-20 (1984) 1647 1651.