__The compution of the current density distribution for a cylindrical
conductor can be carried out with the aid of the finite element method.
When the length of the conductor is large as compared with its width, the
axial dimension of the problem may be neglected and solved in two-dimensional
space. The finite element method is presented from the approach of weighted
residuals. In various applications, the total current passing through the
conductor may be given as a non-periodic function in time. Unless there
are factors to the contrary, the time derivatives may be simply approximated
using finite differences rather than finite elements.__

__The computation of the current density and eddy current losses distribution
in an infinite copper cylinder with a z-directed time varying current is
often an essential part in the design of electrical equipment.__

__The eddy current densities are often among the most complicated to
analyze and difficult to predict. There are many cases in which the conducting
material will have nonisotropic nonlinear anisoltropic characteristics
as well as irregular geometric boundaries. The eddy current losses caused
by nonuniform current density distribution are often significant. A great
deal of effort has gone into the finite element method because of its success
in many area. Some of the most popular books are [5.1], [5.2], [5.3], [5.4],
[5.5], [5.6], and [5.7].__

__The problem of an infinite circular cylinder was recently studied
by Sikora [5.8], however only in the case of sinusoidal current flowing
in the conductor. Moreover, his partial differential equations are written
in terms of the current density and magnetic field intensity vector. Progress
has been made in the solution of eddy currents by expressing the electomagnetic
field in terms of the magnetic vector potential whose curl is the flux
density.__

__The two-dimensional problem using the finite element method approach
will be solved. The currents are confined to the z-direction, but vary
in the radial and circumferential direction.__

__2. The Finite Element Method for a Two-Dimensional Problem__

__In the finite element method, we solve a problem in two-dimensional
real space, over a not necessarily bounded domain __R,
such that (x, y) Œ
R. We approximate this region by ¬, a
closed bounded convex set, whose boundary ?¬
is a piecewise differentiable and rectifiable curve. Over the set ¬,
we consider a set ?n(x,
y) of basis functions which are piecewise continuous polynomials of degree
n over a set of polygons {¡a},
where ¡a is
contained in ¬,
and ja??n(x,
y) vanishes outside ¡a
(i. e. the support of ja
is the polygon ¡a).
We assume that ?n
is a complete set of functions over the class of L2(¬).
The class of L2
functions are such that

f Œ L2(¬) if and only if || f ||2 = ??¬ | f(x, y) |2 dx dy < •.

A set is complete over a class of functions if it can approximate any function in the class to within arbitrary accuracy using the appropriate norm. Since we are only dealing with first degree polynomials, ? will denote ?1.

¬ is restricted in practice to be a set for which the ratio of the maximum width to the minimum width is less than 5 for double precision calculations. A convex set has width d in a direction q if it can be covered by a strip of width d and by no strip of smaller width. A strip of width d in a direction q is the region between two parallel lines a distance d apart making an angle q with respect to the coordinate system. The reason for this restriction is so that long thin regions with sharp corners will not be modeled.

Let Di be any open triangle and Di* its closure contained in ¬. A triangulation of ¬ consists of a set of closed triangles {Di*}i=1,...N such that

(i) Di « Dj = ? (? is the null set) for all 1? i, j ? N, i?j, and

if Di* « Dj* ? ? (i?j), then Di* « Dj* is either a point vertex or

a common side of both triangles,

(ii) K = » i=1,...N Di* is a closed convex set whose boundary ?K is a polygon whose vertices belong to ?¬,

(iii) No interior angle of any Di is less than 10° (this is a practical

requirement because extending this criteria results in an

ill-conditioned system of linear equations),

(iv) a vertex of a triangle Di* is the vertex of an adjacent triangle Dj* and cannot lay on the side of Dj*, and

(v) (Area of (¬-K) / area of ¬) < c, where c < 0.05 in practice.

