Journal of Electrostatics, 26 (1991) 65-79
65
Elsevier
Electrostatic field calculations P. Hammond Department of Electrical Engineering, Southampton University, Southampton S09 5NH, UK (Received May 26, 1990; accepted in revised form September 25, 1990)
Summary Field calculations are required in electrostatics in order to determine the energy and force distributions. These calculations are rendered difficult because the geometry of the field is generally only loosely connected with the shape of the conductors and insulators. The paper approaches the problem by analysing the inherent geometrical properties of the field. These properties are then used to construct a simple computer program based on graphics rather than algebraic equations. One of the useful properties of the method is that it provides upper and lower bounds to the exact field solution.
1. Introduction
Interest in electrostatic fields has a long history which goes back to the investigations of Franklin, Priestley and Cavendish, amongst others, in the middle of the 18th century. Initially these investigations were chiefly experimental and were based on the behaviour of electrostatic friction machines. They led to the idea of electricity as a fluid associated with matter. It was held that this fluid could be stored in capacitors and transferred from one body to another through conductors. Materials were classified as conductors and insulators in accordance with their electrical properties. It was realised that geometrical shape was important and an immediately useful application of this knowledge was Franklin's invention of the pointed lightning conductor, which made a tremendous difference to the safety of tall buildings which had previously been subject to repeated damage in thunderstorms. The discovery of the inverse square law of force between droplets of the electric fluid ( =electric charge) directed attention to the action at a distance as well as the contact action of charge. Gilbert in 1600 had formulated a field theory for magnetic action at a distance, but had thought that electrostatic action was always due to contact between electric particles. Newton, moreover,had consciously avoided the idea of a gravitational field. This meant that the electric field did not attract attention before Faraday's researches in the first half of the 19th century. 0304-3886/91/$03.50 © 1991--Elsevier Science Publishers B.V.
66
Meanwhile the analogy between gravitation and electrostatics based on the inverse square law led to a rapid development of a mathematical theory in terms of the Newtonian potential. The theory was systematised by the use of differential equations, and electrostatic problems came to be regarded as a special class of phenomena governed by Poisson's and Laplace's equations. Before the advent of computers solutions for these equations could be found only in terms of closed functions and this restricted their application to highly symmetrical structures. Although such solutions offered valuable guidance, they could not be applied directly to most engineering problems. With present-day computers these restrictions have been removed and it is now possible to obtain numerical information about electrostatic fields in complicated three-dimensional structures. From the point of view of obtaining solutions to the differential equations the subject can be regarded as needing no further investigation. However, such an approach does not necessarily meet the needs of engineers wishing to design and analyse electrostatic systems of increasing complexity and sensitivity, where two types of problem present themselves. The first type concerns the relationship between the physical behaviour of a system and the mathematical symbols used to describe that behaviour in terms of differential equations. The assertion that the equations govern the phenomena is no immediate help to an engineer, who regards the equations as tools rather than governing principles. What are needed are 'user-friendly' concepts expressed in 'user-friendly' notation and neither Poisson's nor Laplace's equations pass that test. The second type of problem relates to the computer. There may be no difficulty in obtaining numerical information, but this is only a step on the way. Until the handling and control of the information can be undertaken in a 'userfriendly' manner there is little benefit from increasing the computational power, nor does the addition of a post-processor remove the difficulty. What is needed is that the whole computer program should be closely related to the features which are important in the design and operation of an engineering system. Computer experts naturally tend to concentrate on logic and aim to construct programs suitable for use with a large variety of problems. Such universal tools have advantages, but engineers prefer a tool kit in which different spanners fit different nuts. As a result it is common experience that field computation is often relegated to research institutions, whereas engirieers in industry continue to rely on cut-and-dr~ methods reinforced by well-loved archaic formulae. The motive of this article is to adapt mathematical and computational methods to the physical features of electrostatic systems in order to provide userfriendly concepts which lend themselves to user-friendly computation. 2. Physical properties of electrostatic systems
The historical choice of electric charge as the primary quantity in electrical investigations has turned out to be particularly useful in describing the behav-
67
iour of electric current. It concentrates attention on the conductors carrying the current from place to place. In electrostatics the idea of charge is also useful in calling attention to the free charge stored on the surfaces of conductors and the bound charge producing polarisation of dielectrics. But in electrostatics it is not sufficient to think only or chiefly in terms of conducting and insulating materials and their associated electric charges. What matters chiefly is the relative position of the charges and this depends on the shape of the material and the spacing. Moreover the distribution of charge is very sensitive to these geometrical factors and it is generally very difficult to calculate the charge distribution. These considerations explain why the use of the potential is generally preferred in calculations to the use of charge. Electrostatic potential is analogous to height above sea level in gravitation. Every point in space can be given an assigned potential and the entire space can be mapped by equipotential surfaces, which are analogous to contour lines on an ordinary two-dimensional map. The potential gradients can be obtained from the spacing of the equipotentials and this provides important data about possible electric breakdown. For local information throughout a system the potential is a better tool than the charge. Another advantage is that the potential fits the differential equations and the numerical methods based on these equations. However, the potential distributionby itselfdoes not give allthe information required. It is after all not an independent quantity, but a property of electric charge. In terms of measurement it defines the potential energy between all the charges of the complete system and a probe of unit charge. Such a probe could be thought of as a small part of the system and this explains the usefulness of the potential approach to problems of electricbreakdown. But as we have already noticed even in such problems we need not only the value of the potential,but itsgradient and therefore the spacing between the equipotentials as well. This suggests that we need information about volumes of space rather than surfaces. This consideration is reinforced by the need to calculate the energy associated with the electrostaticsystem. Normally energy is the parameter of primary interest to engineers and particularly in problems of design it is important to know both the system energy and itsdistributionthroughout the volume. Neither the charge nor the potential provides this information, because both of these descriptionsof the phenomena concentrate attention on points in space rather than on interconnected volumes. W e need to look for a field theory rather than a particletheory, although we expect that the two types of theory will be closely related. W e have already noted that Gilbert laid the foundation of a magnetic field theory by describing the forcesbetween magnets in terms of a fielddistribution around the magnet. Gilbert's fieldswere of similar spatial dimensions as the magnets themselves, whereas Newton was chiefly interested in the motion of
68
the heavenly bodies, which could be treated as particles because of the vast distances between them. It is likelythat this led Newton to favour action between mass points, whereas Gilbert thought of the fieldas an extension of the magnet in space. In electrostaticswe need something akin to Gilbert's view and we find this in Faraday's idea of electricflux. Faraday focused on the space between charges rather than on the charges themselves. He associated the energy with that space and solved the problem of energy distribution by dividing space not into cubes, rectangles or 'finite elements', but into tubes which connected pairs of positiveand negative charges in such a way as to fillall space (Fig. 1 ). Whereas in the inverse square law calculation isbased on particlesand every particleacts on every other particle, the tubes limit the interaction to pairs of particles on each other. Calculation therefore involves the shape of the tubes rather than sets of equations based on mutual equilibrium between charges. The equations of statics are transposed into a geometrical problem in terms of tubes. Electrostatics becomes a branch of geometry. Oliver Heaviside isreported to have said that 'theengineer must divide space into tubes and not into cubes'. Ifhe did say this,and itsounds a typical Heaviside remark, he hit the nail very squarely on its head. Every tube contains a certain amount of flux,where the word flux is used by analogy with the flow of an incompressible fluid.The inverse square law describes the velocity of such a fluid in relation to its distance from a spherical source. Hence the notion of flux and the division into tubes of flux isequivalent to the notion of action as a distance according to an inverse square law of force. But whereas in the action at a distance model the directions of the forces cross each other in a confusing manner, the idea of tubes of flux presents a geometricalstructure which is easy to visualiseand is decidedly more user-friendly.It is possible to sketch the tubes and to obtain useful visual information by this means. The notion of flux leads easily to the idea of flux density over the crosssection of a tube. This involves subdividing the tubes or making each tube thin enough to allow its cross-section to become a flatsurface and the curved surfaces of the tube perpendicular to this cross-section (Fig. 2). W e are therefore
Fig. I. Tube of flux. Fig. 2. Flux filament.
Fig. 3. Slice of potentials.
69
led to the consideration of thin filaments of flux. Let us now turn our attention from the quantity of flux, which is equal to the charge at each end of a tube, to the potential difference along the tube of flux. In this way we can associate the notion of potential with the flux and therefore with the charge at the ends of the tubes of flux. Since the product of charge and potential difference is energy, we can now associate energy not only with an entire tube but with a filamentary tube throughout its length. In this way we arrive at the distribution of energy in space. One very useful feature of this way of subdividing the energy of a system is that all the subdivisions give positive energy contributions. This should be contrasted with methods of calculation based on action at a distance where the energy depends on the sign of the interactive charges. The energy associated with the tubes of flux is necessarily positive because the direction of the flux is always in the direction of the potential (downhill) gradient. This of course is inherent in the inverse square law constant called the permittivity, which is always positive. The subdivision of space into tubes of flux allows us to associate potential difference and therefore energy with the space. Could we in a similar way start with the potential difference and associate flux with it? This is indeed possible, if we start with slices of space rather than tubes (Fig. 3 ). By a slice we mean a region of space bounded by two equipotential surfaces and a ribbon joining these surfaces at their perimeters. If the surfaces are wide compared with their spacing, the ribbon can be made perpendicular to both surfaces and will therefore be a piece of tube boundary. Hence the volume will be associated with a constant charge and a constant potential difference. Thus each slice will be associated with a positive energy. A particular case of a division into slices arises when the equipotential surfaces are closed. Then the ribbon is unnecessary and the charge is the total enclosed charge rather than a part of such a charge. We now have two possible subdivisions for the calculation of energy. There is no need to use both, although there may be advantages in doing so. Conceptually the tubes are easier to visualise, but the slices have more direct reference to potential gradient. We shall show that in calculation both methods can usefully be combined. It is worth noting that in both methods there is orthogonality between the two types of surfaces, i.e. equipotentials and tube boundaries are at right angles to each other. This property gives useful guidance about the accuracy of field plots if one uses both tubes and slices and ensures their mutual orthogonality. In terms of flux • and potential difference V the energy U of a small region can be written JU=-½ JVJ~
(1)
The negative sign arises from the (arbitrary) definition of potential differ-
70
ences as negative when the gradient is in the direction of the field. The factor 1 occurs when the flux and p.d. are proportional to each other as they are in regions of constant permittivity. If the permittivity varies with the field strength this factor must be replaced by an integral related to the manner in which the energy is changed. In terms of flux density D and electric field strength E we have for a small piece of tube or slice
JtP=D JS
(2)
-JV=EJI
(3)
where gS is a small cross-sectional area and Jl a small length in the direction of the field. A small volume gv can be formed by combining $S and $1. We can write
JU=½ EDbv
(4)
and this defines a volume density of energy ½ED.This expression must be treated with caution, because Jv is not an arbitrary volume independent of D and E. The primary quantities are flux and p.d. and the subdivision of the volume must be carried out in terms of tubes of flux and slices of potential difference. D and E are mathematical limits, but haveno independentphysical properties measureable at a point. Flux always requires area and p.d. always requires length. These geometrical features are guaranteed in the notation of differential forms which treats EJl as a single mathematical entity called a 1form and D JS as a 2-form. The nomenclature is based on the dimensionality of the length. The energy density JU is then a 3-form. This notation is unfamiliar in comparison with the more usual vector notation, but it fits electrostatics far better than the representation of space by a lattice of nodal points. The interested reader is referred to a recent paper by the present author [ 1 ]. As long as the underlying geometrical structure in terms of tubes and slices is kept in mind, the idea of energy density can be used with great benefit. It firmly attaches energy to the entire volume of the system instead of concentrating on charges acting at a distance. Secondly the fact that it is always a positive quantity is a great help in computation as we shall show in the next section. Thirdly energy density is closely related to local field strength and this enables the energy .mlculation to be used to determine the field distribution. We shall show that it is more economical to determine the field via the energy rather than vice versa. Before we leave the discussion of the physical properties of electrostatic systems it is worth noting that the subdivision of systems into tubes or slices can be applied without modification to regions containing polarisable dielectrics. The flux tubes are continuous in a homogeneous dielectric and terminate as before on tbe charges located on conducting surfaces. The energy density re-
71
\,\\\\ Fig. 4. Needle-shaped cavity. Fig. 5. Disc-shaped cavity.
mains given by eqn. (1), but we have to enquire what the local values of E and D are. 'To give a definite answer it is necessary to devise an experiment in which a small probe can be inserted into a cavity in the material. It is important that such a cavity should not disturb the system and one might be tempted to think that any small cavity would meet this requirement, but this cannot be so because the disturbance will depend on the shape of the cavity. It is well-known that only two shapes leave the system undisturbed, namely a needle-shaped hole in the direction of the field (Fig. 4) and a disc-shaped hole perpendicular to the field (Fig. 5). Our previous discussion shows at once that the needleshaped hole is a filamentary tube and the disc-shaped hole a portion of a slice of tl~e field. The needle cavity hardly disturbs the dipole distribution in the dielectric either outside or inside the cavity. The disc by cutting across the dipole chains adds the local polarisation effect to the undisturbed field in the cavity but does not disturb the field of the system because the two charge layers are close together and their effects cancel everywhere outside the cavity. Therefore the needle-shaped cavity will allow E to be measured and the disc will measure ~r E, where ~r is the relative permittivity. This discussion shows that the geometrical approach via tubes and slices is closely linked with measurement. 3. Computation of electrostatic fields
The discussion so far suggests strongly that a geometrical approach to field calculation should lead to economy as well as better control of computation methods. The purpose of the present article is the discussion of the general principles underlying such computation but for a detailed description of a computer package the interested reader is referred to a paper by Sykulski [2 ]. The first step in all computer methods for field calculation is the discretisation of the region of interest. In finite-element programs the individual elements are generally tetrahedra for three-dimensional problems and triangles for two-dimensional p~oblems, because these lend themselves to simple inter-
72
polation procedures for the field between the nodal points. It is interesting that this choice of element is governed by a desire to simplify the field equations. Indeed the entire process in finite-element as well as finite-difference, boundary-element and integral methods is directed to obtain solutions for the local equations and the associated boundary conditions. The motivation for the method of tubes and slices is different because it does not start with equations but with the geometrical shape of the energy distribution. A correct discretisation into either tubes or slices already embodies the solution of the problem and requires hardly any further computational procedures except some simple addition of energies. In particular there is no need for matrix inversion or for iteration. Since those procedures are invariably cumbersome, their absence results in a drastic simplification of the computation. The reason why the discretisation into tubes and slices avoids the solution of simultaneous equations throws a great deal of light on the structure of field problems. The equations used in most methods relate to the field around nodes or smail regions of a system and therefore give information only about the local sources. Since the local field is, however, the result not only of local but also of distant sources such as those on the system boundaries, it is necessary to link the local distributions to the entire system. This naturally leads to a set of simultaneous equations. It is not surprising that the linking of large numbers of pieces of information into a consistent assemblage requires arduous computation. Against this the method of tubes and slices starts with the required linkages, because the tubes and slices already contain information about the complete system. The field geometry embodies not only the local information but the interconnections as well. The equations on the other hand treat the field as a set of particles which have to be glued together and this disregards the geometrical structure of tubes and slices which is the essential feature of the field. The price that has to be paid for the remarkable economy in computation of the geometrical approach is that it is not possible to use finite elements of arbitrary shape. The elements must be either tubes or slices and the computer must be used as a sketch-pad rather than an equation-solver. This has both advantages and disadvantages. The advantages are that the sketch-pad method is particularly appropriate in design work because it allows the designer to see at a glance how the field varies with geometrical and material changes in the system. The method is also ideal for teaching purposes, because it enables the students to obtain visual experience of field distributions. On the other hand the equation-solver approach allows the computation to be undertaken as an independent mathematical task which does not require engineering knowledge of the problem under investigation. The description of the method of tubes and slices as one requiring the sketching of fields may throw doubts on its accuracy. The seeming advantage of algebraic methods is that they produce numerical information which gives
73
apparently total accuracy. The slightly unkind reference to a 'seeming' advantage is prompted by the fact that the information may not be usable. Often recourse is had to a post-processor, which draws a field plot and so arrives at tubes and slices in a round-about manner. Nor must it be thought that the geometrical approach has no quantifiable measure of accuracy. Indeed, the use of both tubes and slices enables the accuracy to be controlled very closely. If the tubes and slices are everywhere orthogonal to each other the field distribution is necessarily correct. A particular example of this is the two-dimensional method of curvilinear squares, where the field can be obtained purely by means of a pencil and a rubber, although readers who have used that method will be aware that quite a lot of skill is required as well. Such a procedure shows that information about the field does not need the solution of equations, nevertheless the method is unsatisfactorybecause it does not give information about the sensitivityof the method to lack of orthogonality. These deficienciescan be readily overcome if the pencil and rubber are replaced by a computer with a visual display. Let us examine the accuracy of a subdivision into tubes of flux.Following an argument given by the famous J.C. Maxwell [3] we can reason as follows: Suppose we constrain the flux by inserting impervious flux boundaries into the system. Ifthese boundaries follow the stream lines (or tubes) correctly they will not disturb the flux, but if they are incorrect they will impede the flux and thereby reduce the flux for a given p.d. The ratio flux/p.d, is the capacitance and an incorrect discretisation into tubes will therefore reduce the capacitance and give a lower bound for its unknown correct value. Now the calculation of the capacitance of a set of flux tubes is a very straight-forward process involving ratios of areas to lengths. A simple program will, therefore, cause the computer to give a number for this capacitance. We now alter the tubes and call up the new capacitance. If the capacitance increases, the alteration has been beneficial. The process can be continued until no appreciable improvement is noted. Of course even then we shall not know how close our answer is, but this difficulty also can be overcome by the use of the additional subdivision into slices. Such a subdivision either leaves the field undisturbed if the slices lie along equipotentials or it reduces the p.d. for a given flux. Hence it provides an upper bound for the unknown capacitance. The slices can be improved to give a reduction of the upper bound, which will approach the correct value. Both tubes and slices used together will provide confidence limits for the true capacitance. Experience shows that the average value of the two capacitances is very stable. Hence it may be unnecessary to narrow the bounds very closely if the total energy is the chief interest. If in addition the local field distribution is important a smaller difference between the bounds may be advisable. Local accuracy can be examined by reference to the orthogonality between the tubes and the slices.
74 4. I l l u s t r a t i o n of the method
Figure 6 shows the cross-section of a tubular capacitor in which the inner electrode is offset from the outer one. The space between the electrodes is filled with a homogeneous material of permittivity e. It is required to find the capacitance per unit length. Symmetry about the horizontal mid-plane enables us to restrict the calculation to half the capacitor as shown in Fig. 7. Figure 8 shows a simple constrained flux distribution which will provide a lower bound for the capacitance
C_=2e i + g +
=7.47e
(5)
Figure 9 shows a simple constrained potential distribution which will give an upper bound C+ =2e
+~+m
=13.07e
(6)
The average value of C_ and C+ is C=10.27e. The confidence limits are _+27.3%. In these very simple hand-calculations we have assumed that all the flux is confined to straight parallel paths and that the potential gradient is uniform and follows straight lines between the electrodes. These very strong constraints can be relaxed slightly while yet retaining the simplicity of a handcalculation.
Fig. 6. Tubular capacitor.
Fig. 8. Simple flux distribution,
Fig. 7. Half-section of capacitor.
75
Fig. 9. Simple potential distribution.
Fig. 10. Modified flux distribution.
Figure 10 shows a possibl~ modified flux distribution for which C =2 e["2+ [
1 + _3+ 1 .5+o5 3 3.5÷o.
+ 17= 7.90~
The modification has improved the lower bound by 5.8%. Figure 11 shows a possible modified potential distribution for which
C+ - 2 e 2×0.5+2__><0.5 + 2×1.5+2×1.5 + 2×2.5+2X2.5 8.5
5.5
17
11
= 2~ [3.339+ 2.226+0.6679] = 12.47~
8.5
5.5
(8)
This has improved the upper bound by 4.6% by treating each of the three sections of potential distribution as two capacitors in series instead of as a single capacitor. This reduces the constraints on the potential distribution by allowing a different gradient in each capacitor. The average of the bounds is now C= 10.19~ with confidence limits + 22.4%. It will be noted that the average value has changed by less than 1%. Further improvements suitable for hand-calculations could be made, but at this stage it is better to resort to a simple computer package such as TAS (Tubes and Slices) [ 2,4 ]. Figure 12 shows a combined discretisation using 20 tubes and 8 slices and Fig. 13 shows the corresponding field map of 10 equally spaced equi-
76
Fie. 11. Modified potential distribution.
3 2 1 0
!
I
o
~
2
s
4
5
6
7
s
9
~o
Fig. 12, TA$ combined tube and slice discretisation 20 tubes and 8 slices. 6
4
2 I
0
2
3
4
5
6
?
8
9
10
Fig. 13, T A S combined tube and potential contours 10 tubes and 10 slices,
potentials and of 10 flux tubes containing equal amounts of flux. The capacitances are C_ -9.47~, C+ = 10.96~ given an average Cffi 10.21 with confidence limits +. 7.3%. The bounds have now been considerably improved, but the average value is almost unchanged. If required, improvements can be made both by altering the shape of the tubes and/or slices and by increasing their number. The TAS program provides almost instantaneous values for the capacitances, so that it is suitable for interactive design work.
77
These results were compared with finite-element calculations. It is not always realised that the normal finite-element programs for electrostatics provide upper bounds rather than the correct solution. It is, however, possible to construct dual finite-element programs to provide lower bounds also [5 ]. The TAS package has an extension which gives finite-element solutions for both kinds of bound using input from the tube and slice calculations and thereby greatly increasing the rate of convergence. For our example the finite-element calculations gave an average value C= 10.14e with confidence limits + 1%. This underlines the stability of the average value which is a feature of the tubes and slices method. Thus the average of the very crude representation of Figs. 8 and 9 is within 1.3% of the finite-element average value. The confidence limits of +_27.3% therefore give a very pessimistic estimate of the accuracy. The improved hand calculation is within 1/2% and the TAS calculation within 3/4% of the finite-element average value. Since the finite-element values are themselves subject to confidence limits of +_1%, it can be safely concluded that the tubes and slices calculations make the complication of finite-element calculation unnecessary if the capacitance of a system is to be calculated. If the field distribution is required there is little to choose in accuracy between a tube/ slice calculation for a relatively small number of tubes and slices and a finite-
5 4 3 2
t 0
0
1
2
3
4
5
6
7
8
9
10
Fig, 14, Finite-element flux and potential contours using 1000 elements. TABLE 1
C_ C+
C,,~ Confidence limits
Simple hand calculation
Improved hand calculation
TAS package
Finite-element highly discretised
7.47~ 13.07~ 10.27e +_,27.3%
7.9c 12.47~ 10.19e ± 22.4%
9.47~ 10,96e 10.21e +_7.3%
10.04~ 10.24~ 10.14e +_=1%
7~4
element calculation using a similar number of elements. For a more accurate field distribution a good compromise is to use a tube/slice solution as the input to a highly-discretised finite-element program. The improved convergence obtained by using the results of the tube/slice method makes it possible to use the finite-element method in an interactive mode on a p.c. Figure 14 shows a field plot for our example using approximately 1000 finite elements, which should be compared with Fig. 13. Clearly the finite-element field looks much more accurate, but the improvement occurs mainly in regions of low energy. Table 1 summarises the capacitance calculations. 5. Conclusions
In electrostatic problems the energy is distributed throughout the system ~nd use has to be made of a field theory to determine the energy distribution. The normal way of doing this is to describe the field in terms of partial differential equations, which then have to be solved throughout the region. In order to make use of compt, ters the region is discretised and the differential equations are replaced by a set of simultaneous algebraic equations. The solution of these produces values of field quantities at the nodes of the discretised region. We have sought to show that an alternative procedure can be developed which describes the field in terms of an energy distribution having geometrical structure. This enables the solution to be found without the use of equations. The structure allows the region to be subdivided into orthogonal systems of tubes
and slices.The capacitance is calculated either from the tubes or from the slices.The Ibrmer provides a lower bound and the latteran upper bound for the correct capacitance. The bounds can be improved by redrawing the tubes or slices.The two bounds provide confidence limits for the capacitance and if these limits are close the local field will be correct. This can be checked by looking at the orthogonality locally. It is claimed that the method has particularvalue for interactivedesign calculations and for teaching. The method shows great economy in computation. If local as well as global accuracy is required,the method can usefullybe combined with a finite-element method.
Acknowledgement The T A S computer program used in the example was written by the author's colleague Dr. J.K. Sykulski. Dr. Sykulski's contribution to the development and application of the method of tubes and slicesis gratefullyacknowledged.
79
References 1 P. Hammond and D. Baldomir, Dual energy methods in electromagnetism using tubes and slices, lEE Proc. 135, A, 3, March 1988, pp. 167-172. 2 J.K. Sykulski, Computer package for calculating electric and magnetic fields exploiting dual energy bounds, IEE Proe. 135, A, 3, March 1988, pp. 145-150. 3 J.C. Maxwell, Electricity and Magnetism, Clarendon Press, Oxford, 1892, Third Edition, Articles 306 and 307. 4 Vector Fields, Kidlington, Oxford, TAS Computer Package. 5 P. Hammond and T.D. Tsiboukis, Dual finite-element calculations for static electric and magnetic fields, IEE Proe. 130, A, 3, May 1983, pp. 105-111.