- Email: [email protected]
Low frequency electromagnetic field computation in three dimensions
COMPUTER METHODS NORTH-HOLLAND
IN APPLIED
MECHANICS
AND ENGINEERING
52 (1985) 653-674
LOW FREQUENCY ELECTROMAGNETIC FIELD COMPUTATION IN THREE D~NSIONS C.W. TROWBRIDGE Ruthe&rd
Appleton Laboratory, Chilton, Didcot, Oxon, U.K. Received 22 October
1984
This paper presents the various differential operator formulations for obtaining numerica solutions to the classical equations of electromagnetics using the finite element method. Numerical examples are presented for 2D static and eddy currents problems. Three-dimensional formulations are discussed for nonlinear static problems using the total and reduced scalar potentials with results from a recent magnet project involving field computation at the design stage and validation after const~ction. Finally a review of some of the formulations for solving 3D eddy current problems are discussed with comparative results against measurement.
1. Introduction Computational modelling is now a basic requirement for the engineer who wishes to design electromagnetic devices. The potential list of end-users of electromagnetic applicable software is enormous, including engineers working on the following applications: equipment and facilities used in physics and engineering requiring magnetic fields of precise shape, for example, charged particle focussing magnets for nuclear physics accelerators, electron beam lithography, and plasma physics containment systems; electromagnetic systems for medical diagnostics such as magnets for nuclear magnetic resonance (NMR) body scanning; probes and evaluation of defects in non-destructive testing equipment; and a multitude of electrical machines and power handling devices of all kinds. A challenging list of problems. In all these cases the designer starts with a performance specification; he then creates from his experience and skill a model which he hopes will go some way toward matching the given specification. In order to test his ideas he needs to predict the performance of his design; first by relatively simple analytic and numerical procedures, and secondly by sophisticated computer techniques, taking account of as many of the design parameters as possible. For many applications two-dimensional simulations are sufficient, but there is an increasing demand for complete three-dimensional solutions in order to minimise the need for building expensive prototypes. Therefore, at the heart of a ‘computer-aided design’ simulation system for an electromagnetic device an algorithm for solving the electromagnetic field equations is required. 00457825/85/$3.30
@ 1985, Elsevier Science Publishers B.V. (North-Holland)
C. W. Trowbridge, Low frequency electromagnetic $eld computation
654
2. Equation for electromagnetic
fields
The classical field equations of Maxwell [l] are sufficient for the range of electromagnetic applications under consideration; these equations together with the nonlinear constitutive laws and appropriate boundary conditions define the field uniquely [2]. The canonical form of these equations are usually written as follows:
(1) (2) div D = p ,
(3)
(compatibility) div I? = 0 ,
(4)
where the vectors J? and H are the electric and magnetic field intensities and the vectors fi and B are the electric and magnetic flux densities respectively. The electric current density and electric charge density are .? and p respectively. The field intensity and flux density quantities are not independent; they are related by the constitutive laws of the material,
B=pH,
(5) (6)
D=&,
where p and E are the permeability and susceptibility respectively. In general both p and E are tensors and ,!Lin particular will be field-dependent. The electric field intensity vector g is further related to the current density J by Ohm’s law,
tTE=J,
(7)
where (T is the conductivity of the material. The above equations can be combined to produce many different forms by eliminating one or more of the dependent fields; for example if p, E and CTare constant, a single equation results for each of the field quantities, i.e. for the magnetic intensity R,
(8) At interfaces between regions with different material constants following continuity conditions: ii*(B2-&)=O, fix(tl;!-HI)=
the field equations
imply the
(9) K,
(IO)
C. W. Trowbridge, Low frequency electromagnetic field computation
A@*-D1)=
q,
655
01)
iix(&BJ=O,
(12)
where ri is the positive normal of the surface separating the two regions (1) and (2), and where K and q are possible surface current and charge densities respectively. In conducting regions (a % .F) for slowly varying time fields the second term on the right-hand side of (8) is negligible; this is the case for time-harmonic problems in which R = H0 exp(-iwt)
(13)
and where o, the angular frequency, is small, such that the wavelength is large compared with the dimensions of the problem. This is equivalent to neglecting the aD/& term in (2), and the resulting subset of equations is known as the quasi-static Maxwell equations, i.e. (2) becomes curl fi = .T.
(14)
Hence the displacement current (Jo/at) term in Maxwell’s equations will be neglected in all that follows since only static and low (power) frequency problems will be considered. At the surface of a conducting region the normal component of current density J, will be assumed to be zero since any tendency for current to flow is opposed by the electric charges accumulating at the surface. In fact, taking the divergence of (2) and using (4) the equation of continuity results, divj=
--@fat.
It is shown in the standard texts [3] that (15) implies an short relaxation times for conductors; even for a poor microsecond. In the low frequency limit approximation instantaneously to the conductor surface; so for interior
(1%
exponential decay of charge with very conductor like water it is less than a it is assumed that free charges migrate points
divj=O,
(16)
and on the surfaces
&=o.
(17)
3. Magnetic vector potentials Equation
(4) implies the existence of a vector potential
A,
B=curlA, which by Helmholtz’s theorem [4] will be unique if the divergence is specified together with the tangential component on the boundary. Since in the time-dependent case both electric and
C. W. Trowbridge, Low frequency electromagnetic field computation
656
magnetic fields are coupled, (1) and (18) imply E=-($+W),
(19)
where V is an electric scalar potential; equation in A-V results,
hence by combining
(l), (2) and (5) an equilibrium
curlbcurlA=-g($+DV);
(20)
also, from (16)
(21) Equations (20) and (21) are the defining equations for a quasi-static conductivity u and permeability ,s. It will be sufficiently general embedded inside a region of free space L& which may be infinite. The defining equation for region 0, is simpler, since u = 0,
field in a region 0, of to consider 0, to be
curlLcurlA=O; P
(22) is possible since curl H = 0 and a magnetic
indeed, further simplification may be introduced,
scalar potential
B = -grad 9,
$
(23)
which by (4) satisfies
v2+ = 0.
(24)
On the interface r12 between 0, and 0, the standard continuity must be satisfied, @“)I
=
@“)2,
(ml
=
(fix
conditions,
i.e., (9) and (lo),
(25)
*
Thus the scalar $ must be coupled to the vector A in such a way as to satisfy (25). Also at the surface f,, (17) must be satisfied, which from (19) implies
aA” av
u(
-+-& at
Hence, it appears
>
=o.
that (20) and (21) for region a,,
(26) and (24) for region J2,, together
with
C. W. Trowbridge, Low frequency electromagnetic field computation
657
interface and coupling conditions, (25) and (26) and the external boundary conditions on region a,, completely define the problem. The choice of divergence, or gauge, on A will be considered later, see Section 5.2.4.
4. Static limiting case There are three limiting cases to consider - electrostatics, The electrostatic case is straightforward, i.e. from (1)
magnetostatics
and current flow.
curl E = 0, which implies the existence of a scalar potential
V,
E = -vv, and from the divergence
(27) condition
(3),
(28)
v*Evv=p.
Hence, electrostatic problems will satisfy Poisson’s equation (28) and the electrostatic potential V will be unique in a domain L! provided V (Dirichlet condition) or aV/dn (Neumann condition) are prescribed on the boundary-the solution of a pure Neumann problem (cYV/&I defined on the boundary only) will be unique if V is defined at one point at least. For magnetostatic problems the conductors are usually prescribed, carrying a current density x. Hence the vector potential A could be used; i.e., from (14) and (18) curl L curl A = Js .
(2%
P
Equation (29) is most useful in two dimensions since in this case A has only a single component; in three dimensions there is the problem of gauge to contend with as well. Furthermore, an A formulation would require the conductors to be discretised - a formidable problem in three dimensions. It is more economical to explore the possibilities of scalar potentials, of which there are two kinds. The first kind is usually termed the total scalarpotential and is only valid for the conductor-free regions of the problem, i.e., curl fi = 0 ;
(30)
fi=
(31)
then -V$,
where $ is the total scalar potential. two parts [5],
The second kind arises by partitioning
the field H into
(32)
C. W. Trowbridge, L.ow frequency electromagnetic fit&i computation
658
where gS is defined by (33) which implies that
and (34) in which 4 is often called the reduced scalar potential (4) Cc,satisfies
[6] and is valid everywhere
v-pvt#J=o
[5]. From
(35)
and C$satisfies V.pV#
= div,uHS,
(36)
with suitable boundary conditions. Finally, for problems of electrical conduction (16) and (7)
the definition of a scalar potential follows from
divJ=Vcr-B=V.aVV=O.
(37)
Thus the potentials for all three types of problem therefore all use the same algorithms.
satisfy the quasi-Poisson
equation
and can
are straightforward.
In this
5. Numerical solutions 5.1. Two-dimensional problems Numerical solutions of the field equations for two dimensions case (20), for the single component vector potential, reduces to
and finite difference and finite element methods have been extensively used to solve (38)-the literature is very rich in examples [7,8,9]. Two examples of two-dimensional solutions are shown in Figs. 1 and 2, in which the package PE2D (Poissons Equation in 2 Dimensions) [lo] has been used to solve a plane magnetostatic problem
and an eddy current
problem
with axial symmetry.
C. W. Trowbridge, Low frequency electromagnetic jield computation
659
3t
I-
I-
I
-6.OOC 3
Fig. l(a). Permanent magnet motor. Finite element mesh generated processor [lo]. There are 2508 first-order triangles in the mesh.
5.2. Three-dimensional
by using the PE2D computer
code pre-
static problems
It is clear that solutions for electrostatic and current flow problems can readily be obtained by the use of scalar potentials and standard finite element methods- there are numerous examples in the literature [ll]. The magnetostatics case is not so clear; e.g. in Section 4 two kinds of scalar potential were defined but it can easily be seen that the presence of conductors carrying electric current in a problem poses difficulties in using a total scalar potential since it is not single-valued. The reduced potential which avoids this difficulty, has been used with some success, especially for predicting fields in the free space regions 16,121. The use of the reduced potential was investigated at Rutherford Appleton Laboratory [13] with results that indicated that, under many circumstances, erroneous fields could be computed in highly permeable regions. The reasons for this were identified with the defining
C. W. Trowbridge, Low frequency electromagnetic field computation
660 Cl.OO(lr
4.001
2
~ 001
G.OOl
I
3
!
3
I
1
-2.i)Ol
-4.OOI
-i-1
~ clril
_ i -5.000
I
0.000
Fig. l(b). Permanent magnet motor. Finite element solution displayed by contours of constant This example is the subject of a patent held by Prof. K. Binns of Liverpool University.
5.000 vector potential.
equation itself, i.e. (32), which may be expressed as
in which it is obvious that if both E?m and V4 are of similar magnitude, VC#Jwill have to be computed to the same order of accuracy as r?, if cancellation errors are to be avoided. In fact, E;i, is often calculated by analytic or by numerical quadrature to high accuracy and can have a strong spatial variation. Furthermore, errors can be exacerbated when 8 is computed by (5) in a nonlinear iteration for predicting saturation effects in iron. For these reasons it makes good sense to calculate the difference quantity direct by use of the total scalar potential. This can easily be achieved for a very large range of problems by subdividing the domain of the problem into subdomains and selecting the total potential for those regions that do not contain conductors and the reduced potential for the remainder.
C. W. Trowbridge, Low frequency electromagnetic field computation
4 II
II
II 0 0
661
i
2ll.o00-
10.000-
L \
-I
Fig. Z(a). Axi-symmetric eddy currents. The picture shows contours of induced electric current density in a molten metal experiment for investigating forces in an induction furnace. The experiment consists of mercury in a stainless steel pot surrounded by a 50 Hz coil. The problem was solved using PE2D code [lo] with a mesh consisting of 1813 second-order triangles.
Thus, for a two-region
problem the defining equations
are as follows: for region 1,
for region 2,
There is a saving in cost as well because the source field E?, does not appear in (40); in fact, when the two regions are coupled together the source term reappears for points ‘on the interface. The two regions are coupled by use of the standard interface conditions, i.e. (25), which when re-expressed in terms of potentials are
(41)
C. W. Trowbridge, Low frequency electromagnetic field computation
662
4o.oot
-
i 3o.ooc
-_ -
20.00[
-
1o.ooc
-
o.oot
-t
Fig. 2(b). Axi-symmetric
- i
eddy currents.
Lines of constant
(rA,) giving the rest&ant
($),=(iG-$),.
by quadrature
flux.
(42)
The explicit relation between 4 and 1,5is obtained by integrating interface surface r12, i.e.,
With H, determined
magnetic
from the Biot-Savart
law,
(42) along a contour over the
C. W. Trowbridge, Low frequency electromagnetic field computation
663
Permanent magnet sources can be taken into account by adding the coercive field to (32) and by using the appropriate constitutive relationship [14). Again a very efficient ~go~thm can be obtained by applying the Galerkin method to (39) and (40) and coupling the two regions together by (43) [13]. A general code for threedimensional static fields called TOSCA (TWO SCAlar potentials) has been developed at Rutherford Appleton Laboratory, based on this approach, with several comparisons made against measurement [ 141. In fact the literature of viable 3D computer implementations even for the statics case is not very extensive, see [14,15,16] for a selection. The practice of adopting successful computer codes from other ‘fields’ to electro-magnetics has been reported from time to time [17,18] this should be done with great care and with regard to the cancellation errors arising from the exclusive use of the reduced scalar potential-see (32) and comments above. That serious errors arise can be seen in Fig. 3 for results of an experiment carried out at RAL. In Figs. 4 and 5 successful 3D results using the TOSCA code are shown in which an electromagnetic for a particle accelerator was optimised and designed before, and subsequently verified experimentally after construction. 5.3. free-dimensional
eddy current problems
Provided the basic assumptions for neglecting the displacement current term in (2) are valid, then solutions to eddy currents problems should be attainable by solving the quasi-static subset of Maxwell’s equations-see Section 2 and 3. However, there are many possibilities for selecting a solution variable. 53.1. Total field vector R The total field vector I? is attractive since in conducting regions there are three unknowns at each discretised point whilst in free space when conductivity is zero a single scalar is
Fig. 3(a). Cancellation errors using reduced scalar potential. The test model consisted of a 3D slice of an infinitely long C magnet powered by a pair of coils. The inside and outside radii of the C-core were 5.0 and 8.0cm respectively. Eight-node hexahedra were used, the model being equivalent to a 2D system with 400 nodes. The picture shows iron discretisation only.
C. W. Trowbridge, Low frequency electromagnetic jeld computation
664
.IL8> 0
0
-
Ii
180 Fig. 3(b). Cancellation errors using reduced scalar potential. The results are plotted on a semi-circle of 6.125 cm around the origin, the path intersecting element centroids. The full line shows the results using the total scalar potential which has the correct behaviour (accuracy <5%), the circles are for reduced potential with solution evaluated at element centroids, the crosses for solution values not at element centroids.
sufficient (see (23) and (24)). Furthermore, this approach will result in calculation of the fields directly, i.e., there will be no loss in accuracy caused by the numerical differentiation of potentials; also, possibly of greater importance, there will be no loss of accuracy from field cancellation due to the subtraction of partial fields. However, a difficulty does exist at the interface between conducting regions of differing permeability p. For example, for constant p and u (8) reduces to
v2H=pu$.
(45)
Equation (45) in Cartesian coordinates separates out into equations for each component, thus applying a weighted residual method [ll] to the ‘x’ component equation, the residual at the ith node is given by
and after integrating
the first term by parts, (46)
C. W. Trowbridge, Low frequency electromagnetic field computation
665
Fig. 4. SNS accelerator extraction magnet EBV2. Computer model of end region showing finite element mesh used in the TOSCA [14] code.
The surface integral will not cancel automatically in (46) for an element edge in the interface separating two regions, since the coupling of the integrands will not in general satisfy the continuity condition (9). A further problem arises also at any sharp corner where the normal is not defined. These difficulties can be avoided by the introduction of edge variables-see for example [ 191.
5.3.2. T-Llformulation This formulation involves the partitioning fi=
of the total H field into two parts, i.e.,
T-gradR,
where F is called the electric vector potential from (14),
(47) and L! is the magnetic
curl rT = .7 ; hence F has the same ‘curl’ as fi but need not have the same divergence.
scalar potential.
Thus
(48) In fact from (1) (5)
666
C. W. Trowbridge, Low frequency electromagnetic field computation
--_-FLL-.-L._ 9o EVB2
Colrulatict~~
f. Meosurement
at 1835.53A
X
Fig. 5. SNS accelerator extraction magnet EBV2 results. Contours of differences between calculated and measured fields.
(6), (7) and (8) the defining equation
in p is given by
curl curl F = -pus (aF-Vfi), x
(49)
div F = V20
(501
if fi, 0 are constant. Extensive work, using this method, has been carried out by the group at GEC Stafford [20]; however, there is a problem at the interface where the normal component of T is discontinuous and an extra term has to be added to the equations to satisfy (9) but, unlike the H formulation, the conducting regions can be readily coupled to the magnetic scalar potential.
C. W. Trowbridge, Low frequency electromagnetic $eld computation
667
Care must also be taken to minimise the cancellation error introduced by (47) itself. The question of uniqueness of T is also unclear since the Helmholtz theorem requires both the ‘curl’ and the ‘divergence’ to be defined, whereas here the divergence condition is only indirectly defined. However unique solutions are readily obtained in practice. 5.3.3. R-lj%formulution Since the T vector in the F-0 method is not continuous at the interface an alternative closely related approach [21] is to define a vector R such that
in the region L?,, where the conductivity Here, from (1) and (14),
is non-zero,
and J/ is a magnetic
but
scalar potential.
curl curl R = up $ (R - Vr(r)
(52)
and also, from (4) and (51) divp(E-V$)=O. Now, if in the surrounding R=: -VJ/ with
(53) volume
JzZ which contains
0,
the conductivity
VjXV$=O. In order to couple the two regions together,
is zero, then
(54) (9) and (10) must be satisfied, i.e.,
(56) Now, if in (55) and (56) it is specified that t(r and ~(~~/~~) are each continuous, then it follows that Z?t and @,, are zero on all r,,. Thus l? and $J satisfy (52) and (53) in R; $ satisfies (54) in J& and R’ = 0 on fi2. This approach offers the simplest formulation in non-conducting regions and appears as a natural extension of the two-scalar approach used in the statics case (see Section 5.2). However, as these equations stand the functions R and Yr,are not unique. Many different solutions that satisfy (52) and (53) could be found by adding part of the gradient of J/ to R. However, it was predicted that the finite element discretisation would force a unique solution [21] and reasonable results were obtained. The advantage of this technique is that the difficult interface problem disappears completely although the cancellations problem in the interior of a conducting region still remains.
C. W. Trowbridge, LQWfrequency electromagnetic field computation
668
5.3.4. A-4 formulation
The vector potential A defined by (18) has already been introduced, see Section 3, and has been used very effectively for two-dimensional problems. The extension to 3D appears, however, to involve a further unknown electric potential inside conducting regions, see (20), and will result in 4 unknowns per mesh point. Nevertheless application of the weighted residual method to (20) will result in the correct interface condition between regions of different permeability. Suppose, for example, that a Galerkin weighted residual method is applied to (20) i.e., (57)
where w is a suitable vector weight function. In order to ensure C, continuity across element boundaries it is necessary to apply the vector form of Green’s theorem to (57) [2], i.e.,
and it can be seen that the surface term in (58) is just the tangential component of fi which must be continuous by-virtue of (10). The divergence of A has not yet been specified, although by taking the divergence of both sides of (20) implies divrr($+Vv)=O, cf. (21). Recently a successful attempt to eliminate the electrical potential in conducting regions has been reported [22]; this is achieved by introducing a modified vector potential A* that satisfies curl A* = B and divaA*=O;
60)
hence from (59)
aA* -= at
e+vv,
(61)
at
and (20) simplifies to 1 curl -curl P
aA* A* = -u ~ at
*
(62)
This is the same result that would have been obtained by simply neglecting the electric scalar
C. W. Trowbridge, Low frequency electromagnetic field computation
669
V from the outset as has often been done in the past. However, here this procedure is justified by (68), and furthermore, the modified vector A* has been shown to be unique for the time harmonic case in [22]. Uniqueness can be shown generally as follows: From the vector form of Green’s theorem and by (62)
I(
L[curl~*i2+~~*@&)dR=~ fi P
(A*xVxA*).fidr. r
(63)
Now if A* is not unique there will be at least two solutions AT and A: and accordingly the difference
satisfies (62) and (63). At the boundary AT = A*2, and hence both SA* and curl A* = 0 at the boundary. Furthermore, at time t = 0, A* = 0 everywhere, and so by integrating (63) with respect to time,
i.e.,
The left-hand side of (64) is positive for all SA*, therefore it follows that SA* = 0 for all values of t. Thus AT = /iz and A* is unique. Also in [22] it was shown the conducting region can be coupled to a non-conducting region defined by a single scalar potential 4. This formulation meets most of the difficulties; i.e., it has the minimum number of variables, well-behaved interface conditions, and a symmetric system of equations. However, the possible loss of accuracy in obtaining fields from potentials still remain a difficulty although there exists the possibility of improving the field solutions by resolving the problem curl A* = B by a Galerkin method to improve the B values. Furthermore the behaviour of the formulation as g + 0 is problematical. 53.5. Comparative results In this section four eddy current formulations are compared with experimental measurement; namely, (a) ‘T-L! method as developed by GEC Stafford [20], (b) ‘R-3/’ and ‘A-4’ methods as developed by RAL, and (c) a circuit analogue method developed at Bath University [23]. This latter group also conceived, designed and constructed the test rig for the experiments, the so-called ‘Bath Cube’ experiment [23]. See Fig. 6 for a schematic drawing of the experiment. In the experiment the fields at the surface of a conducting block inserted into the gap of an electromagnet producing a 50 Hz field were measured. The results for Case 2 are in Table 1 and show the vertical field and phase as a function of position. In the last two columns the RMS error and number of degrees of freedom used are given. All the details of the results are given in the literature cited above except some new results obtained at RAL for the A-4 formulation which are in Fig. 7 for Case 2.
670
C. W. Dowbridge, Low frequency electromagnetic field computation
Table 1 Comparison (Case 2)
of 3D eddy current formulations
x=0
and methods against measurements
x=3
,X=.5
x=7
x=9
from the ‘Bath Cube’ experiment
x=11
Method
F&is error
Degrees of freedom
1.08
203
77
r,
3 E ‘r: ot w” z Lzl c cl
77 z f
Modulus (mT) Error Phase (“) Error
-10
-80
-60
Modulus (mT) Error Phase (“) Error
15.6 0.6 -8 2
8.8 1.5 -89 -9
8.2 0.4 -73 -13
Modulus (mT) Error Phase (“) Error
14.1 -0.9 -11 -1
6.5 -0.8 -65 15
Modulus (mT) Error Phase (“) Error
13.7 -1.3 -14 -4
Modulus (mT) Error Phase (“) Error
12 -3 -27 -17
F-7 IL?4 Modulus (mT) Error -z Phase (“) $ Error z
15
18.5 3.5 -2 8
7.3
7.8
13.6
12.0
7.0
0
0
14.2 0.6 -6 -4
13 1 -2 -2
6.8 -2 -1 -1
7.1 -0.7 -51 9
12.6 -1 -7 3
11 -1 -1.5 -1.5
6.3 -0.7 -1.0 -1.0
6.7 -0.6 -82 -2
7.1 -0.7 -68 -8
12.1 -1.5 -10 0
11.1 -0.9 -2 -2
6.2 -0.8 -1 -1
7.6 -0.3 -80 0
7.4 -0.4 -75 -15
if.7 -1.9 -14 -4
11.1 -0.8 -0.5 -0.5
5.8 -1.2 -0.6 -0.6
6.4 -0.9 -71 9
6.7 -1.1 -67 -7
13.0 1.0 -0.6 -0.6
6.9 -0.1 -0.4 -0.4
-10
15.9 2.3 -1 9
6.46
0.84
888
6.76
0.94
659
3.57
1.5
242
8.7
1.71
112
6.27
The results for the A-d) formulation as displayed in Fig. 7 indicate some convergence toward the shape of the measured results as the discretisation is refined but it will be noted that there is still a residual error of the order of 10%. At this time there is no information available on mesh refinement experiments for the other approaches. Substantial improvements are to be expected if the mesh is refined locally where the errors are largest and also if the external magnet is modelled more accurately. A general point worth making is that all of the techniques agree roughly, both with each other and with the measurements despite the crudity
C. W. Trowbridge, Low frequency electromagnetic jieti computation
Section
671
(Cases 1 & 2)
-r 1 0.22
nqm
Case 1
:p:
Case 2
Case 3
Plan View of Blocks Fig. 6. ‘Bath Cube’ eddy current experiment.
of the modelling. The real test of these methods will arise when more complex 3D problems are analysed and accordingly future work will be in this direction; it is also hoped to introduce the new adaptive and hierarchial element schemes using ‘a posteriori’ error estimates in order to achieve automatic mesh refinement and economy [24].
6. Conclusions The present technology of computational modelling, for electromagnetic systems, allows good predictions to be made for two- and three-dimensional statics problems, and for two-dimensional eddy current problems. The accuracy achieved will depend upon the degree of discretisation and on the knowledge available for the material properties involved; also the degree of discretisation useable depends upon the computing power available. In practice, today, the quality of the results from computer codes still depends critically upon the expertise available; it is a common experience that users with knowledge of both computational modelling and electromagnetics achieve the best results. Consideration of integral methods for low frequency problems was outside the scope of this paper; however, it is clear from the literature that, for a wide range of problems, the use of the integral operator formulations can be very effective, see [30]. The three-dimensional eddy current problem is still the subject of extensive research; and some preliminary results from ‘work in progress’ in the U.K. are reported in this paper. This work is only a small part of an international community undertaking investigations into methods, see [25-291; and so there is hope that the reliable, robust, efficient algorithms will be available shortly.
672
C. W. Trowbridge, Low freqwncy electromagnetic j&U computation
B, (MT)
0
0
“f-y----;w
14 X (CM)
0
NUMERICAL
2637 DEGREES OF
_
-5
FREEDOM
652 DEGREES OF FREEDOM
-i0
242 DEGREES OF
,,I,
FREEDOM
EXPERIMENTAL
BATH
CUBE -
B, ALONG LINE
Fig. 7(a).Comparisons versus measurement.
of ‘A*-+’
CASE 2
Y = 5 CM,
Z = 0,2
method with measurement
CM
for the ‘Bath Cube’ Experiment.
Moduh~s and phase
The author is indebted to colleagues John Simkin and Cris Emson for information on their recent work on three-dimensional eddy currents and to Mrs. Pam Peisley for preparing the typescript and iHustrations.
673
C. W. Trowbridge, Low frequency electromagnetic field computation
ERROR (MT)
I
x
I
_
-4
(CM)
2637 DEGREES
OF
652
DEGREES FREEDOM
OF
242
DEGREES FREEDOM
OF
FREEDOM
_-_7 ERROR
(DEGREES) ,I,,
x 10
-2 1'
BATH ERRORS
Fig. 7(b). Comparisons of ‘A*+ phase versus measurement.
method
ALONG
CUBE LINE
with measurement
-
CASE
2
Y = 5 CM, z = 0.2 CM
for the ‘Bath Cube’ Experiment.
Errors
in modulus
References [l] [2] [3] [4] [5]
J.C. Maxwell, A dynamical theory of the electromagnetic field, Roy. Sot. Trans. 155 (1864) 459-512. J.A. Stratton, Electromagnetic Theory (McGraw-Hill, New York, 1941) 486. J.A. Stratton, Electromagnetic Theory (McGraw-Hill, New York, 1941) 15. Alan M. Portis, Electromagnetic Fields, Sources and Media (Wiley, New York, 1978) 638. J. Simkin and C.W. Trowbridge, Which potential?, RL-78-001/B, 1978.
and
674
C. W. Trowbridge, Low frequency electromagnetic field computation
[6] A.G. Armstrong, C.J. Collie, J.S. Simkin and C.W. Trowbridge, The solution of three-dimensional magnetostatic problems using scalar potentials, Proc. COMPUMAG Conf., Grenoble (1978) 1.2. [7] A.A. Winslow, Numerical solution of the quasi-linear Poisson equation in a non-uniform triangular mesh, J. Comput. Phys. 1 (1971) 149-172. [S] M.W.K. Chari and P. Silvester, Finite element analysis of magnetically saturated dc machines, IEEE Trans. Power Apparat. Syst. 90 (1971) 2362. [9] S.J. Polak, A. Wachters and A. de Beer, An account of the use of the finite element method for magnetostatics, Proc. COMPUMAG Conf., Oxford, 1976. [lo] C.S. Biddlecombe, N.J. Diserens, C. Riley and J. Simkin, PE2D user guide. Version 6.3, Rutherford Laboratory, Rept. RL-81-089, 1983. [ll] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, Maidenhead, 3rd ed., 1977). [12] O.C. Zienkiewicz, J. Lyness and D.R.J. Owen, Three-dimensional magnetic field determination using a scalar potential-a finite element solution, IEEE Trans. Magnetics 13 (1977) 1649-1656. [13] J. Simkin and C.W. Trowbridge, On the use of the total scalar potential in the numerical solution of field problems in electromagnets, Internat. J. Numer. Meths. Engrg. 14 (1979) 4234!0. [ 141 J. Simkin and C.W. Trowbridge, Three-dimensional non-linear electromagnetic held computations using scalar potentials, Proc. Inst. Elec. Engrg. 127 B(6) (1980). [15] W. Wolff and W. Muller, General numerical solution of the magnetostatic equations, Wiss. Ber. AEGTelefunken 49(3) (1976) 77-86. [16] A. de Beer, S.J. Polak, A.J.H. Wachters and J.S. van Welij, The use of Paddy for the solution of 3-D magnetostatic problems, Proc. COMPUMAG, Chicago 1981; IEEE Trans. Magnetics 18(2) (1982) 617. [17] T.W. McDaniel, R.B. Fernandez, R.R. Root and R.B. Anderson, An accurate scalar potential finite element method for linear, two-dimensional magnetostatic problems, Internat. J. Numer. Meths. Engrg. 19 (1983) 725-737. [18] F.E. Baker, S.H. Brown, J.R. Brauer and T.R. Gerhardt, Comparison of magnetic fields computed by finite element and classical series methods, Internat. J. Numer. Meths. Engrg. 19 (1983) 271-280. [19] A. Bossavit and J.C. Verne, A fixed FEM-BIEM method to solve 3-D eddy current problems, Proc. COMPUMAG, Chicago, IL, 1981; IEEE Trans. Magnetics 18(2) (1982) 431. [20] T.W. Preston and A.B.J. Reece, Solution of three-dimensional eddy current problems: the T-L! method, Proc. COMPUMAG Conf., Chicago, IL, 1981; IEEE Trans. Magnetics 18 (2) (1982) 486. [21] C.S. Biddlecombe, E.A. Heighway, J. Simkin and C.W. Trowbridge, Methods for eddy current computation in three dimensions, Proc. COMPUMAG Conf., Chicago, IL, 1981; IEEE Trans. Magnetics 18(2) (1982) 492. 1221 C.R.I. Emson and J. Simkin, An optimal method for 3-D eddy currents, Proc. COMPUMAG Conf., Genoa, 1983; IEEE Trans. Magnetics 19(6) (1983). [23] J.J. Balchin and J.A.M. Davidson, Numerical method for calculating magnetic flux and eddy current distribution in three dimensions, IEE Proc. A 127(l) (1980) 6.53. [24] D.W. Kelly, S.R. Gago, O.C. Zienkiewicz and I. Babushka, A posteriori error analysis and adaptive processes in the finite element method: Parts 1 and 2, Intemat. J. Numer. Meths. Engrg. 19 (1983) 1593. [25] J. Simkin, ed., COMPUMAG Conference Proceedings (Rutherford Appleton Laboratory, Oxford, 1976). [26] J.C. Sabonnadiere, ed., COMPUMAG Conference Proceedings (Laboratoire d’Electrotechnique de Grenoble, ERA 524 CNRS, Grenoble, 1979). [27] COMPUMAG Conference Proceedings, Chicago, IL, 1981, IEEE Trans. Magnetics 18(2) (1982). [28] COMPUMAG Conference Proceedings, Genoa, 1983, IEEE Trans. Magnetics 19(6) (1983). [29] COMPUMAG Conference Proceedings, Fort Collins, CO, 1985. [30] C.W. Trowbridge, Applications of integral equation methods for the numerical solution of magnetostatic and eddy current problems, Paper presented to Internat. Conf. Numerical Methods in Electrical and Magnetic Field Problems, Santa Margherita, Italy, 1976; see also Rutherford Laboratory Rept. RL-76-071, 1976, and Finite Elements in Electrical and Magnetic Field Problems (Wiley, Chichester, 1979) 191-213.
IN APPLIED
MECHANICS
AND ENGINEERING
52 (1985) 653-674
LOW FREQUENCY ELECTROMAGNETIC FIELD COMPUTATION IN THREE D~NSIONS C.W. TROWBRIDGE Ruthe&rd
Appleton Laboratory, Chilton, Didcot, Oxon, U.K. Received 22 October
1984
This paper presents the various differential operator formulations for obtaining numerica solutions to the classical equations of electromagnetics using the finite element method. Numerical examples are presented for 2D static and eddy currents problems. Three-dimensional formulations are discussed for nonlinear static problems using the total and reduced scalar potentials with results from a recent magnet project involving field computation at the design stage and validation after const~ction. Finally a review of some of the formulations for solving 3D eddy current problems are discussed with comparative results against measurement.
1. Introduction Computational modelling is now a basic requirement for the engineer who wishes to design electromagnetic devices. The potential list of end-users of electromagnetic applicable software is enormous, including engineers working on the following applications: equipment and facilities used in physics and engineering requiring magnetic fields of precise shape, for example, charged particle focussing magnets for nuclear physics accelerators, electron beam lithography, and plasma physics containment systems; electromagnetic systems for medical diagnostics such as magnets for nuclear magnetic resonance (NMR) body scanning; probes and evaluation of defects in non-destructive testing equipment; and a multitude of electrical machines and power handling devices of all kinds. A challenging list of problems. In all these cases the designer starts with a performance specification; he then creates from his experience and skill a model which he hopes will go some way toward matching the given specification. In order to test his ideas he needs to predict the performance of his design; first by relatively simple analytic and numerical procedures, and secondly by sophisticated computer techniques, taking account of as many of the design parameters as possible. For many applications two-dimensional simulations are sufficient, but there is an increasing demand for complete three-dimensional solutions in order to minimise the need for building expensive prototypes. Therefore, at the heart of a ‘computer-aided design’ simulation system for an electromagnetic device an algorithm for solving the electromagnetic field equations is required. 00457825/85/$3.30
@ 1985, Elsevier Science Publishers B.V. (North-Holland)
C. W. Trowbridge, Low frequency electromagnetic $eld computation
654
2. Equation for electromagnetic
fields
The classical field equations of Maxwell [l] are sufficient for the range of electromagnetic applications under consideration; these equations together with the nonlinear constitutive laws and appropriate boundary conditions define the field uniquely [2]. The canonical form of these equations are usually written as follows:
(1) (2) div D = p ,
(3)
(compatibility) div I? = 0 ,
(4)
where the vectors J? and H are the electric and magnetic field intensities and the vectors fi and B are the electric and magnetic flux densities respectively. The electric current density and electric charge density are .? and p respectively. The field intensity and flux density quantities are not independent; they are related by the constitutive laws of the material,
B=pH,
(5) (6)
D=&,
where p and E are the permeability and susceptibility respectively. In general both p and E are tensors and ,!Lin particular will be field-dependent. The electric field intensity vector g is further related to the current density J by Ohm’s law,
tTE=J,
(7)
where (T is the conductivity of the material. The above equations can be combined to produce many different forms by eliminating one or more of the dependent fields; for example if p, E and CTare constant, a single equation results for each of the field quantities, i.e. for the magnetic intensity R,
(8) At interfaces between regions with different material constants following continuity conditions: ii*(B2-&)=O, fix(tl;!-HI)=
the field equations
imply the
(9) K,
(IO)
C. W. Trowbridge, Low frequency electromagnetic field computation
A@*-D1)=
q,
655
01)
iix(&BJ=O,
(12)
where ri is the positive normal of the surface separating the two regions (1) and (2), and where K and q are possible surface current and charge densities respectively. In conducting regions (a % .F) for slowly varying time fields the second term on the right-hand side of (8) is negligible; this is the case for time-harmonic problems in which R = H0 exp(-iwt)
(13)
and where o, the angular frequency, is small, such that the wavelength is large compared with the dimensions of the problem. This is equivalent to neglecting the aD/& term in (2), and the resulting subset of equations is known as the quasi-static Maxwell equations, i.e. (2) becomes curl fi = .T.
(14)
Hence the displacement current (Jo/at) term in Maxwell’s equations will be neglected in all that follows since only static and low (power) frequency problems will be considered. At the surface of a conducting region the normal component of current density J, will be assumed to be zero since any tendency for current to flow is opposed by the electric charges accumulating at the surface. In fact, taking the divergence of (2) and using (4) the equation of continuity results, divj=
--@fat.
It is shown in the standard texts [3] that (15) implies an short relaxation times for conductors; even for a poor microsecond. In the low frequency limit approximation instantaneously to the conductor surface; so for interior
(1%
exponential decay of charge with very conductor like water it is less than a it is assumed that free charges migrate points
divj=O,
(16)
and on the surfaces
&=o.
(17)
3. Magnetic vector potentials Equation
(4) implies the existence of a vector potential
A,
B=curlA, which by Helmholtz’s theorem [4] will be unique if the divergence is specified together with the tangential component on the boundary. Since in the time-dependent case both electric and
C. W. Trowbridge, Low frequency electromagnetic field computation
656
magnetic fields are coupled, (1) and (18) imply E=-($+W),
(19)
where V is an electric scalar potential; equation in A-V results,
hence by combining
(l), (2) and (5) an equilibrium
curlbcurlA=-g($+DV);
(20)
also, from (16)
(21) Equations (20) and (21) are the defining equations for a quasi-static conductivity u and permeability ,s. It will be sufficiently general embedded inside a region of free space L& which may be infinite. The defining equation for region 0, is simpler, since u = 0,
field in a region 0, of to consider 0, to be
curlLcurlA=O; P
(22) is possible since curl H = 0 and a magnetic
indeed, further simplification may be introduced,
scalar potential
B = -grad 9,
$
(23)
which by (4) satisfies
v2+ = 0.
(24)
On the interface r12 between 0, and 0, the standard continuity must be satisfied, @“)I
=
@“)2,
(ml
=
(fix
conditions,
i.e., (9) and (lo),
(25)
*
Thus the scalar $ must be coupled to the vector A in such a way as to satisfy (25). Also at the surface f,, (17) must be satisfied, which from (19) implies
aA” av
u(
-+-& at
Hence, it appears
>
=o.
that (20) and (21) for region a,,
(26) and (24) for region J2,, together
with
C. W. Trowbridge, Low frequency electromagnetic field computation
657
interface and coupling conditions, (25) and (26) and the external boundary conditions on region a,, completely define the problem. The choice of divergence, or gauge, on A will be considered later, see Section 5.2.4.
4. Static limiting case There are three limiting cases to consider - electrostatics, The electrostatic case is straightforward, i.e. from (1)
magnetostatics
and current flow.
curl E = 0, which implies the existence of a scalar potential
V,
E = -vv, and from the divergence
(27) condition
(3),
(28)
v*Evv=p.
Hence, electrostatic problems will satisfy Poisson’s equation (28) and the electrostatic potential V will be unique in a domain L! provided V (Dirichlet condition) or aV/dn (Neumann condition) are prescribed on the boundary-the solution of a pure Neumann problem (cYV/&I defined on the boundary only) will be unique if V is defined at one point at least. For magnetostatic problems the conductors are usually prescribed, carrying a current density x. Hence the vector potential A could be used; i.e., from (14) and (18) curl L curl A = Js .
(2%
P
Equation (29) is most useful in two dimensions since in this case A has only a single component; in three dimensions there is the problem of gauge to contend with as well. Furthermore, an A formulation would require the conductors to be discretised - a formidable problem in three dimensions. It is more economical to explore the possibilities of scalar potentials, of which there are two kinds. The first kind is usually termed the total scalarpotential and is only valid for the conductor-free regions of the problem, i.e., curl fi = 0 ;
(30)
fi=
(31)
then -V$,
where $ is the total scalar potential. two parts [5],
The second kind arises by partitioning
the field H into
(32)
C. W. Trowbridge, L.ow frequency electromagnetic fit&i computation
658
where gS is defined by (33) which implies that
and (34) in which 4 is often called the reduced scalar potential (4) Cc,satisfies
[6] and is valid everywhere
v-pvt#J=o
[5]. From
(35)
and C$satisfies V.pV#
= div,uHS,
(36)
with suitable boundary conditions. Finally, for problems of electrical conduction (16) and (7)
the definition of a scalar potential follows from
divJ=Vcr-B=V.aVV=O.
(37)
Thus the potentials for all three types of problem therefore all use the same algorithms.
satisfy the quasi-Poisson
equation
and can
are straightforward.
In this
5. Numerical solutions 5.1. Two-dimensional problems Numerical solutions of the field equations for two dimensions case (20), for the single component vector potential, reduces to
and finite difference and finite element methods have been extensively used to solve (38)-the literature is very rich in examples [7,8,9]. Two examples of two-dimensional solutions are shown in Figs. 1 and 2, in which the package PE2D (Poissons Equation in 2 Dimensions) [lo] has been used to solve a plane magnetostatic problem
and an eddy current
problem
with axial symmetry.
C. W. Trowbridge, Low frequency electromagnetic jield computation
659
3t
I-
I-
I
-6.OOC 3
Fig. l(a). Permanent magnet motor. Finite element mesh generated processor [lo]. There are 2508 first-order triangles in the mesh.
5.2. Three-dimensional
by using the PE2D computer
code pre-
static problems
It is clear that solutions for electrostatic and current flow problems can readily be obtained by the use of scalar potentials and standard finite element methods- there are numerous examples in the literature [ll]. The magnetostatics case is not so clear; e.g. in Section 4 two kinds of scalar potential were defined but it can easily be seen that the presence of conductors carrying electric current in a problem poses difficulties in using a total scalar potential since it is not single-valued. The reduced potential which avoids this difficulty, has been used with some success, especially for predicting fields in the free space regions 16,121. The use of the reduced potential was investigated at Rutherford Appleton Laboratory [13] with results that indicated that, under many circumstances, erroneous fields could be computed in highly permeable regions. The reasons for this were identified with the defining
C. W. Trowbridge, Low frequency electromagnetic field computation
660 Cl.OO(lr
4.001
2
~ 001
G.OOl
I
3
!
3
I
1
-2.i)Ol
-4.OOI
-i-1
~ clril
_ i -5.000
I
0.000
Fig. l(b). Permanent magnet motor. Finite element solution displayed by contours of constant This example is the subject of a patent held by Prof. K. Binns of Liverpool University.
5.000 vector potential.
equation itself, i.e. (32), which may be expressed as
in which it is obvious that if both E?m and V4 are of similar magnitude, VC#Jwill have to be computed to the same order of accuracy as r?, if cancellation errors are to be avoided. In fact, E;i, is often calculated by analytic or by numerical quadrature to high accuracy and can have a strong spatial variation. Furthermore, errors can be exacerbated when 8 is computed by (5) in a nonlinear iteration for predicting saturation effects in iron. For these reasons it makes good sense to calculate the difference quantity direct by use of the total scalar potential. This can easily be achieved for a very large range of problems by subdividing the domain of the problem into subdomains and selecting the total potential for those regions that do not contain conductors and the reduced potential for the remainder.
C. W. Trowbridge, Low frequency electromagnetic field computation
4 II
II
II 0 0
661
i
2ll.o00-
10.000-
L \
-I
Fig. Z(a). Axi-symmetric eddy currents. The picture shows contours of induced electric current density in a molten metal experiment for investigating forces in an induction furnace. The experiment consists of mercury in a stainless steel pot surrounded by a 50 Hz coil. The problem was solved using PE2D code [lo] with a mesh consisting of 1813 second-order triangles.
Thus, for a two-region
problem the defining equations
are as follows: for region 1,
for region 2,
There is a saving in cost as well because the source field E?, does not appear in (40); in fact, when the two regions are coupled together the source term reappears for points ‘on the interface. The two regions are coupled by use of the standard interface conditions, i.e. (25), which when re-expressed in terms of potentials are
(41)
C. W. Trowbridge, Low frequency electromagnetic field computation
662
4o.oot
-
i 3o.ooc
-_ -
20.00[
-
1o.ooc
-
o.oot
-t
Fig. 2(b). Axi-symmetric
- i
eddy currents.
Lines of constant
(rA,) giving the rest&ant
($),=(iG-$),.
by quadrature
flux.
(42)
The explicit relation between 4 and 1,5is obtained by integrating interface surface r12, i.e.,
With H, determined
magnetic
from the Biot-Savart
law,
(42) along a contour over the
C. W. Trowbridge, Low frequency electromagnetic field computation
663
Permanent magnet sources can be taken into account by adding the coercive field to (32) and by using the appropriate constitutive relationship [14). Again a very efficient ~go~thm can be obtained by applying the Galerkin method to (39) and (40) and coupling the two regions together by (43) [13]. A general code for threedimensional static fields called TOSCA (TWO SCAlar potentials) has been developed at Rutherford Appleton Laboratory, based on this approach, with several comparisons made against measurement [ 141. In fact the literature of viable 3D computer implementations even for the statics case is not very extensive, see [14,15,16] for a selection. The practice of adopting successful computer codes from other ‘fields’ to electro-magnetics has been reported from time to time [17,18] this should be done with great care and with regard to the cancellation errors arising from the exclusive use of the reduced scalar potential-see (32) and comments above. That serious errors arise can be seen in Fig. 3 for results of an experiment carried out at RAL. In Figs. 4 and 5 successful 3D results using the TOSCA code are shown in which an electromagnetic for a particle accelerator was optimised and designed before, and subsequently verified experimentally after construction. 5.3. free-dimensional
eddy current problems
Provided the basic assumptions for neglecting the displacement current term in (2) are valid, then solutions to eddy currents problems should be attainable by solving the quasi-static subset of Maxwell’s equations-see Section 2 and 3. However, there are many possibilities for selecting a solution variable. 53.1. Total field vector R The total field vector I? is attractive since in conducting regions there are three unknowns at each discretised point whilst in free space when conductivity is zero a single scalar is
Fig. 3(a). Cancellation errors using reduced scalar potential. The test model consisted of a 3D slice of an infinitely long C magnet powered by a pair of coils. The inside and outside radii of the C-core were 5.0 and 8.0cm respectively. Eight-node hexahedra were used, the model being equivalent to a 2D system with 400 nodes. The picture shows iron discretisation only.
C. W. Trowbridge, Low frequency electromagnetic jeld computation
664
.IL8> 0
0
-
Ii
180 Fig. 3(b). Cancellation errors using reduced scalar potential. The results are plotted on a semi-circle of 6.125 cm around the origin, the path intersecting element centroids. The full line shows the results using the total scalar potential which has the correct behaviour (accuracy <5%), the circles are for reduced potential with solution evaluated at element centroids, the crosses for solution values not at element centroids.
sufficient (see (23) and (24)). Furthermore, this approach will result in calculation of the fields directly, i.e., there will be no loss in accuracy caused by the numerical differentiation of potentials; also, possibly of greater importance, there will be no loss of accuracy from field cancellation due to the subtraction of partial fields. However, a difficulty does exist at the interface between conducting regions of differing permeability p. For example, for constant p and u (8) reduces to
v2H=pu$.
(45)
Equation (45) in Cartesian coordinates separates out into equations for each component, thus applying a weighted residual method [ll] to the ‘x’ component equation, the residual at the ith node is given by
and after integrating
the first term by parts, (46)
C. W. Trowbridge, Low frequency electromagnetic field computation
665
Fig. 4. SNS accelerator extraction magnet EBV2. Computer model of end region showing finite element mesh used in the TOSCA [14] code.
The surface integral will not cancel automatically in (46) for an element edge in the interface separating two regions, since the coupling of the integrands will not in general satisfy the continuity condition (9). A further problem arises also at any sharp corner where the normal is not defined. These difficulties can be avoided by the introduction of edge variables-see for example [ 191.
5.3.2. T-Llformulation This formulation involves the partitioning fi=
of the total H field into two parts, i.e.,
T-gradR,
where F is called the electric vector potential from (14),
(47) and L! is the magnetic
curl rT = .7 ; hence F has the same ‘curl’ as fi but need not have the same divergence.
scalar potential.
Thus
(48) In fact from (1) (5)
666
C. W. Trowbridge, Low frequency electromagnetic field computation
--_-FLL-.-L._ 9o EVB2
Colrulatict~~
f. Meosurement
at 1835.53A
X
Fig. 5. SNS accelerator extraction magnet EBV2 results. Contours of differences between calculated and measured fields.
(6), (7) and (8) the defining equation
in p is given by
curl curl F = -pus (aF-Vfi), x
(49)
div F = V20
(501
if fi, 0 are constant. Extensive work, using this method, has been carried out by the group at GEC Stafford [20]; however, there is a problem at the interface where the normal component of T is discontinuous and an extra term has to be added to the equations to satisfy (9) but, unlike the H formulation, the conducting regions can be readily coupled to the magnetic scalar potential.
C. W. Trowbridge, Low frequency electromagnetic $eld computation
667
Care must also be taken to minimise the cancellation error introduced by (47) itself. The question of uniqueness of T is also unclear since the Helmholtz theorem requires both the ‘curl’ and the ‘divergence’ to be defined, whereas here the divergence condition is only indirectly defined. However unique solutions are readily obtained in practice. 5.3.3. R-lj%formulution Since the T vector in the F-0 method is not continuous at the interface an alternative closely related approach [21] is to define a vector R such that
in the region L?,, where the conductivity Here, from (1) and (14),
is non-zero,
and J/ is a magnetic
but
scalar potential.
curl curl R = up $ (R - Vr(r)
(52)
and also, from (4) and (51) divp(E-V$)=O. Now, if in the surrounding R=: -VJ/ with
(53) volume
JzZ which contains
0,
the conductivity
VjXV$=O. In order to couple the two regions together,
is zero, then
(54) (9) and (10) must be satisfied, i.e.,
(56) Now, if in (55) and (56) it is specified that t(r and ~(~~/~~) are each continuous, then it follows that Z?t and @,, are zero on all r,,. Thus l? and $J satisfy (52) and (53) in R; $ satisfies (54) in J& and R’ = 0 on fi2. This approach offers the simplest formulation in non-conducting regions and appears as a natural extension of the two-scalar approach used in the statics case (see Section 5.2). However, as these equations stand the functions R and Yr,are not unique. Many different solutions that satisfy (52) and (53) could be found by adding part of the gradient of J/ to R. However, it was predicted that the finite element discretisation would force a unique solution [21] and reasonable results were obtained. The advantage of this technique is that the difficult interface problem disappears completely although the cancellations problem in the interior of a conducting region still remains.
C. W. Trowbridge, LQWfrequency electromagnetic field computation
668
5.3.4. A-4 formulation
The vector potential A defined by (18) has already been introduced, see Section 3, and has been used very effectively for two-dimensional problems. The extension to 3D appears, however, to involve a further unknown electric potential inside conducting regions, see (20), and will result in 4 unknowns per mesh point. Nevertheless application of the weighted residual method to (20) will result in the correct interface condition between regions of different permeability. Suppose, for example, that a Galerkin weighted residual method is applied to (20) i.e., (57)
where w is a suitable vector weight function. In order to ensure C, continuity across element boundaries it is necessary to apply the vector form of Green’s theorem to (57) [2], i.e.,
and it can be seen that the surface term in (58) is just the tangential component of fi which must be continuous by-virtue of (10). The divergence of A has not yet been specified, although by taking the divergence of both sides of (20) implies divrr($+Vv)=O, cf. (21). Recently a successful attempt to eliminate the electrical potential in conducting regions has been reported [22]; this is achieved by introducing a modified vector potential A* that satisfies curl A* = B and divaA*=O;
60)
hence from (59)
aA* -= at
e+vv,
(61)
at
and (20) simplifies to 1 curl -curl P
aA* A* = -u ~ at
*
(62)
This is the same result that would have been obtained by simply neglecting the electric scalar
C. W. Trowbridge, Low frequency electromagnetic field computation
669
V from the outset as has often been done in the past. However, here this procedure is justified by (68), and furthermore, the modified vector A* has been shown to be unique for the time harmonic case in [22]. Uniqueness can be shown generally as follows: From the vector form of Green’s theorem and by (62)
I(
L[curl~*i2+~~*@&)dR=~ fi P
(A*xVxA*).fidr. r
(63)
Now if A* is not unique there will be at least two solutions AT and A: and accordingly the difference
satisfies (62) and (63). At the boundary AT = A*2, and hence both SA* and curl A* = 0 at the boundary. Furthermore, at time t = 0, A* = 0 everywhere, and so by integrating (63) with respect to time,
i.e.,
The left-hand side of (64) is positive for all SA*, therefore it follows that SA* = 0 for all values of t. Thus AT = /iz and A* is unique. Also in [22] it was shown the conducting region can be coupled to a non-conducting region defined by a single scalar potential 4. This formulation meets most of the difficulties; i.e., it has the minimum number of variables, well-behaved interface conditions, and a symmetric system of equations. However, the possible loss of accuracy in obtaining fields from potentials still remain a difficulty although there exists the possibility of improving the field solutions by resolving the problem curl A* = B by a Galerkin method to improve the B values. Furthermore the behaviour of the formulation as g + 0 is problematical. 53.5. Comparative results In this section four eddy current formulations are compared with experimental measurement; namely, (a) ‘T-L! method as developed by GEC Stafford [20], (b) ‘R-3/’ and ‘A-4’ methods as developed by RAL, and (c) a circuit analogue method developed at Bath University [23]. This latter group also conceived, designed and constructed the test rig for the experiments, the so-called ‘Bath Cube’ experiment [23]. See Fig. 6 for a schematic drawing of the experiment. In the experiment the fields at the surface of a conducting block inserted into the gap of an electromagnet producing a 50 Hz field were measured. The results for Case 2 are in Table 1 and show the vertical field and phase as a function of position. In the last two columns the RMS error and number of degrees of freedom used are given. All the details of the results are given in the literature cited above except some new results obtained at RAL for the A-4 formulation which are in Fig. 7 for Case 2.
670
C. W. Dowbridge, Low frequency electromagnetic field computation
Table 1 Comparison (Case 2)
of 3D eddy current formulations
x=0
and methods against measurements
x=3
,X=.5
x=7
x=9
from the ‘Bath Cube’ experiment
x=11
Method
F&is error
Degrees of freedom
1.08
203
77
r,
3 E ‘r: ot w” z Lzl c cl
77 z f
Modulus (mT) Error Phase (“) Error
-10
-80
-60
Modulus (mT) Error Phase (“) Error
15.6 0.6 -8 2
8.8 1.5 -89 -9
8.2 0.4 -73 -13
Modulus (mT) Error Phase (“) Error
14.1 -0.9 -11 -1
6.5 -0.8 -65 15
Modulus (mT) Error Phase (“) Error
13.7 -1.3 -14 -4
Modulus (mT) Error Phase (“) Error
12 -3 -27 -17
F-7 IL?4 Modulus (mT) Error -z Phase (“) $ Error z
15
18.5 3.5 -2 8
7.3
7.8
13.6
12.0
7.0
0
0
14.2 0.6 -6 -4
13 1 -2 -2
6.8 -2 -1 -1
7.1 -0.7 -51 9
12.6 -1 -7 3
11 -1 -1.5 -1.5
6.3 -0.7 -1.0 -1.0
6.7 -0.6 -82 -2
7.1 -0.7 -68 -8
12.1 -1.5 -10 0
11.1 -0.9 -2 -2
6.2 -0.8 -1 -1
7.6 -0.3 -80 0
7.4 -0.4 -75 -15
if.7 -1.9 -14 -4
11.1 -0.8 -0.5 -0.5
5.8 -1.2 -0.6 -0.6
6.4 -0.9 -71 9
6.7 -1.1 -67 -7
13.0 1.0 -0.6 -0.6
6.9 -0.1 -0.4 -0.4
-10
15.9 2.3 -1 9
6.46
0.84
888
6.76
0.94
659
3.57
1.5
242
8.7
1.71
112
6.27
The results for the A-d) formulation as displayed in Fig. 7 indicate some convergence toward the shape of the measured results as the discretisation is refined but it will be noted that there is still a residual error of the order of 10%. At this time there is no information available on mesh refinement experiments for the other approaches. Substantial improvements are to be expected if the mesh is refined locally where the errors are largest and also if the external magnet is modelled more accurately. A general point worth making is that all of the techniques agree roughly, both with each other and with the measurements despite the crudity
C. W. Trowbridge, Low frequency electromagnetic jieti computation
Section
671
(Cases 1 & 2)
-r 1 0.22
nqm
Case 1
:p:
Case 2
Case 3
Plan View of Blocks Fig. 6. ‘Bath Cube’ eddy current experiment.
of the modelling. The real test of these methods will arise when more complex 3D problems are analysed and accordingly future work will be in this direction; it is also hoped to introduce the new adaptive and hierarchial element schemes using ‘a posteriori’ error estimates in order to achieve automatic mesh refinement and economy [24].
6. Conclusions The present technology of computational modelling, for electromagnetic systems, allows good predictions to be made for two- and three-dimensional statics problems, and for two-dimensional eddy current problems. The accuracy achieved will depend upon the degree of discretisation and on the knowledge available for the material properties involved; also the degree of discretisation useable depends upon the computing power available. In practice, today, the quality of the results from computer codes still depends critically upon the expertise available; it is a common experience that users with knowledge of both computational modelling and electromagnetics achieve the best results. Consideration of integral methods for low frequency problems was outside the scope of this paper; however, it is clear from the literature that, for a wide range of problems, the use of the integral operator formulations can be very effective, see [30]. The three-dimensional eddy current problem is still the subject of extensive research; and some preliminary results from ‘work in progress’ in the U.K. are reported in this paper. This work is only a small part of an international community undertaking investigations into methods, see [25-291; and so there is hope that the reliable, robust, efficient algorithms will be available shortly.
672
C. W. Trowbridge, Low freqwncy electromagnetic j&U computation
B, (MT)
0
0
“f-y----;w
14 X (CM)
0
NUMERICAL
2637 DEGREES OF
_
-5
FREEDOM
652 DEGREES OF FREEDOM
-i0
242 DEGREES OF
,,I,
FREEDOM
EXPERIMENTAL
BATH
CUBE -
B, ALONG LINE
Fig. 7(a).Comparisons versus measurement.
of ‘A*-+’
CASE 2
Y = 5 CM,
Z = 0,2
method with measurement
CM
for the ‘Bath Cube’ Experiment.
Moduh~s and phase
The author is indebted to colleagues John Simkin and Cris Emson for information on their recent work on three-dimensional eddy currents and to Mrs. Pam Peisley for preparing the typescript and iHustrations.
673
C. W. Trowbridge, Low frequency electromagnetic field computation
ERROR (MT)
I
x
I
_
-4
(CM)
2637 DEGREES
OF
652
DEGREES FREEDOM
OF
242
DEGREES FREEDOM
OF
FREEDOM
_-_7 ERROR
(DEGREES) ,I,,
x 10
-2 1'
BATH ERRORS
Fig. 7(b). Comparisons of ‘A*+ phase versus measurement.
method
ALONG
CUBE LINE
with measurement
-
CASE
2
Y = 5 CM, z = 0.2 CM
for the ‘Bath Cube’ Experiment.
Errors
in modulus
References [l] [2] [3] [4] [5]
J.C. Maxwell, A dynamical theory of the electromagnetic field, Roy. Sot. Trans. 155 (1864) 459-512. J.A. Stratton, Electromagnetic Theory (McGraw-Hill, New York, 1941) 486. J.A. Stratton, Electromagnetic Theory (McGraw-Hill, New York, 1941) 15. Alan M. Portis, Electromagnetic Fields, Sources and Media (Wiley, New York, 1978) 638. J. Simkin and C.W. Trowbridge, Which potential?, RL-78-001/B, 1978.
and
674
C. W. Trowbridge, Low frequency electromagnetic field computation
[6] A.G. Armstrong, C.J. Collie, J.S. Simkin and C.W. Trowbridge, The solution of three-dimensional magnetostatic problems using scalar potentials, Proc. COMPUMAG Conf., Grenoble (1978) 1.2. [7] A.A. Winslow, Numerical solution of the quasi-linear Poisson equation in a non-uniform triangular mesh, J. Comput. Phys. 1 (1971) 149-172. [S] M.W.K. Chari and P. Silvester, Finite element analysis of magnetically saturated dc machines, IEEE Trans. Power Apparat. Syst. 90 (1971) 2362. [9] S.J. Polak, A. Wachters and A. de Beer, An account of the use of the finite element method for magnetostatics, Proc. COMPUMAG Conf., Oxford, 1976. [lo] C.S. Biddlecombe, N.J. Diserens, C. Riley and J. Simkin, PE2D user guide. Version 6.3, Rutherford Laboratory, Rept. RL-81-089, 1983. [ll] O.C. Zienkiewicz, The Finite Element Method (McGraw-Hill, Maidenhead, 3rd ed., 1977). [12] O.C. Zienkiewicz, J. Lyness and D.R.J. Owen, Three-dimensional magnetic field determination using a scalar potential-a finite element solution, IEEE Trans. Magnetics 13 (1977) 1649-1656. [13] J. Simkin and C.W. Trowbridge, On the use of the total scalar potential in the numerical solution of field problems in electromagnets, Internat. J. Numer. Meths. Engrg. 14 (1979) 4234!0. [ 141 J. Simkin and C.W. Trowbridge, Three-dimensional non-linear electromagnetic held computations using scalar potentials, Proc. Inst. Elec. Engrg. 127 B(6) (1980). [15] W. Wolff and W. Muller, General numerical solution of the magnetostatic equations, Wiss. Ber. AEGTelefunken 49(3) (1976) 77-86. [16] A. de Beer, S.J. Polak, A.J.H. Wachters and J.S. van Welij, The use of Paddy for the solution of 3-D magnetostatic problems, Proc. COMPUMAG, Chicago 1981; IEEE Trans. Magnetics 18(2) (1982) 617. [17] T.W. McDaniel, R.B. Fernandez, R.R. Root and R.B. Anderson, An accurate scalar potential finite element method for linear, two-dimensional magnetostatic problems, Internat. J. Numer. Meths. Engrg. 19 (1983) 725-737. [18] F.E. Baker, S.H. Brown, J.R. Brauer and T.R. Gerhardt, Comparison of magnetic fields computed by finite element and classical series methods, Internat. J. Numer. Meths. Engrg. 19 (1983) 271-280. [19] A. Bossavit and J.C. Verne, A fixed FEM-BIEM method to solve 3-D eddy current problems, Proc. COMPUMAG, Chicago, IL, 1981; IEEE Trans. Magnetics 18(2) (1982) 431. [20] T.W. Preston and A.B.J. Reece, Solution of three-dimensional eddy current problems: the T-L! method, Proc. COMPUMAG Conf., Chicago, IL, 1981; IEEE Trans. Magnetics 18 (2) (1982) 486. [21] C.S. Biddlecombe, E.A. Heighway, J. Simkin and C.W. Trowbridge, Methods for eddy current computation in three dimensions, Proc. COMPUMAG Conf., Chicago, IL, 1981; IEEE Trans. Magnetics 18(2) (1982) 492. 1221 C.R.I. Emson and J. Simkin, An optimal method for 3-D eddy currents, Proc. COMPUMAG Conf., Genoa, 1983; IEEE Trans. Magnetics 19(6) (1983). [23] J.J. Balchin and J.A.M. Davidson, Numerical method for calculating magnetic flux and eddy current distribution in three dimensions, IEE Proc. A 127(l) (1980) 6.53. [24] D.W. Kelly, S.R. Gago, O.C. Zienkiewicz and I. Babushka, A posteriori error analysis and adaptive processes in the finite element method: Parts 1 and 2, Intemat. J. Numer. Meths. Engrg. 19 (1983) 1593. [25] J. Simkin, ed., COMPUMAG Conference Proceedings (Rutherford Appleton Laboratory, Oxford, 1976). [26] J.C. Sabonnadiere, ed., COMPUMAG Conference Proceedings (Laboratoire d’Electrotechnique de Grenoble, ERA 524 CNRS, Grenoble, 1979). [27] COMPUMAG Conference Proceedings, Chicago, IL, 1981, IEEE Trans. Magnetics 18(2) (1982). [28] COMPUMAG Conference Proceedings, Genoa, 1983, IEEE Trans. Magnetics 19(6) (1983). [29] COMPUMAG Conference Proceedings, Fort Collins, CO, 1985. [30] C.W. Trowbridge, Applications of integral equation methods for the numerical solution of magnetostatic and eddy current problems, Paper presented to Internat. Conf. Numerical Methods in Electrical and Magnetic Field Problems, Santa Margherita, Italy, 1976; see also Rutherford Laboratory Rept. RL-76-071, 1976, and Finite Elements in Electrical and Magnetic Field Problems (Wiley, Chichester, 1979) 191-213.