The reader is referred to Ciarlet [5.9] for details. The vertices of the triangles are the nodes, and the triangles will be the finite elements. On the elements, polynomial interpolations will be used to represent the values of the function between the nodes because polynomials are easy to calculate and manipulate. In particular, we will use linear interpolations on the elements.

In one dimension, the first order element is chosen as a continuous piecewise linear function between adjacent nodes as illustrated below. In a two-dimensional space domain, the most commonly used finite elements are the first order triangular finite elements. The term "first order" refers to the order of the polynomial approximation, while "triangular" refers to the shape of the elements. The nodal points are the vertices of the triangular element.

Fig. 5.1 Examples of elements in a grid.

Given a nodal point, pi, we define a weighting function to be a function whose value is 1 at the nodal point, pi, and 0 at all the other nodal points. On each element, the function will be linear and will match the values at the nodes of the element. This defines the function jj over K. The weighting function associated with a node looks like a pyramid centered at that node and whose base is the polygon ¡. This is illustrated both for 1 dimensional and 2 dimensional cases below.

ji(x) = Ì (x-xi+1)/(xi - xi+1) xi ? x ? xi+1

Ó 0 elsewhere

The shape of these function gives rise to
the name *hat* or *chapeau* or *pyramid* functions [5.10]

We have a weighting function jj centered at each of the nodes, pj. If each node is assigned a value aj, then SMj=1 aj jj gives a piecewise linear surface made up of many triangular faces.

Consider the scalar partial differential equation in A(x, y, t)

(5.1) n0—2A = - Js(x, y, t) + s?A/?t (x, y) ŒR, t Œ (0, ?)

subject to the current constraint

(5.2) C(t) ?F s dx dy - ?F s ?A/?t dx dy = I(t)

with boundary conditions

A(x, y, t)
= Y(x, y) (x,
y) ???RD,
t Œ
[0, ?) ?A(x, y, t)/?n = 0 (x, y) ???RN,**
**t Œ
[0, ?)

A(x, y, 0) = A0(x, y) (x, y) ŒR,

where

A(x, y, t) is the z-component of the magnetic vector potential

which is assumed to be continuous and twice differentiable,

Js(x, y, t) = s C(t), the source current density,

C(t) is an unknown for which to solve and must satisfy (5.2),

(see the section on finite differences),

I(t) is the given total current in the conductor,

a piecewise linear continuous function,

Y(x, y) is a prescribed function on ?RD,

representing Dirichlet boundary conditions,

A0(x, y) is the initial condition for A,

?R is the boundary of the domain, ?R = ?RD » ?RN ,

where the Neumann boundary condition holds on ?RN and

the Dirichlet boundary condtion holds on ?RD,

F is the cross section of the conductor perpendicular

to the z-axis,

s(x, y, T(x, y, t)) is the conductivity,

s varies with the temperature T

n0 is the reluctivity of free space, and

n denotes the outward normal to the boundary curve ?RN.

We assume the solution A(x, y, t) exists and approximate it by

(5.3) a(x, y, t) = y(x, y) + Â ai(t) ji(x, y)

i=1

We define the residual

(5.4) r (x, y, t) = n0—2a(x, y, t) + Js(x, y, t) - s ?a(x, y, t)/?t.

For fixed t, we will consider how to minimize r(x, y, t) by choosing the coefficients ai(t). Since the class ? is complete in L2, it can be shown that there exists a sequence of rk(x, y, t) with suitable jk,j(x, y) such that

(5.5) ??R rk(x, y, t) wi(x, y) dx dy -> 0 i=1, ..., •

where k is an index on the ever increasingly fine meshes and

the weighting or test functions wi are chosen from ?. However, for practical purposes, we will choose a finite set of the wi, say {wi}, 1?i?N, and require

(5.6) ??R r(x, y, t) wi(x, y) dx dy = 0 i = 1, ..., N.

This procedure is known as *the method of
weighted residuals. *In essence, what we are doing is defining that
each projection onto our weighting function is zero. It is equivalent to
requiring that r(x, y, t) vanishes in a finite dimensional sub-space of
L2(R).
Since ? is
complete and has a countable dense subset, it is theoretically possible
to achieve any level of accuracy in the L2
norm by choosing a sufficient number of representative members of the set
{wi}, see Rudin [5.11]

