DC Ion Flow Fields : Uniqueness of Solution and Application of the Charge Simulation Method byJ.N.SHENG,Z.YANU?ZdB. Xi&
Jiaotong
University,
L.QIN
The People’s
Republic
of China
and G.GELA Department
of Electrical
Engineering,
Ohio State
University,
Columbus,
Ohio 43210, U.S.A.
ABSTRACT : Some theoretical and practical (computational) questions are addressed related to the solution of DC ion $owjeld problems, such as those of overhead HVDC transmission lines, electrostatic precipitators and DC power equipment with gaseous insulation. Ions are produced by corona activities at energized parts. These ions lead to formation of space charge whichjlls the surrounding space and infIuences the electricjeld distribution. The arrangement of the space charge itself is, at the same time, governed by the electric field in space. The problem is therefore nonlinear, and iterative approaches are normally employed to compute the field quantities. The question of uniqueness of solution is addressed$rst. Then, conditions neededfor convergence of iterative algorithms are studied. A new and very successful algorithm is presented. An important feature of the proposed algorithm is the use of the Charge Simulation Methodfor solution of Poisson fields. Numerical examples are given to illustrate the presented concepts.
I. Introduction
A difficult yet fascinating problem in the subject of electrostatics is furnished by the ion flow field associated with DC corona. Unlike with AC excitation where reversal of field direction takes place every half cycle, ions produced by DC corona activities are set by the local DC electric field, into steady motion away from the conductor in the corona. The space around the conductors becomes filled with these ions and the resultant field is no longer harmonic, i.e. it does not satisfy the (homogeneous) Laplace’s equations. Furthermore, the drift and spatial distribution of the ions (the space charge) are governed locally by the electric field, which itself consists of the contributions provided by the geometry of the conductors with the applied voltages, and by the corona-generated space charge. The loop is therefore closed: the local electric field (and the space potential) depends on the coronagenerated space charge, while the latter distributes itself under the action of the electric field.
0 TheFranklin lnstituteOOl&W32/88 $3.00+0.00
315
.
J. N. Sheng et al.
Common examples of such ion flow fields are overhead HVDC transmission lines in corona, electrostatic precipitators and similar devices, as well as DC power equipment with gaseous insulation. Attempts to formulate computational approaches for DC ion flow field problems date back many decades. Initially, for overhead HVDC lines at least, attention was focused not on the field quantities, but primarily on power loss due to corona. This task does not necessarily require detailed knowledge of field and space charge distributions, and approximate methods were developed (1, 2) which center around calculations of only global terminal parameters such as the corona current. Recent concern with environmental impact of HVDC lines has, however, highlighted the need to study directly the field quantities, particularly the electric field strength, and the space charge and current densities at ground level (3). To complicate matters further, effects of wind which can disperse the space charge over a considerable area, must also be included. Suitable mathematical models of the ion flow field follow from physical understanding of the problem. It is immediately observed however, that the problem is strongly nonlinear. Furthermore, the space potential (or its gradient, the electric field) and the space charge are closely interdependent. Analytical solutions cannot be derived except for a few special cases. The problem is normally solved numerically through the use of a series of simplifying assumptions. Also, most commonly, iterative procedures are employed where estimates of the space potential and space charge are improved in each successive iteration. We address some theoretical and practical numerical aspects of obtaining meaningful solutions to the problem of DC ion flow fields, with particular reference to overhead HVDC transmission lines. First, the question of solution uniqueness is discussed. It is shown that, for simple cases at least, rigorous general proofs of uniqueness can be provided. Secondly, convergence of the iterative solution process is investigated. It is again demonstrated, rigorously for simple cases, that carefully constructed iterative approaches will be convergent. Thirdly, preferred numerical computational methods are discussed and justified, and the recommended successful iterative algorithm for general DC ion flow field problems is presented. Finally, numerical examples are given for overhead HVDC transmission line geometries.
II. DC Corona in Air The phenomenon of corona in air, manifests itself as a steady unidirectional flow of ionic species, away from the conductors in corona. While details of corona processes are complicated, the fundamental mechanisms can be described as follows for the purposes of this paper. When the voltage applied to a conductor (with respect to the ground, for instance) exceeds the corona onset voltage, electrons in the highest field region (i.e. close to the conductor) acquire sufficient kinetic energy from their travel along the lines of electric flux to undergo ionizing collisions upon impact with air particles. During an ionizing collision, an electron is knocked Q-the air particle. Hence, the ionizing collisions result in creation of additional electrons and positive ions. For the case of positive conductor voltage, both the ionizing electrons and the electrons Journal
316
of the Franklin Pcrgamon
lnsiituie Press plc
DC Ion Flow Fields
FIG. 1. Flux lines in a bipolar configuration. resulting from the ionizing collisions are attracted to the anode and neutralized there. The positive ions are repelled by the anode into low-field regions. For the case of negative conductor voltage, the resulting positive ions are attracted back to the cathode and neutralized there, while the ionizing electrons and those resulting from the ionizing collision are repelled away from the cathode. However, these electrons are quickly captured by air particles to form negative ions, which in turn are repelled by the cathode into low-field regions. From this rather simplistic description of corona phenomena, the following important observations are derived : -the polarity of ions produced by corona activities near a conductor is the same as the polarity of the voltage applied to that conductor, -the ions (the space charge) are repelled from the conductor where they were produced and fill the encompassing space. For unipolar cases, such as a conductor near a grounded plane, ions of only one polarity (coincident with that of the conductor) are present. For bipolar cases, such as a conductor with positive voltage and another one with negative voltage near a grounded plane, both positive and negative ions may exist in some regions of space. The situation is depicted in Fig. 1. Positive ions dominate in the indicated region between the positive conductor and the plane, while mostly negative ions occupy the corresponding space near the negative conductor. Ions of both polarities diffuse into and recombine in the bipolar region between the two conductors. The presence of space charge near a conductor changes the electric field there. From detailed analysis by Kaptsov (4), it is evident that, after corona sets in, the electric field strength at the edge of the ionization layer is practically independent of the applied voltage when the latter exceeds (in magnitude) the corona onset level. This provides a needed but difficult boundary condition which must be implemented in the solution of the ion flow field problem. The region of active ionization (where the ions are generated through ionizing collisions of electrons and air particles) is very thin. Its dimensions are usually neglected in comparison with other geometrical parameters of the problem. Therefore, the above mentioned boundary condition is imposed at the surface of the conductor in corona. At the Vol
325, No.
Printed
3, pp. 315 ~334, ,988
in Great
Britain
317
J. N. Sheng et al. same time, it should be pointed out that the density of space charge at the surface of a conductor in corona is normally not available in practical examples of DC ion flow fields. Finally, temporal aspects of the problem need not be considered. If the geometry and the applied voltage remain fixed, a quasi-steady-state is established where electrical and physical quantities may be functions of position, but they are independent of time.
III. Mathematical For unipolar equation :
Description
of DC Ion Flow Fields
fields in air, the space potential
&I satisfies the familiar
Poisson’s
where p is the space charge density; a,, is the permittivity of free space. The corona
current
density
J satisfies the continuity
equation
V.J=O.
(2)
The current density J is produced by (steady) migration the action of the local electric field E : J = pkE where k is the ion mobility. Diffusion and wind effects are neglected
:
of space charge p under
(3)
in Eq. (3). Also,
E = -V4.
(4)
Equations (l)-(4) comprise the mathematical model of unipolar corona and show that rj and p are interdependent. The model can also be expressed as a thirdorder nonlinear partial differential equation in 4 only : v * (V’fjV~)
= 0.
(5)
The boundary conditions normally specified include the voltages applied to the conductors, with the ground taken as reference. These boundary conditions are of Dirichlet type. An additional boundary condition is also needed. For transmission line studies, this normally is of the Neumann type and is provided by the assumption (4) that the potential gradient at the conductor in corona remains at the corona onset value E,,. For the general bipolar case, including ion diffusion and wind effects, the governing equations are
J, = p+k+E-D+Vp+ +P+W 318
(7) Journal of the Franklin lnst~tute Pergamon Press plc
DC Ion Flow Fields J_ = p_k_E+DpVp_
--p-w
V* J,
= -Rip+pP/e
V* J-
= + Rip+p-/e
(8) (9) (10)
where p+ and p are the densities in space of the positive and negative ions ; J, and JP are the positive and negative ion current densities ; k, and k_ are the mobilities of positive and negative ions ; D, and DP are the diffusion coefficients of positive and negative ions ; w is the wind velocity ; Ri is the recombination coefficient of ions ; e is the electron charge. The continuity equation also holds, where J = J, + JP for the bipolar case. Of course, when p _ (or p +) is set to zero in Eqs (6)-( lo), one obtains equations that govern the ion flow field for the positive (or negative) unipolar case, including diffusion and wind effects. The presence of the gradient terms Vp, and Vp- in Eqs (7) and (8) presents the need for additional boundary conditions. Careful analysis shows that a suitable choice is provided by homogeneous Neumann-type boundary conditions for p+ and pP at the conductors in corona, i.e. Vp, = Vp- = 0 at the coronating electrodes. This selection of boundary conditions is justified as follows. Considering as an example the simple case of unipolar corona in the coaxial cylindrical geometry, the expression for p is of the form (11)
P=V/;i+Rg’rj where A and B are constants ; r, is the radius of the inner conductor r is the radial coordinate.
;
Near the inner conductor, i.e. for r z r,, p is insensitive to radial location and remains approximately independent of r , . Consequently, close to the conductor in corona, dpldr % 0. This conclusion also applies for more complicated cases (5). For example, for the geometry of a conductor r, = 0.00056 m in radius, placed H = 0.25 m above a grounded plane (H/r = 446.4), Fig. 2 shows profiles of p as a function of position along the vertical line of symmetry through the conductor center. To facilitate comparisons, the abscissa in Fig. 2 is normalized to the conductor-plane distance H-r,, while the ordinate is expressed as per unit of the space charge density pr at the conductor surface for the respective profile. Wind and diffusion effects are not included. Figure 2 clearly shows that for nominal transmission line operating voltages T/ which are close to the corona onset voltage V,,, p near the conductor surface is indeed insensitive to location and the homogeneous boundary condition, Vp = 0 at the conductor surface, is satisfactory. The situation worsens considerably when corona intensity increases (v >> V,,). This is also true when the problem geometry Vol. 325, No. 3, pp. 315-334, 1988 Printed in Great Britam
319
J. N. Sheng et al.
D. Normalized
FIG. 2. Normalized profiles of p along line of symmetry, H/r = 446.4
0
I
I 0.2
I 0.4
I 0.6
I 0.6
I
1.0
0. Normalized
FIG. 3. Normalized profiles of p along line of symmetry, H/r = 20. deviates significantly (H/r = 20) from that of typical overhead transmission lines, see Fig. 3. For normal HVDC ljnes however, the just-discussed boundary condition on p is satisfactory (5) and permits consideration of ion diffusion effects in Eqs (7) and (8). Except for very simple geometries (infinite parallel plates, concentric spheres, infinite coaxial cylinders) the equations governing the ion flow field problem cannot be solved analytically. Numerical solutions for practical cases, such as overhead HVDC transmission lines, typically include the following simplifying assumptions :
(4
the thickness of the active ionization layer around the conductor in corona is neglected ; @I for applied voltages above corona onset level, the surface potential gradient of the coronating conductor equals the onset value E,,, and is independent of the applied voltage ;
320
Journal ofthe Franklin
lnst~tutc Pergamon Prcsq plc
DC Ion Flow Fields
(4 (4
mobilities k, and k_ of positive and negative ions, respectively, are constant, independent of the local electric field strength, but k, may be different from k- ; quasi-steady-state conditions exist, i.e. corona activities and wind velocity are independent of time.
Assumption (a) is valid for practical geometries. Assumption (b) originally proposed by Kaptsov (4), has been investigated thoroughly and confirmed to be satisfactory (6, 7), especially when the applied voltage is not very high above the corona onset level. Additional complications introduced by variable ion mobilities are normally not warranted. Therefore, assumption (c) is retained in this paper. An equivalent average value is used for k, and a different one for k_. Assumption (d) permits omission of high frequency components and electromagnetic field interactions, and is valid for stable corona conditions and normal wind speeds. IV. Uniqueness of Ion Flow Field Solution The question of uniqueness is first investigated for the simple case of unipolar corona with ion diffusion and wind effects neglected, see Eq. (5). The ideas presented are extended to more complicated cases. The mathematical properties of Eq. (5) were investigated by Atten (8), who showed that the solution for q!~is unique and the problem is properly posed in the special case when the boundary conditions consist of: -specified potentials at the conductor (applied voltage) and ground (reference zero voltage) ; -specified space charge density pe at the emitter, i.e. the conductor surface.
or
The above set of boundary conditions is not normally available in practice (3). Rather, the value of the electric field at the surface of the conductor in corona, instead of the space charge density, is provided by Kaptsov’s assumption. In order to show uniqueness of solution for Eq. (5) in the broader case which includes Atten’s and practical situations, consider a volume R bounded by surface S consisting of (m + 1) surface sections, Fig. 4. Boundary sections S,, S2,. . . , S,
FIG. 4. Volume R bounded Vol. 325, No. 3, pp 315-334, Printed in Great Britain
by surface S.
1988
321
J. N. Sheng et al. are the conductor surfaces, and S, is the outside boundary (including possibly ground) on which potential is taken as zero. If the solution to Eq. (5) in R is not unique, then at least two potential functions 4, and @2 exist which are dissimilar in R but both satisfy the given boundary conditions. The difference o of the two functions, w = 4, - cj2, satisfies the homogeneous boundary conditions and Eq. (5) in R, i.e. V * (V’oVo)
= 0
inR
o=O Using Green’s second identity
(12)
0nS.
yields
(W*Y+VV(D*V’I’)dK! sR
=
is
cDV’I’.dS.
(13)
The direction of dS is outward and normal to the surface S, and 0 and Y are arbitrary functions which possess continuous second partial derivatives. Setting 0 = o and VY = V*oVo in Eq. (13) gives
s
V20(Vw)* dQ =
(144
WV*mgdS.
R
The integral over surface S consists of contributions from all parts S,, S2,. . . , S, and S,. For infinite or semi-infinite regions R, So extends to infinity, but integrals are still bounded since 4 (and therefore, o) decreases monotonically with the distance r from conductors. Furthermore, since w = 0 on boundaries, the righthand side of Eq. (14b) is zero, i.e. n
J
V20(Vo)*
dQ = 0.
R
In Eq. (14b), the term V20 is related to the space charge density p (Poisson’s equation), which is nonzero but does not change sign in R (unipolar corona). Therefore, Eq. (14b) holds only if (Vu)* = 0 everywhere. This in turn means that o is constant everywhere. But, o = 0 on the boundaries ; hence, w = 0 everywhere. Consequently, 4 1 = C#I*and the solution of Eq. (5) in !A is shown to be unique. It is interesting to observe that no restrictions were placed in the above discussion of Eq. (5) on the kinds of boundary conditions. Therefore, in general, the solution is unique when values of potential 4, its derivative V4 and/or space charge density p are specified on the boundaries. The same conditions may be specified on all boundaries, or different conditions may be specified on different parts. Of course, the total number of conditions must correspond to the order of the equation. However, if values of only V# and p are given, then the solution will be unique within an additive constant provided the boundaries remain equipotential surfaces, i.e. 0 = 0 on S. The above discussion cannot be easily extended to the more difficult general case of unipolar corona including ion diffusion and wind effects. Equation (5) would then be of the form : Journal
322
of the Franklin
Institute plc
Pcrgamon Press
DC Ion Flow Fields V*(V’$4’c$)-F=
where
0
W>
F = k’, [V * (DVp) - VQw)] ; lJ D is the ion diffusion
coefficient.
It should be observed that 4 and p are not independent of each other. Rather, as discussed in the Introduction, the space potential 4 depends on the coronagenerated space charge p, while the latter distributes itself under the action of the total local electric field E, where E = -V4. For this reason, a general argument in favor of uniqueness of solution to Eq. (15a) is not straightforward. However, if one assumes for the moment that p is known at some point in the computational process, then Eq. (15a) can be rewritten as V . (V’@‘cj)
= F
(15b)
where the right-hand side is known. In that case, the natural solution, which satisfies Eq. (15b) with F = 0, i.e. as in Eq. (5), is unique for reasons outlined earlier. The particular solution is totally determined by the function F and is also unique for a given F (and p). The available Dirichlet boundary conditions applied to the total solution (natural plus particular) serve to determine the unknown constants. Conversely, when one assumes that 4 is known in a particular computational step, then it is possible to show that p can be determined uniquely from Eq. (15~) : ; [V * (DVp-V
* (pw))] = F1
(154
where F, = V. (V’cjVcp) is known. This suggests a two-step iterative process which yields unique solution for the respective quantities (4 or p) at each step. The iterative process is discussed subsequently. Uniqueness of solution for the bipolar case is, unfortunately, a still more difficult task. In the positive and negative regions, see Fig. 1, the space is occupied by primarily unipolar space charge of respective polarity. Therefore, arguments pertinent to the unipolar case apply. For the bipolar region between the conductors in Fig. 1, uniqueness proofs are not available. Note, however, that the bipolar region approaches but does not include the grounded plane, see Fig. 1. Consequently, uniqueness results for the respective unipolar cases are sufficient when interest is focused on the ion flow field quantities at the plane. The general discussions regarding uniqueness of solution and kinds of boundary conditions needed to show that successful computational approaches are possible for the practical situations where values of 4 (applied voltages) and of V4 (from Kaptsov’s assumption) are given as boundaries. As often happens, several approaches have been attempted with success, even without these theoretical results. Vol 325, No. 3, pp. 315-334, Printed in Great Britain
1988
323
J. N. Sheng et al. V. Iterative Numerical Solutions and Their Convergence The equations governing the ion flow field clearly reproduce the physical interactions between the field quantities. The equations can be solved analytically only in very simple unipolar cases (9). For all other situations, numerical procedures must be used. Due to the inherent coupling between the field quantities, the most successful numerical approaches have been iterative in nature. For example (l& 12), space potential 4 is computed from Poisson’s equation [Eq. (1) or (6)] with some assumed or previously computed values of the space charge. The iteration cycle is closed by computing the space charge from the continuity Equation (2) using the most recent values of d. One must ask whether such an iterative process is convergent. This question is investigated in this paper, first for the simplest case of unipolar DC corona without ion diffusion and wind effects. Experience gained from this simple exercise is extended to more complicated situations, including the general bipolar case with diffusion and wind contributions. Keeping in mind that unique solutions to the ion flow field problem exist, as demonstrated in the previous section, convergence of the iterative process is established in three separate steps. Firstly, it is pointed out that, during each iteration cycle, the solution for 4 from Eq. (1) or (6) is unique for previously computed values of the space charge. Secondly, solutions for the space charge in each iteration step have been proven to be unique. Finally, it is argued that, with unique solutions in each iteration cycle, the process should converge to the solution shown to be unique earlier. (a) Uniqueness of solution for 4 An equation which can be used to compute 4 is the Poisson’s equation [Eq. (1) (or (6) for bipolar fields)] when space charge density p is available. Uniqueness properties of Poisson’s equation are well known and need not be reported here. It is also known that calculations of 4 from the continuity equation, if needed, yield a unique result as well. Therefore, 4 can be computed uniquely. (b) Uniqueness of solution for p Considering the simple unipolar stitution of Eq. (3) into (2) gives
case, Eqs (l)-(4),
again
as an example,
v - (pV4) = 0.
sub-
(16)
At this point, 4 is assumed to be known and only one boundary condition on p is needed. This boundary condition is pe, the value of p at the surface of the conductor on corona, assumed to be available from previous calculations. If the solution of Eq. (16) for p is not unique in R, then there exist at least two functions, p, and p *, which are dissimilar in R but satisfy Eq. (16) and the boundary condition pe. The difference 0 of the two functions, 0 = p , - p2, satisfies Eq. (16) and homogeneous boundary conditions : V*(OV4)=0 O=O
in R (17)
0nS. Journal
324
of the Franklm Pergamon
Institute Press plc
DC Ion Flow Fields Using applying
Green’s second Eq. (17), gives
identity,
sR
Eq. (13), and setting
OVO V@dn
=
ss
O’gds
@ = 0, VY = OVb,
= 0.
and
(18)
The right-hand side of Eq. (18) is zero since 0 satisfies the homogeneous boundary conditions. Also, 0 decreases monotonically with distance, so the integral is bounded. Expansion of Eq. (17) yields vo*v4 Application
= -ov*(p.
(19)
of Eq. (19) to (18) gives
s
02V2d dR = 0.
(20)
n
Since V’4 is related to p (p # 0) and V’4 does not change sign in Q (unipolar corona), Eq. (20) can be satisfied only if 0 = 0 everywhere. Consequently, it is shown that p, = p2 everywhere, i.e. the solution of Eq. (16) for p is unique. For the unipolar case including ion diffusion and wind effects, Eq. (16) is modified to V.(kpV4+DVp-pw)
= 0.
(21)
The presence in Eq. (21) of the term Vp, related to the diffusion component, introduces the need for an additional boundary condition, as discussed earlier. This boundary condition is given as vp = 0
(22)
at the conductor surface. Careful analysis (10)of Eq. (21) indicates that the ion flow field problem is properly posed and uniqueness of solution can be demonstrated (without the wind term at least) when values of both the potential Q, and space charge p are given at the conductor and ground. Clearly, this is not the case here. Potential values are specified at the required boundaries, but space charge values are not. Instead, related boundary conditions [Eq. (22), Kaptsov’s assumption] are available, which do not render the problem properly posed. Nevertheless, physical understanding and observation of corona indicates that unique solution exists, even though mathematical proofs have not been formulated. For the bipolar DC ion flow field case where ion recombination is an important process, again, no general proof of uniqueness is available. When attention is confined to ion flow field quantities at the ground plane (see Fig. 1) comments on uniqueness for the respective unipolar regions are applicable. It is also clear that the Poisson’s equation yields a unique answer for p, but such calculations are usually not appealing from the numerical viewpoint, since second derivatives of 4 would have to be computed numerically (a Finite Element Method using third-order spline functions may, nevertheless, provide satisfactory results). Vol. 325, No. 3, pp. 315-334, Prmted m Great Britain
1988
325
J. N. Sheng et al. (c) Convergence of the iterative process Having studied the equations for 4 and p separately in each iteration, one final detail needs to be resolved to assure convergence of the entire iterative process. In most practical cases, the boundary conditions consist of applied voltages and the value of the electric field strength E,, at the surface of conductor(s) in corona (Kaptsov’s assumption). Use of the applied voltages is sufficient to compute 4 from Poisson’s equation with p available from the previous iteration. However, computation of a new p from the continuity equation to complete the iteration cycle, requires the knowledge of pr, the space charge density at the surface of conductor(s) in corona. This is of course, the information assumed to be known by Atten (S), but not available explicitly in normal practical situations. The boundary condition normally known in practice is, as pointed out earlier, provided by Kaptsov’s assumption. A disparity therefore exists, which has occupied researchers for many decades (l-4, 10-16). The difficulty may be resolved by using the value of the electric field E,,, provided by Kaptsov’s assumption to adjust previously computed per based on the following guidelines (5). Physical understanding of the ionized field leads to the conclusion that pTmp, the computed value of pe at a point on the conductor surface, is too electric field strength at that great (in magnitude) when EyOmp,the computed vicinity, is smaller than E,,,. Therefore, ~~~~~ should be decreased (in magnitude). Conversely, the magnitude of prmp should be increased if Epmp > E,,. The details of adjusting prmr are expected to be problem-dependent and to influence the rate of convergence of the iterative process. Several successful correction methods have been developed (5, 11-14).
VI. A New Solution Algorithm The new solution algorithm proposed in this paper is iterative in nature and follows the general outline presented in the previous sections. The region of interest D, is subdivided into small (triangular) elements using an automated grid generating program. Applied voltages are defined at the boundaries and E,, is computed for conductors in corona. For unbounded regions, such as is the case for overhead HVDC transmission lines, the region of interest is bounded by an artificial auxiliary surface surrounding the line. Measurements under existing lines show that the ion flow currents produced by corona on line conductors are concentrated mainly in the area under the line wires. In particular, these currents are very small in regions far from the line, even in the presence of normal winds. For this reason, it is possible to place the artificial auxiliary surface at a convenient location around the line, and to consider in the calculations only the effects of space charges within this surface. The iterative process of calculations is started by assigning initial values to p+ and p_ at the grid nodes to be used in Poisson’s equation. An important modification, however, is introduced at this point to facilitate improved numerical performance of the method. In traditional approaches, Poisson’s equation is normally solved for the space potential 4, using previously obtained values of the space charge density p. A method such as the Finite Difference Method (FDM) or Journal
326
of the Franklin Pergamon
Institute Press plc
DC Ion Flow Fields the Finite Element Method (FEM) would be utilized to compute 4. At the same time, it is observed that computation of p from the continuity equation requires the knowledge of not the space potential 4, but of its spatial derivative (gradient) E = (- V4). Such a procedure is prone to significant numerical errors and subsequent difficulties, particularly close to the conductor surface where variation of the potential to be differentiated is very rapid with position. In order to avoid the above outlined problems, the Charge Simulation Method (CSM) is used in this paper to compute the needed components of the electric field strength vector E directly from Poisson’s equation in each iteration step. The obtained values of E are of course unique since uniqueness proofs discussed for 4 also apply to E. The new values of E are then used to compute new space charge density in preparation for the next iteration step. The Weighted Residual Method (WRM), applied to Eq. (23), is utilized for this purpose. Equation (23) is derived by inserting Eq. (7) into (9) for calculation of the positive space charge density, i.e.
v-[(k+E+w)~+ --D+Vp+l+(Kp-le)p+ = 0.
(23)
A similar procedure is followed for the negative space charge density p _, using Eqs (8) and (10). The flowchart of the iterative algorithm proposed in this paper appears in Fig. 5 for the case of overhead HVDC transmission lines. The CSM and the WRM are described next.
VII. The Charge Simulation Method
(CSM)
for Poisson Fields
In the presence of corona the electric field strength E in space consists of two components (17), the electrostatic field E,, and the field E, due to the coronagenerated space charges. The electrostatic field E,, is produced by the voltages applied to the line conductors, with respect to ground. This field E,, can be computed using the CSM and the Dirichlet boundary conditions on the potentials at the conductor surfaces, and zero potential at the ground. According to the CSM approach, the conductors and their images (used to represent the influence of the ground) are replaced by their simulating charges located outside the region of interest (i.e. within the conductors and their images). The values of these simulating charges are computed using the potentials at the boundaries. Then, the electrostatic field E,, is obtained as the summation of the contributions from all the simulating charges. In particular, the component E,, does not depend on the corona-generated space charges. The second component E,, can be computed easily from the point form of Gauss’ law (18). For example, the contribution dEcT at point T in Q provided by the simulating line space charge dq, at node U is : d&r
= U/271~1)(dqulrdru,
(24)
where rUT is the distance between T and node U ; ruT is the unit vector from point T to node U. Vol. 325, No. 3, pp. 315-334, Printed in Great Britain
1988
327
J. N. Sheng et al. Begin 1) Input:
Line configumtion I
Grid autosubdivision I Input:
7,
initial
Line
value
‘.
lteratlons
voltage of p,
V; wind velocity
at conductor
Permissible
Culculatep,
crmrs
surface,
i+; K, RI, D contmlled
of ions; no. of
E, ,f,
E of ion flow
field
by CSM
I output:
E, J.p
’ FIG.
2‘
5. Flowchart of the iterative algorithm.
Superposition of the contributions from space charges at all nodes gives the electric field strength vector EcT at point Tin Q, due to all corona-generated space charges. The total electric field strength ET at point Tin R can then be obtained as the superposition of the just calculated components EcT and the electrostatic field E,, at the same point T,using simultaneously the simulating charges (and their images) both in the conductors and in the space outside the conductors. The space charges p+ and p_ are updated in each iteration. Since numerical differentiation is not required, the computed results are of greater accuracy than would be possible with FEM.
328
Journal of the Franklin lnstltute Pergamon Press plc
DC Ion FIow Fields A further advantage of the CSM is the fact that E can be calculated directly at all points in the region R, not only at the grid nodes. With other methods (FEM, for example) values of E are available only at the nodes and numerical interpolation would be required to obtain information at other locations.
VIII. The Method of Weighted Residuals
(MWR)
for Calculation of p
To complete the iteration cycle, the space charge distributions p + and p at grid nodes must be computed from the knowledge of the electric field strength vector E. The method of weighted residuals applied to Eq. (23) is utilized for this purpose using the triangular grid. For example, for Eq. (23), the space charge density p+ at any point within the element (triangle) k can be approximated as
(25) where rick)is the number of nodes in element k (i.e. 3) ; p T!. is the value of p + at node i of element k ; Nik) is the shape function i in element k. The approximation c+ to p+ is obtained by collecting contributions elements of the grid. As 3+ is not the exact solution of Eq. (23), a non-zero R will exist, where
from all residual
R = V*[(k+E+w)F+-D+VjY+]+(Rip”_/e)P+.
(26)
Setting the integral of the residual R over Q to zero in Eq. (26) leads to a set of algebraic equations which can be solved for the nodal values of ,i?+, i.e. jQWRd,=jn
where W is Through a collection linear shape used defines
W{V*[(k+E+w)F+
the weighting function Eq. (27), the continuous of contributions from functions were found the weighting function W(k) =
-D+Vp+]+(R,?_/e)p”+}dS2
= 0
(27)
for each element. global field problem is expressed in terms of individual elements. Triangular elements and to be sufficient. The element matching method components as : 0
outside element k
1
inside element k.
In Eqs (26) and (27), p_ is the most recent available approximation to p _. A procedure similar to that outlined above is applied for calculation of new estimates of p-. The values of the electric field strength E, at the conductor surfaces, available from the previous iteration are also compared with the number E,, given by assumption (b). Deviations of the computed values of E, from Eon are used to Vol. 325, No. 3, pp. 315-334. Printed in Great Britan
1988
329
J. N. Sheng et al. increase or decrease p+ and p_ appropriately to speed up convergence. (See previous discussion.) It should be mentioned that computations of d+ and fi_ followed the frontal approach, starting from the conductor surfaces and propagating outward. The rate of convergence was found to be quite good. Furthermore, not only assumption (b), i.e. the value of E,, at the surface of the conductor in corona, but also other related boundary conditions (corona onset voltages, charge densities at conductor surfaces) can be employed in the outlined computational method.
IX. Examples The excellent computational accuracy of the proposed method is first demonstrated with the aid of the simple coaxial cylindrical geometry, for which analytical closed-form solutions exist. Then, calculated results for a unipolar line configuration are compared with scale-model measurements. Finally, the method is applied to two full-scale bipolar line cases for which experimental data are available. Calculations were performed on Honeywell DFS8/52 and Prime 550 computers. Both the computational time and the number of iterations required for convergence were less than those for other published methods (15, 19).
Example 1 : coaxial cylindrical geometry The coaxial cylindrical geometry consists of an inner conductor 0.1 cm in radius energized to 25 kV with respect to a grounded outer cylinder of 1.794 cm radius. Using the Peek’s formula, the computed corona onset gradient E,, at the surface of the inner conductor is 58.85 kV cm-‘. The computed corona onset voltage V,, is 16.99 kV. The region of interest between the inner and outer electrodes is automatically divided into 480 elements with 260 interior nodes. Figure 6 shows
‘0-P i
86-i
I i :
4-i FIG. 6. Relative computational
330
o---o FEM -0
CSM
r (cm) errors in E for the coaxial cylindrical geometry.
Journal of the Franklm lnst~tute Prrgamon Press plc
DC Ion Flow Fields the relative errors in E with respect to the analytical solution when only the traditional FEM is used, and when the recommended CSM is employed. The vastly superior performance of the recommended CSM over the FEM is clearly evident in Fig. 6, particularly in the area of main interest, i.e. close to the inner conductor where the field strength is high. Near the outer cylinder, in the low field strength area, unsuitably large element size leads to significant errors. Hence, strategic node location particularly near boundaries, seems to be of importance in the CSM. Example 2 : unipolar DC line model The unipolar DC line model consists of a thin conductor 1.78 mm in diameter suspended 27.8 cm above a ground plane. The applied voltage is +45 kV DC, the measured corona onset voltage is Y,, = 35 kV and the measured positive ion mobility (14) is k, = 1.2 cm* V s-‘. Separate experiments were performed with wind speeds w of 0 and 3.1 m s- ’ ; the wind direction was perpendicular to the model line. Calculations were conducted using parameter values corresponding to experimental conditions. Figure 7 shows the measured and calculated lateral profiles, at ground level, of the electric field strength E (normalized to its maximum measured value E,,,,,) and of current density J (similarly normalized to J,J. Agreement between measured and calculated data is very good and points to the excellent accuracy of the proposed computational method. Unlike in other published approaches (2), the accuracy does not deteriorate in the presence of wind, see Fig. 7. The effect of lateral wind is, of course, to distort the profiles of E and J and to shift them in a downwind direction, see Fig. 7. The proposed computational
- -
+
0 +
-2
-I
ca1cu10ted,w=3. I rn s-’ Measured,w=O Measurfxl.w=3.1 0
m f’ mrl
I
2
L/H
FIG. 7. Measured
Vol. 325, No. 3, pp. 3lM34, Printed in Great Britain
and calculated
-2
-I
0
I
I 2
L/H
lateral ground-level DC line model.
profiles of E and J under a unipolar
1988
331
J. N. Sheng et al.
--w=3,lms-1 0
I -0.4
I
I
I
I
-0.2
0
0.2
0.4
L (ml
FIG. 8. Calculated equipotential
contours for a unipolar DC line model.
technique also permits accurate and detailed investigation of other quantities and their sensitivity to wind. Equipotential surfaces around the coronating conductors, see Fig. 8, can be calculated directly from the applied voltage and the calculated space charge using the CSM.
Example 3 : full-scaIe bipolar f 600 kV DC line In Fig. 9, the computed lateral ground-level profile of E is compared with available (20) measured data and previous calculations for a + 600 kV DC transmission line. The good agreement between the results computed in this paper and the measured data in Fig. 9 confirms the accuracy of the presented method. The corona parameters used for calculations in this paper were taken from Ref. (20).
X. Conclusions The discussions presented in this paper confirm that unique solutions exist for the nonlinear ion flow field problem. For simple cases, rigorous proofs are formulated. Analysis of more complicated cases reveals that, due to the nature of the available boundary conditions, the problem is not always properly posed. Nevertheless, iterative algorithms with good convergence properties can be developed for general cases. Physical understanding of the problem is of immense importance in formulating the algorithms and supplements the analytical treatment. The recommended algorithm presented in this paper is constructed to utilize the useful features of some existing approaches, such as the Charge Simulation Method and the Method of Weighted Residuals. The algorithm is relatively simple to implement and is computationally efficient. It also appears to be robust and suitable 332
Journalof the Franklin Pergamon
Institute Press plc
DC Ion Flow Fields Calculated:
Measured:
FIG. 9. Comparison
of computed and measured lateral ground-level profiles of E under a f 600 kV bipolar DC line (20).
for the solution of a wide range of both unipolar and bipolar ion flow field problems. The numerical examples included in the paper confirm the good engineering accuracy of the method. References (1) V. I. Popkov, “On the theory of unipolar DC corona”, &ektrichestvo, No. 1, pp. 3348, 1949. (English technical translation TT 1093 by National Research Council, Canada.) (2) V. I. Popkov, “The theory of bipolar corona on conductors”, Izv. Akad. Nauk SSR Otdelenie Tekhni. Nauk, Vol. 4, pp. 433-448, 1948. (English technical translation TT 1074 by National Research Council, Canada.) (3) W. Janischewskyj, P. Sarma Maruvada and G. Gela, “Corona losses and ionized fields of HVDC transmission lines”, Paper 3609 CIGRE 1982, Paris, l-9 Sept. 1982. (4) N. A. Kaptsov, “Electrical Processes in Gases and in Vacuum” (in Russian), OGIZ, Moscow, 1947. (5) G. Gela, “Computation of ionized fields associated with unipolar DC transmission systems”, Ph.D. Thesis, Department of Electrical Engineering, University of Toronto, Apr. 1980. (6) R. T. Waters, T. E. Richard and W. B. Stark, “Electrical field measurements in DC corona discharge”, ZEE Conf Publ. No. 90, pp. 188-190, 1982. (7) M. Khalifa and M. Abdel-Salam, “Calculating the surface fields of conductors in corona”, Proc. IEE, Vol. 120, pp. 15741575, 1973. Vol. 325. No. 3, pp 315-334, Printed in Great Britain
1988
333
J. N. Sheng et al. (8) P. Atten, “Etude mathematique du probleme du champ electrique affect& par un flux permanent d’ions unipolaires et application a la theorie de la sonde froide”, Doctoral Thesis, Faculty of Science, University of Grenoble, France, Nov. 1969. (9) N. N. Tikhodeev, “Differential equation of unipolar corona and its integration in simpler cases” Zh. T&h. Fiz. (in Russian), Vol. 25, No. 8, pp. 1449-1457, 1955. (10) N. J. Felici, “Recent advances in the analysis of DC ionized electric fields”, Direct Curr., Part I, pp. 252-260, Sept. 1963 ; Part II, pp. 278-287, Oct. 1963. (11) P. Sarma Maruvada and W. Janischewskyj, “Analysis of corona losses on DC transmission lines, Part I : unipolar lines”, IEEE Trans. Pwr Appar. Syst., Vol. PAS-88. No. 10, pp. 1476-1491, 1969. (12) P. Sarma Maruvada and W. Janischewskyj, “Analysis of corona losses on DC transmission lines, Part II : bipolar lines”, IEEE Trans. Pwr Appar. Syst., Vol. PAS-88, No. 10, p. 14761491, 1969. (13) C. Wen and Z. Yan, “Calculation of ion flow fields in air under HVDC transmission lines”, Int. Conf. Props and Applic. Dielectric Materials, Xi&, China, pp. 754757, Jun. 1985. (14) B. L. Qin, J. N. Sheng, Z. Yan and G. Gela, “Accurate calculation of ion flow field under HVDC bipolar transmission lines”, Paper 86TD513-6, ZEEE/PES 1986 Transmission and Distribution Conf. and Exposition, Anaheim, California, 1419 Sept. 1986. (15) M. Abdel-Salam, M. Farghally and S. Abdel-Satter, “Finite element solution of monopolar corona equations”, IEEE Trans. Elect. Znsul., Vol. EI-18, pp. 110-l 19, 1983. (16) T. Takuma, T. Ikeda and T. Kawamoto, “Calculation of ion flow fields of HVDC transmission lines by the finite element method”, IEEE Trans. Pwr Appar. Syst., Vol. PAS-loo, pp. 48024810, 1981. (17) J. N. Sheng, Z. Yan and B. L. Qin, “Uniqueness of the solution of DC ion flow field”, Internal research report (in Chinese), Department of Electrical Engineering, Xi&n Jiaotong University, China, Oct. 1985. (18) L. M. Magid, “Electromagnetic Field, Energy and Waves” John Wiley, New York, 1972. (19) W. Janischewskyj and G. Gela, “Finite element solution for electric fields of coronating DC transmission lines”, IEEE Trans. Pwr Appar. Syst., Vol. PAS-98, pp. lOOO_ 1012, 1979. (20) “Transmission Line Reference Book HVDC to + 600 kV”, Electric Power Research Institute, Palo Alto, California, 1982.
334
Journal
of the Franklm Pergamon
Institute Press plc