1. Introduction

The computation 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 anistropic 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 electromagnetic 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(&not;) if and only if || f ||2 = ??&not; | 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.

&not; 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 &not;. A triangulation of &not; consists of a set of closed triangles {Di*}i=1,...N such that

(i) Di &laquo; Dj = ? (? is the null set) for all 1£ i, j £ N, i£ j, and

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

a common side of both triangles,

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

(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 (&not;-K) / area of &not;) < 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 &iexcl;. This is illustrated both for 1 dimensional and 2 dimensional cases below.

Fig 5.2 Pyramid functions

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 &raquo; ¶ RN ,

where the Neumann boundary condition holds on ¶ RN and

the Dirichlet boundary condition 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.

where the set of functions ji(x, y) are piecewise linear polynomials, ai(t) are coefficients at a given time, and y(x, y) is chosen to satisfy the boundary conditions. The functions a(x, y, t) does not satisfy the Neumann boundary condition. The functions ji are called the interpolation, trial, basis, or shape functions in different applications.

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 &not;. Writing the function a(x, y, t) out from the definition (5.3), we have

(5.10) ò ò &not; n0 [( a1 ¶ j1/¶ x + . . . + aN ¶ jN/¶ x ) ?ji/?x

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

= ò ò &not; ( 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 = ò ò &not; n0 ( ¶ jj/¶ x ¶ ji/¶ x + ¶ jj/¶ y ¶ ji/¶ y ) dx dy

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

Ei = ò ò &not; ( s¶ a/¶ t ) ji dx dy

Ui = ò ò &not; ( -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 &ordf; (an+1 - an) / Dt,

equation (5.11) is rewritten as,

(5.20) S an+1 = D (an+1 - an) / Dt + cn+1 V

where

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

parts

D relates the conductivity 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) dan+1 = cn+1V - S an.

To solve this equation, we note that we do not have the value of cn+1 explicitly. We put dan+1 = da0, n+1 + cn+1da*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) da*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

The magnetic vector potential at the next time step is given from

(5.25) an+1 = an + da0,n+1 + cn+1 da* 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(&not;)}1/2 = O(Dr)+O(Dt)

n=1

where || ||H1(&not;) is the norm in the Hilbert space over &not; 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 Schenectady, 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 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 underlying 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 dependence 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, Schenectady, 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 Electromagnetic 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).