The way in which the weighting functions are chosen gives rise to various different methods, such as the moment, Galerkin, subdomain, or collocation method. For example, when the collocation method is chosen, the weighting function is the Dirac delta function, and the scheme is similar to the finite difference scheme.

In the *Galerkin method*, the weighting
functions are chosen to be

(5.7) wi = ji , i = 1, 2, ... N,

the same as the basis functions. Rewriting r(x, y, t) in terms of a(x, y, t), equation (5.6) becomes

(5.8) ??R n0 —2a(x, y, t) ji(x, y) dx dy

= ??R [- J(x, y, t) + s a(x, y, t)/?t ] ji(x, y) dx dy

i=1, 2, . . . , N.

Integrating the left hand side by parts to obtain symmetry between the functions a and ji, we obtain

(5.9)??R n0 (?a/?x ?ji/?x + ?a/?y ?ji/?y ) dx dy

+ ??RD (?a/?x ji + ?a/?y ji ) dl

+ ??RN (?a/?x ji + ?a/?y ji ) dl

= ??R ( s?a/?t - Js(x, y, t) ) ji dx dy i=1, ..., N.

The second integral on the left hand side vanishes because over the path ?RD, the residual r(x, y, t) is zero since a(x, y, t) satisfies the boundary conditions exactly. The integral over the path ?RN is set equal to zero explicitly in the weighted residual formulation or in a "weak" sense, by requiring that

???RN(?a/?x + ?a/?y)ji dl = 0.

This results in an additional constraint which must be imposed, see Zienkiewicz [5.12] However, it has been shown, see Ciarlet [5.13], that this last integral does not have to be imposed as a separate equation to be enforces but rather arises naturally in the finite element formulation. We now restrict the region which we will model the magnetic vector potential from R to its approximation ¬. Writing the function a(x, y, t) out from the definition (5.3), we have

(5.10) ??¬ n0 [( a1 ?j1/?x + . . . + aN ?jN/?x ) ?ji/?x

+ ( a1 ?j1/?y + . . . + aN ?jN/?y ) ?ji/?y ] dx dy

= ??¬ ( s?a/?t - Js(x, y, t) ) ji dx dy for i = 1, ... N.

The above equation can be written in matrix form as

(5.11) **S** __a__ = __E__ + __U__

where

Sij = ??¬ n0 ( ?jj/?x ?ji/?x + ?jj/?y ?ji/?y ) dx dy

__a__ = (a1,
a2,
. . . , aj, . . . , aN)t

Ei = ??¬ ( s?a/?t ) ji dx dy

Ui = ??¬ ( -Js(x, y, t) ) ji dx dy

1 ?i, j ? N

For each element, we compute the local contribution to the problem. The reason for this is that if two nodes do not lie on a common finite element, then the product of their weighting functions is zero. This fact is easily seen if we recall that the weighting function is zero outside the polygon formed by its nearest neighbors. This contribution is described as entries in a local matrix Se which will then be embedded in a global matrix Sg.

(5.12) Sije = ? n0 ( ?jj/?x ?ji/?x + ?jj/?y ?ji/?y ) dx dy

e

where e denotes an arbitrary element.

It can be shown (see Silvester [5.14]) that the local matrix Se will be of

the form:

| bi2 + ci2 bibj + cicj bibm + cicm |

(5.13) Se = n0 /4 D | bibj + cicj bj2 + cj2 bjbm + cjcm |

| bibm + cicm bjbm + cjcm bm2 + cm2 |

where

(5.12) (xi, yi) is the coordinate of the ith node of element e,

bi = ym - yj

ci = xm - xj in cyclic modulo 3, with the permutation (i, j, m),

De is the area of the triangle e,

i, j, and m are indices on the vertices of the triangle e arranged in

