Trowbridge, C.W. “Three-Dimensional Analysis”
The Electrical Engineering Handbook
Ed. Richard C. Dorf
Boca Raton: CRC Press LLC, 2000
44
Three-Dimensional
Analysis
44.1Introduction
44.2The Field Equations
Low Frequency Fields?Statics Limit
44.3Numerical Methods
Finite Elements?Edge Elements?Integral Methods
44.4Modern Design Environment
44.1 Introduction
The three-dimensional analysis of electromagnetic fields requires the use of numerical techniques exploiting
the best available computer systems. The well-found laboratory will have at its disposal a range of machines
that allow interactive data processing with access to software that provides geometric modeling tools and has
sufficient central processing unit (CPU) power to solve the large (>10,000) set of algebraic equations involved,
see Section 44.4.
44.2 The Field Equations
The classical equations governing the physical behavior of electromagnetic fields over the frequency range dc
to light are Maxwell’s equations. These equations relate the magnetic flux density (B), the electric field intensity
(E), the magnetic field intensity (H), and the electric field displacement (D) with the electric charge density
(r) and electric current density (J). The field vectors are not independent since they are further related by the
material constitutive properties: B = mH, D = eE, and J = sE where m, e, and s are the material permeability,
permittivity, and conductivity, respectively. In practice these quantities may often be field dependent, and
furthermore, some materials will exhibit both anisotropic and hysteretic effects. It should be strongly stated
that accurate knowledge of the material properties is one of the most important factors in obtaining reliable
simulations.
Because the flux density vector satisfies a zero divergence condition (div B = 0), it can be expressed in terms
of a magnetic vector potential A, i.e., B = curl A, and it follows from Faraday’s law that E = –(]A/]t + ?V),
where V is the electric scalar potential. Neither A nor V is completely defined since the gradient of an arbitrary
scalar function can be added to A and the time derivative of the same function can be subtracted from V
without affecting the physical quantities E and B. These changes to A and V are the gauge transformations,
and uniqueness is usually ensured by specifying the divergence of A and sufficient boundary conditions. If
? { A = –(msV + me]V/]t) (Lorentz gauge) is selected, then the field equations in terms of A and V are:
C. W. Trowbridge
Vector Fields, Inc.
? 2000 by CRC Press LLC
(44.1)
where s has been assumed piecewise constant. This choice of gauge decouples the vector potential from the
scalar potential. For the important class of two-dimensional problems there will be only one component of A
parallel to the excitation current density. For fields involving time, at least two types can be distinguished: the
time harmonic (ac) case in which the fields are periodic at a given frequency v, i.e., A = A
o
exp(jvt), and the
transient case in which the time dependence is arbitrary.
Low Frequency Fields
In the important class of problems belonging to the low frequency limit, i.e., eddy current effects at power
frequencies, the second derivative terms with respect to time (wave terms) in Eq. (44.1) vanish. This approxi-
mation is valid if the dimensions of the material regions are small compared with the wavelength of the prescribed
fields. In such circumstances the displacement current term in Maxwell’s equations is small compared to the
free current density and there will be no radiation [Stratton, 1941]. In this case, while a full vector field solution
is necessary in the conducting regions, in free space regions, where s = 0 and curl H = J
s
Eqs. (44.1) can be
replaced by ?
2
c = 0, where c is a scalar potential defined by H = –?c. The scalar and vector field regions are
coupled together by the standard interface conditions of continuity of normal flux (B) and tangential field (H).
Statics Limit
In the statics limit (dc) the time-dependent terms in Eq. (44.1) vanish, and the field can be described entirely
by the Poisson equation in terms of a single component scalar potential, which will be economic from the
numerical point of view. In this case the defining equation is
? { m?f 5 ? { mH
s
(44.2)
where f is known as the reduced magnetic scalar potential with H = H
s
– ?f, and H
s
the source field given
by the Biot Savart law. Some care is needed in solving Eq. (44.2) numerically, in practice, as H
s
will often be
calculated to a higher accuracy than f. For instance, in regions with high permeability (e.g., ferromagnetic
materials), the total field intensity H tends to a small quantity which can lead to significant errors due to
cancellation between grad f and H
s
, depending upon how the computations are carried out. One approach
that avoids this difficulty is to use the total scalar potential c in regions that have zero current density [Simkin
and Trowbridge, 1979], i.e., where H = –?c and H
c
is the coercive field for the material c satisfies
? { m?c 5 ? { mH
c
(44.3)
Again, the two regions are coupled together by the standard interface condition that results, in this case, in a
potential “jump” obtained by integrating the tangential continuity condition, i.e.,
(44.4)
where G is the contour delineating the two regions that must not intersect a current-carrying region; otherwise
the definition of c will be violated.
?′ ?′ + + = ? ?×
é
?
ê
ê
ù
?
ú
ú
+=?×
11
2
2
2
2
m
s
?
?
?
?
m
m
?
?
ms
?
?
A
A
A
A
t
t
V
t
V
t
V
e
e
fy=+
ò
H
s
d G
G
0
? 2000 by CRC Press LLC
44.3 Numerical Methods
Numerical solutions for the field equations are now routine for a large number of problems encountered in
magnet design; these include, for example, two-dimensional models taking into account nonlinear, anisotropic,
and even hysteretic effects. Their use for complete three-dimensional models is not so widespread because of
the escalating computer resources needed as the problem size increases. Nevertheless, 3-D solutions for non-
linear statics devices are regularly obtained in industry, and time-dependent solutions are rapidly becoming
more cost effective as computer hardware architectures develop.
Finite Elements
This increasing use of computer-based solutions has come about largely because of the generality of the finite
element method (FEM). In this method, the problem space is subdivided (discretized) into finite regions
(elements) over which the solution is assumed to follow a simple local approximating trial function (shape
functions). In the simplest situation, a particular element could be a plane hexahedra defined by its eight vertices
or nodes and a solution of Eq. (44.3) may be approximated by
(44.5)
Because a hexahedra has eight nodes it is natural to select a bilinear trial function with eight parameters; see
Fig. 44.1 for other examples. The functions N
i
are called the local shape functions and the parameters U
i
are
the solution values at the nodes. The finite elements can be integrated into a numerical model for the whole
problem space either by (a) the variational method in which the total energy of the system is expressed in terms
of the finite element trial functions and then minimized to determine the best solution or (b) the weighted
residual method in which the formal error (residual), arising by substituting the trial functions into the defining
equation, is weighted by a suitably chosen function and then integrated over the problem domain. The best fit
for the trial function parameters can then be obtained by equating the integral to zero. Both methods lead to
a set of algebraic equations and are equivalent if the weighting functions are chosen to be the trial functions
(Galerkin’s method [Zienkiewicz, 1990]). At the element level, the residual R
i
is given by
(44.6)
where Q (RHS) denotes the sources of Eqs. (44.2) or (44.3). The integrals can be readily evaluated and assembled
for the whole problem by superposition, taking account of the boundary conditions and removing the redun-
dancy at shared nodes. At interelement boundaries in a region of particular potential [reduced Eq. (44.2) or
total Eq. (44.3)] the solution is forced to be continuous, but the normal flux (i.e., m]U/]n) will only be
continuous in a weak sense, that is to say the discontinuity will depend upon the mesh refinement.
The FEM provides a systematic technique for replacing the continuum field equations with a set of discrete
algebraic equations that can be solved by standard procedures. In Fig. 44.2 a typical field map is shown for a
permanent magnet machine modeled by a computer simulator that can take into account nonlinearity and
permanently magnetized materials. Although hysteresis effects can be included, the computational resources
required can be prohibitive because of the vector nature of magnetization. The magnetic material must be
characterized by a large number of measurements to take account of the minor loops, and from these the
convolution integrals necessary to obtain the constitutive relationships can be evaluated [Mayergoyz, 1990].
These characteristics must then be followed through time; this can be implemented by solving at a discrete set
of time points, given the initial conditions in the material.
yaaaaaaaa?=++++ + + + =
?
u x y z xy yz zx xyz NU
ii12 345 6 7 8
RNNdUNQd
iijji
=??
é
?
ê
ê
ù
?
ú
ú
+
òò
mW W
elem elem
? 2000 by CRC Press LLC
Although the FEM is widely used by industry for electromagnetic problems covering the entire frequency
range, there are many situations where special methods are more effective. This is particularly the case for high-
frequency problems, e.g., millimeter and microwave integrated circuit structures where integral equation
techniques and such procedures as transmission line modeling (TLM), spectral domain approach, method of
lines, and wire grid methods are often preferred [Itoh, 1989] (see Chapter 43).
Edge Elements
Using potentials and nodal finite elements (see Fig. 44.1) rather than field components directly has the advantage
that difficulties arising from field discontinuities at material interfaces can be avoided. However, if the element
basis functions [see Eq. (44.5)] are expressed in terms of the field (H, say) constrained along an element edge,
then tangential field continuity is enforced [Bossavit, 1988]. The field equations [Eq. (44.1)] in terms of the
field intensity for the low frequency limit reduce to
(44.7)
FIGURE 44.1 Three-dimensional second-order Isoparametric finite elements. (a) Left, master tetrahedron, 10 nodes in
local space (j, h, z); right, actual tetrahedron, 10 nodes in global space (x, y, z). (b) Left, master hexahedron, 20 nodes in
local space (j, h, z); right, actual hexahedron, 10 nodes in global space (x, y, z).
?′?′ + =H
H
s
?m
?
()
t
0
? 2000 by CRC Press LLC
and a suitable edge variable basis function form for solving this equation by finite elements using tetrahedral
elements is
h(r) = a + b ′ r (44.8)
where r is the position vector and a and b, respectively, are vectors dependent on the geometry of the element.
The basis function expansion is given by
H =
(
h
e
(r)H
e
(44.9)
where h
e
is the vector basis function for edge e, and H
e
is the value of the field along an element edge (see
Fig. 44.3). The functions, Eqs. (44.8) and (44.9), have the property of being divergence free, and most important
they ensure that the tangential component of H is continuous while allowing for the possibility of a discontinuity
in the normal component. In nonconducting regions where the field can be economically modeled by a scalar
potential, standard nodal elements can be used. At the interface the edge elements couple exactly with the nodal
elements.
Integral Methods
An alternative procedure is to solve the field equations in their integral form, see also Chapter 43. For example,
in magnetostatics, the magnetization vector M given by M = (m – 1)H can be used instead of H to derive an
integral equation over all ferromagnetic domains of the problem, i.e.,
FIGURE 44.2 Permanent magnet motor.
? 2000 by CRC Press LLC
(44.10)
where R is the distance between the source and field point, respectively. For problems with linear materials
Eq. (44.10) reduces to integrations over the bounding surfaces of materials in terms of the magnetic scalar
potential, i.e.,
(44.11)
Equation (44.11) can be solved numerically by the boundary element method (BEM) in which the active
surfaces are discretized into elements. The advantages of integral formulations compared to the standard
differential approach using finite elements are (a) only active regions need to be discretized, (b) the far field
boundary condition is automatically taken into account, and (c) the fields recovered from the solution are
usually very smooth. Unfortunately, the computational costs rapidly escalate as the problem size increases
because of the complexity of the system coefficients and because the resulting matrix is fully populated, whereas
in the differential approach the coefficients are simple and the matrix is sparse, allowing the exploitation of
fast equation solution methods.
44.4 Modern Design Environment
The most common system used in software packages is one in which the pre-processor includes data input,
model building, and mesh (element) generation. Although fully automated meshing is now a practical possibility
it needs to be combined with error estimation in order to allow the generation of optimal meshes. This approach
is now common for 2-D systems and is available in many 3-D systems. Figure 44.4 illustrates a field simulation
environment in which the solution processor includes an adaptive mesh generator controlled by a posteriori
error estimation. This avoids the costly and essentially heuristic task of mesh generation which, in the past,
had to be performed by the designer. Furthermore a modern system should have solid modeling capabilities
driven by parametric data allowing the user to specify the appropriate engineering quantities, e.g., in the case
of a solid cylinder the radius and length are all that is needed to specify the geometry at some predefined
location. The system should also be supported by a database which is compliant with evolving standards such
as STEP (Standard for the Exchange of Product data-ISO 10303 [Owen, 1993]) thus allowing data communi-
cation between other systems.
The environment illustrated in Fig. 44.4 also shows tools for automatic optimization that are now becoming
feasible in industrial design applications. Both deterministic and stochastic methods for minimizing constrained
objective functions of the design space have been developed for electromagnetic applications (For a review see
Russenschuck, 1996). It must be emphasized, however, that the use of optimizing methods is only part of the
total problem of design. For example, the process of automatic synthesis based on design rules and engineering
FIGURE 44.3Edge variables for a tetrahedron element, h
1
= (F
2
– F
1
)/l.
Mr Hr M()( ) (=- ￠-?×?
?
è
?
?
?
÷
é
?
ê
ê
ù
?
ú
ú
ò
m
p
1
1
s
R
d)
4
1
W
W
4
11
pf
?f
?
f
?
?
=- -
?
è
?
?
?
÷ò
Rn
R
n
d
/
G
G
? 2000 by CRC Press LLC
knowledge may provide engineers with a complementary methodology to assist in the creativity that is the
essence of design (see Lowther, 1996). An example of a finite element solution displayed on a modern work-
station for a three-dimensional eddy current problem with a moving conductor retarded by a c-core electro-
magnet is shown in Fig. 44.5. The post-processor indicates the magnitude and direction of the induced eddy
FIGURE 44.4 Electromagnetic design environment.
FIGURE 44.5 Moving copper strip, retarded by an electromagnet.
? 2000 by CRC Press LLC
currents in the moving copper strip by the solid cones and the magnetic field by the gray scaled contours. This
modest problem was modeled with 5,664 degrees of freedom (No. of equations) and needed 198 cp seconds
running on a workstation rated at SPECfp92 96.5. However most industrial problems will require many more
degrees of freedom, and typically a non-linear magnetostatic problem with 200,000 equations needed
75 mintues and 25 Mbytes of RAM on the same machine. The resources needed for transient non-linear
problems will be far greater.
Defining Terms
Biot Savart law:
where R is the distance from the source point to the field point.
Interface conditions: (B
2
– B
1
) { n = 0
(D
2
– D
1
) { n = v
(H
2
– H
1
) 3 n = K
(E
2
– E
1
) 3 n = 0
where K and v are the surface current and charge densities, respectively.
Maxwell’s equations:
Related Topics
35.1 Maxwell Equations ? 45.1 Introduction ? 45.3 Analytical Issues in Developing a Computer Model
References
A. Bossavit, “Rationale for ‘edge elements’ in 3-D field computation,” IEEE Trans. on Magnetics, vol. 24, no. 1,
p. 74, 1988.
I. Tasuo (Ed.), Numerical Techniques for Microwave and Millimeter-Wave Passive Structures, New York: John
Wiley, 1989.
D. A. Lowther, “Knowledge-based and numerical optimization techniques for the design of electromagnetic
devices,” IJ Num. Mod., vol. 9, no. 1,2, pp. 35–44, 1996.
I. Mayergoyz, Mathematical Models of Hysteresis, New York: Springer-Verlag, 1990.
J. Owen, STEP An Introduction, Information Geometers Ltd, 47 Sockers Avenue, Winchester, UK, 1993.
S. Russenschuck, “Synthesis, inverse problems and optimization in computational electromagnetics,” IJ Num.
Mod., vol. 9, no. 1,2, pp. 45–58, 1996.
P. P. Silvester and R. L. Ferrari, Finite Elements for Electrical Engineers, 2nd ed., Cambridge: Cambridge University
Press, 1990.
HJ
ss
R
d=′?
ò
1
4
1
p
W
W
?× =
?× =
?′ =-
?′ = +
D (Gauss's law)
B
E (Faraday's law)
H J (Ampere's law + displacement current)
r
?
?
?
?
0
B
D
t
t
? 2000 by CRC Press LLC
J. Simkin and C. W. Trowbridge, “On the use of the total scalar potential in the numerical solution of field
problems in electromagnetics,” IJNME, vol. 14, p. 423, 1979.
J. A. Stratton, Electromagnetic Theory, New York: McGraw Hill, 1941.
O. C. Zienkiewicz, The Finite Element Method, 3rd ed., New York: McGraw Hill, 1990.
Further Information
Conferences on Computation of Electromagnetic Fields Compumag Proceedings:
Oxford, UK, 1976 (Ed. J. Simkin) Rutherford Appleton Laboratory, Chilton, Oxon, UK.
Grenoble, France, 1979 (Ed. J. C. Sabonnadiere) ERA 524 CNRS, Grenoble, France.
Chicago, USA, 1981 (Ed. L. Turner) IEEE Trans. Mag. 18 (2) 1982.
Genoa, Italy, 1983 (Ed. G. Molinari) IEEE Trans. Mag. 19 (6) 1983.
Fort Collins, USA, 1985 (Ed. W. Lord) IEEE Trans. Mag. 21 (6) 1985.
Graz, Austria, 1987 (Ed. K. Richter) IEEE Trans. Mag. 24 (1) 1988.
Tokyo, Japan, 1989 (Ed. K. Miya) IEEE Trans. Mag. 26 (2) 1990.
Sorrento, Italy, 1991 (Ed. R. Martone) IEEE Trans. Mag. 28 (2) 1992.
Miami, USA, 1993 (Ed. D. A. Lowther) IEEE Trans. Mag. 30 (5) 1994.
Berlin, Germany, 1995 (Ed. O Biro) IIEEE Trans. Mag. To be published May 1996.
Conference on Electromagnetic Field Computation, CEFC Proceedings
Washington, USA, 1988 (Ed. I. Mayergoyz), IEEE Trans. Mag. 25 (4), 1989.
Toronto, Canada, 1990 (Ed. J. Lavers), IEEE Trans. Mag. 27 (5), 1991.
Claremont, USA, 1992 (Ed. S. R. Hoole), IEEE Trans. Mag.
Aix-les-Bains, France, 1994 (Ed. J. C. Sabonnadiere), IEEE Trans. Mag. 31 (3), 1995.
Special Issue on Computational Magnetics (Eds. J. Sykulski and P. Silvester)
Int. Jour. Of Num. Mod. J. Wiley, vol. 9, no. 1 & 2, 1996.
? 2000 by CRC Press LLC