clockwise order.

The value of a(x, y, t) within the element e is given as

(5.14) ae(x, y, t) = ai(t) ji(x, y) + aj(t) jj(x, y) + am(t) jm(x, y)

namely this is a linear combination of the values on the surrounding vertices. Similarly we have

(5.15) ?ae(x, y, t)/?t =

?ai(t)/?t ji(x, y) + ?aj(t)/?t jj(x, y) + ?am(t)/?t jm(x, y)

where i, j, and m are indices on the nodes (vertices) of the element (triangle) e.

We make the following approximation

(5.16) Ue = ??e Js(x, y, t) dx dy = ??es(x, y) C(t) dx dy ? se C(t) ?e

where C(t) is as defined above (5.2). se is the electrical conductivity of the material filling the element e.

(5.17) Ee = ??e s ?a/?t dx dy ???e se ?a/?t dx dy =

1/3 se (?ai/?t + ?aj/?t + ?am/?t) ?e

For every element e, the associated local matrix Se relates the nodes to the element. The global matrix Sg relates all the nodes and elements to model the geometry of our problem. We form a table that relates the element numbers with the global numbering of the nodes. Using this information, it is possible to embed the matrix Se into the matrix Sg. Each entry in the local matrix is added to an entry in the global matrix so that it multiplies the appropriate entry in the unknown vector (as it did in the local matrix). The reader is referred to Akin [5.10] or Zienkiewicz [5.12] for worked out examples of this procedure.

There are two types of boundary conditions
with which we will be concerned. The Dirichlet boundary condition A = AD
is enforced by assigning specific values for certain components of the
vector __a__. One row for each such ai can be removed
from the matrix equation. The value of the function

a(x, y, t) assigned by this boundary condition is included in the y term. It is of the form

(5.18) y(x, y) = S a'i(x, y) ji(x, y),

where a'i(x, y) is independent of time for our problem. The Neumann boundary condition ?A/?n = 0, is imposed implicitly in the weighted residual formulation.

It is possible to use finite elements over the time variable. Kamar [5.15] shows that for certain choices of the weighting function in time, a finite element formulation in time reduces to a finite difference formulation. However, we will use finite differences in time as there is no specific advantages for using finite elements. Using the backwards difference approximation

(5.19) ?__a__/?t ª
(__a__n+1 - __a__n) / Dt,

equation (5.11) is rewritten as,

(5.20)** **S __a__n+1
= D (__a__n+1 - __a__n)
/ Dt + cn+1
__V__

where

S accounts for the n0—2 term after it has been integrated by

parts

D relates the conductivities to correspond with the proper terms for the magnetic vector potential, E = D (an+1 - an) / Dt

__V__ relates the position of the source
current term with the proper element, __V__ = (1/cn+1)
__U__

an+1 is a column vector which holds the nodal values of the

magnetic vector potential at time tn+1,

cn+1 is a scalar constant for which to solve, C(tn+1) = cn+1.

This equation along with the current constraint (5.2) is to be solved. At every time step, Bedrosian [5.16] chooses to use the substitution dan+1 = an+1 - an to get

(5.21) (S - D/Dt)
d__a__n+1
= cn+1__V____ - S a__n.

To solve this equation, we note that we do
not have the value of cn+1 explicitly. We put d__a__n+1
= d__a__0,
n+1
+ cn+1d__a*__1,
n+1 ,
and rewrite the above equation as

(5.22) (S - D/Dt) da0, n+1 = - S an and

(5.23) (S - D/Dt)
d__a*__1,
n+1
= c*n+1 __V__

where c*n+1=1. The advantage of this form is that it is possible to separate the change in the magnetic vector potential from two different influences. The effect described by equation (5.22) is the "free running solution" that describes the change in the magnetic vector potential as if it were set to some initial value (namely that what it was at the time = tn) and looked at a time Dt later (time = tn+1) without the effect of the source term. The effect described by equation (5.23) is the response of the magnetic vector potential due to induced currents by an electric field applied over the conductor. We treat equation (5.20) as if it were linear over each time step. The absolute magnitude of cn+1 is not important in obtaining dan because we scale this effect in order to get the correct total dan+1. Now the cn+1 must be calculated to give the correct value of the source current over the conductor. This is given by rearranging the current constraint equation

(5.24) C(tn+1) = —————————————————————————————

?F s dx dy - ?F s da1,n+1/Dt dx dy

(5.25) __a__n+1
= __a__n + d__a__0,n+1
+ cn+1 d__a*
__1,n+1

Now we update all the quantities (namely the temperature, the conductivity and the heat capacity) and repeat the process over the next time step. This procedure is identical to that described in the finite difference formulation.

Along with the development of the finite element algorithm, several authors have tried to investigate the convergence properties of these methods. The most important work in this area has been done by Babuska [5.17]. Several other texts have been written on this subject and the interested reader is referred to Oden and Carey [5.18], Thomeé [5.19], and Oden and Reddy[5.20]

For the case of a sinusoidal current I(t), Zlamal [5.21] gives an error estimate for a similar formulation of our problem

N

(5.26) {Dt S ||an - An||2H1(¬)}1/2 = O(Dr)+O(Dt)

n=1

where || ||H1(¬) is the norm in the Hilbert space over ¬ of the function

A(x, y, t) and its first spatial derivatives and n is an index on time. Unfortunately, the details would require much more elaboration and are presented in the reference.

The described algorithm has been written as the computer program TDTEMP by Bedrosian at the General Electric Research and Development Corporation in Schnectady, New York. This program was used to obtain all the finite element data. The presentation here is different in that his formulation is based upon the the ideas from functionals and the method of variationals while-as the presented formulation uses the method of weighted residuals. One advantage of the method of weighted residuals is that it is easier to describe to a non-specialist to the field. The approach via functionals requires identifying the interesting partial differential equation as the Euler equation to a variational problem. Another approach is to identify the energy stored in the system with the desirable functional. On the other hand, to fully understand why the method of weighted residuals works requires an understanding of the completeness of L2. Moreover, further theoretical study would require substantially more functional analysis, especially when one wishes to understand much of the work on error analysis. The underlaying mathematics occurs sooner in the development of the variational approach.

The matrix equations are solved using pre-conditioned conjugate gradient methods. The subroutine to solve the linear system was written by Silvester at the Imperial College in England. A copy is available though the electrical engineering department at McGill University. This is the standard way in which to solve large sparse systems of equation in many unknowns. This topic is not the main study here, and the interested reader is referred to Steward [5.22] or Golub and Van Loan [5.23] for more details.

Hence if we know the material properties, s and n, and the dimensions of the conductor, this model can yield useful results of the two-dimensional current density and magnetic vector potential to a fairly accurate degree in order to predict the performance of engineering devices.

__3. Conclusion__

In this presentation, two-dimensional transient fields and eddy current problems are described by parabolic differential equations in terms of the magnetic vector potential. The numerical solution is achieved by application of the finite element method, where space was discretized with first order triangles and time was discretized with finite differences. This procedure allows the consideration of field and temperature dependence of the electrical conductivity and can be used to model arbitrary regions.

The calculations made show the usefulness of the finite element method in determining the above field quantities. Graphs are provided which illustrate the dependance of the data upon the grid size and the time steps.

Fig. 5.3 Finite element calculations at 0.5 msec.

1) Analytical solution

2) FE @ Dt = 0.500 msec, Dr = R/16

3) FE @ Dt = 0.125 msec, Dr = R/16

4) FE @ Dt = 0.500 msec, Dr = R/32

5) FE @ Dt = 0.125 msec, Dr = R/32

Fig. 5.4 Finite element calculations at 0.5 msec. (close up)

1) Analytical solution

2) FE @ Dt = 0.500 msec, Dr = R/16

3) FE @ Dt = 0.125 msec, Dr = R/16

4) FE @ Dt = 0.500 msec, Dr = R/32

5) FE @ Dt = 0.125 msec, Dr = R/32

Fig. 5.5 Finite element calculations at 0.5 msec.

1) Analytical solution

2) FE @ Dt = 0.125 msec, Dr = R/16

3) FE @ Dt = 0.125 msec, Dr = R/32

Fig. 5.6 Finite element calculations at 0.5 msec. (close up)

1) Analytical solution

2) FE @ Dt = 0.125 msec, Dr = R/16

3) FE @ Dt = 0.125 msec, Dr = R/32

Fig. 5.7 Finite element calculations at 1.0 msec.

1) Analytical solution

2) FE @ Dt = 0.500 msec, Dr = R/16

3) FE @ Dt = 0.125 msec, Dr = R/16

4) FE @ Dt = 0.500 msec, Dr = R/32

5) FE @ Dt = 0.125 msec, Dr = R/32

Fig. 5.8 Finite element calculations at 1.0 msec. (close up)

1) Analytical solution

2) FE @ Dt = 0.500 msec, Dr = R/16

3) FE @ Dt = 0.125 msec, Dr = R/16

4) FE @ Dt = 0.500 msec, Dr = R/32

5) FE @ Dt = 0.125 msec, Dr = R/32

Fig. 5.9 The magnetic vector potential vs. the radial distance.

Finite element calculations at 0.5 msec.

1) FE @ Dt = 0.500 msec, Dr = R/32 (------ ------)

2) FE @ Dt = 0.250 msec, Dr = R/32 (- - - - - - - - )

3) FE @ Dt = 0.125 msec, Dr = R/32 (--- --- ---)

Fig. 5.10 The magnetic vector potential vs. the radial distance.

Finite element calculations at 1.0 msec.

1) FE @ Dt = 0.500 msec, Dr = R/32 (------ ------)

2) FE @ Dt = 0.250 msec, Dr = R/32 (- - - - - - - - )

3) FE @ Dt = 0.125 msec, Dr = R/32 (--- --- ---)

__References.__

__[5.1] Tandon, S.C. and Chari, M.V.K. "Transient Solution of the Diffusion
Equation by the Finite Element Method," Journal of Applied Physics
(53), 3, March 1981, pp. 2431-2432.__

__[5.2] Becher, E. B. and Carey, G. F. and Oden, J. T. , Finite
Elements:__

*An Introduction, Vol. I*, Prentice - Hall, (1981)

__[5.2] Konrad A., Chari, M.V.K., and Csendes, Z.J. "New Finite Elements,"
IEEE MAG.-18, No. 2, March 1982, pp. 450-453 (refer to eqs 6-11
on p. 451).__

__[5.3] Tandon, S.C., Armor, A.F., and Chari, M.V.K., "Nonlinear Transient
Finite Element Field Computation for Electrical Machines and Devices,"
IEEE Transactions on Power Apparatus and Systems. Vol. PAS-102,
No. 5, May 1983__

__[5.4] Jain, M. K., Numerical Solution of Differential Equations,
2nd Edition, A Halsted Press Book, John Wiley and Sons, (1984)__

__[5.4] Garg, V. K., Weiss, J., and Del Vecchio, R. M., "A Parametric
Study of Electromagnetic Launcher Rails Using the WEMAP System." Workshop
on Electromagnetic Field Computation. IEEE Region 1, Schnectady, New York
12345, October 20 & 21, 1986.__

__[5.5] Fletcher, C. A. J., Computational Galerkin Methods.
Springer Verlag,, New York (1984)__

__[5.5] Weiss, J. and Garg, V. K. "Finite Element Solution of Transient
Eddy Current Problems in Multiply-Excited Magnetic Systems." Westinghouse
Report 1986.__

__[5.6] Weiss, J. and Csendes, Z. J. "A One-Step Finite Element Method
for Multiconductor Skin Effect Problems," IEEE PAS Winter Meeting, NY 82WM102-2,
1982.__

__[5.7] Melkes, F. and Zlamal, M., "Numerical Solution of Nonlinear
Quasi-Stationary Magnetic Fields," International Journal for Numerical
Methods in Engineering, Vol. 19, (1983), pp.1053-1062__

__[5.8] Starzynski, J., Sikora, J., and Wincenciak, S., "Using the
Finite Element Method to Teach to Theory of Electomagnetic Fields in the
Classroom," Elektrotechnicky Obzor, Prague, Vol. 75 (1986), No. 5 -6, pp.
270-273.__

__[5.8] Bathe K.J. and Wilson E.L., Numerical Methods in Finite
Element Analysis. Prentice-Hall, Englewood Cliffs, New Jersey (1976).__

__[5.9] Ciarlet, Ph. G., The Finite Element Method for Elliptic
Problems. North-Holland Publ. Co. Amsterdam-New York (1978).__

__[5.9] Ikeda T., Maximum Principles in Finite Element Models for
Convection-Diffusion Phenomena. North-Holland Publ. Co., Amsterdam
New York (1983).__

[5.10] Akin, J. E. , *Application
and Implementation of Finite Element Methods*. Academic Press, Inc.,
(1982).

[5.11] Rudin, W., *Real and
Complex Analysis*, 2nd Edition, McGraw-Hill, NY,
(1974).

[5.11] Webb, John, "Lecture Notes for Finite Elements in Electromagnetics," McGill University, Montreal, (1986).

[5.12] Strang, G. and Fix, G.
*An Analysis of the Finite Element Method*. Prentice-Hall, Inc. Englewood
Cliffs, New Jersey (1973).

[5.12] Zienkiewicz, *The Finite
Element Method, *McGraw Hill, London, 3rd ed., (1977)

[5.13] Ciarlet, Ph. G., *Numerical
Analysis of the Finite Element Method*, Les Presses de l'university
de Montreal, C.P. 6128, Succ. "A", Que., Canada H3C-3J7, (1976).

[5.14] Silvester P.P. and Ferrari,
R.L., *Finite Elements for Electrical Engineers*, Cambridge University
Press, (1983).

[5.15] Kamar, A.M. , "Solution
fo Nonlinear Eddy Current Problems Using Residual Finite Element Methods
for Space and Time Discretization," *IEEE Transactions on Magnetics*,
Vol. MAG-19, No. 5, September 1983.

[5.16] Bedrosian, G., Private communication on TDTEMP, Nov 1984, (see Perry)

[5.17] Bieterman, M. and Babuska, I., "The Finite Element Method for Parabolic Equations: I. A Poseriori Error Estimations," Numerische Mathematik 40 pp.39-371 (1982)

[5.17] Zlamal, M., "Finite Element
Solution of Quasistationary Nonlinear Magnetic Fields," *R.A.I.R.O.-Numerical
Analysis*, Vol. 16. No. 2 (1982) pp. 161-191.

[5.18] Oden, J. T. and Carey,
G. F., *Finite Elements : Mathematical Aspects, Vol. IV*, Prentice
- Hall, (1983)

[5.19] Thomee, V., *Galerkin
Finite Element Methods for Parabolic Problems*, Springer-Verlag, (1984).

[5.20] Oden, J. T. & Reddy,
J. N., *An Introduction to the Mathematical Theory of Finite Elements*,
John Wiley and Sons, 1976.

[5.21] Strang, G. *Introduction
to Applied Mathematics*, Wellesley-Cambridge Press, Wellesley, Mass.,
(1986).

[5.22] Zienkiewicz, O. C. and
Morgan, K. *Finite Elements and Approximations*, John Wiley and Sons,
New York (1983)

[5.22] Steward, G. W., Introduction to Matrix Computations, Academic Press, NY (1973)

[5.23] Golub, G. H. and Van Loan, C.F., Matrix Computations, The John Hopkins University Press (1983).