Trowbridge, C W."Three-Dimensional Analysis The Electrical Engineering Handbook Ed. Richard C. Dorf Boca raton crc Press llc. 2000
Trowbridge, C.W. “Three-Dimensional Analysis” The Electrical Engineering Handbook Ed. Richard C. Dorf Boca Raton: CRC Press LLC, 2000
44 Three-Dimensional Analy SIS 44.1 Introduction 44.2 The Field equations Low Frequency Fields. Statics Limit 44.3 Numerical Methods C. W. Trowbridge Finite Elements. Edge Elements.Integral Methods Vector Fields. Inc 44.4 Modern Design Environment 44.1 Introduction The three-dimensional analysis of electromagnetic fields the use of numerical techniques exploiting he 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 Maxwells 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 (p)and electric current density (). The field vectors are not independent since they are further related by the material constitutive properties: B=HH,D=EE, and j= oE where u, E, and o 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 Because the flux density vector satisfies a zero divergence condition(div B=0), it can be expressed in terms f a magnetic vector potential A, ie,B=curl A, and it follows from Faraday s law that E=-(aA/at+ vv), 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 V. A=-(uoV+ HEaV/at)(Lorentz gauge)is selected, then the field equations in terms of A and V are c 2000 by CRC Press LLC
© 2000 by CRC Press LLC 44 Three-Dimensional Analysis 44.1 Introduction 44.2 The Field Equations Low Frequency Fields • Statics Limit 44.3 Numerical Methods Finite Elements • Edge Elements • Integral Methods 44.4 Modern 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
VxV×A+σ V.A + uo where o 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 o, i. e, A=A. exp(jot), 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 Maxwells 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 o =0 and curl H=Js Eqs. (44. 1)can be eplaced by V2i=0, where t is a scalar potential defined by H=-vd. The scalar and vector field regions are oupled 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 V·μVd=V·pH (44.2) where d is known as the reduced magnetic scalar potential with H= H,-va, and H, the source field given by the biot Savart law. Some needed in solving Eq.(44.2)numerically, in practice, as H, will often be calculated to a higher accuracy than 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 and H,, depending upon how the computations are carried out. One approach that avoids this difficulty is to use the total scalar potential y in regions that have zero current density [Simkin and Trowbridge, 1979, i.e., where H=-Vo and H is the coercive field for the material ap satisfies v·μψ=V·μH 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 φ=y+|Hdr (444) where I is the contour delineating the two regions that must not intersect a current-carrying region; otherwise the definition of v will be violated c 2000 by CRC Press LLC
© 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 = Ao 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 approximation 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 = Js 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 ¹ • mHs (44.2) where f is known as the reduced magnetic scalar potential with H = Hs – ¹f, and Hs the source field given by the Biot Savart law. Some care is needed in solving Eq. (44.2) numerically, in practice, as Hs 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 Hs , 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 Hc is the coercive field for the material c satisfies ¹ • m¹c 5 ¹ • mHc (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. —¥ —¥ + + = — —× È Î Í Í ˘ ˚ ˙ ˙ + = —×— 1 1 2 2 2 2 m s ¶ ¶ ¶ ¶ m m ¶ ¶ ms ¶ ¶ A A A A t t V t V t V e e f y = + Ú Hs dG G 0
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, anisotropi and even hysteretic effects. Their use for complete three-dimensional models is not so widespread because of ne 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 regio (elements)over which the solution is assumed to follow a simple local approximating trial function(shap 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 a, a,+ a3y+a,2+asxy a yz+ a,zx t a NU;(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, are called the local shape functions and the parameters U, are lements 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 Galerkins method [Zienkiewicz, 1990)). At the element level, the residual R, is given by VNμVNd92U1+|NQd (44.6) clem 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 ancy 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. uaU/an)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 equired can be prohibitive because of the vector nature of magnetization. The magnetic material must be haracterized 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 c 2000 by CRC Press LLC
© 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 nonlinear 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 Ni are called the local shape functions and the parameters Ui 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 Ri 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 redundancy 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. y ª u = a1 2 + a x + a3 y + a4z + a5 xy + a6 yz + a7zx + a8 xyz = Â Ni Ui R N N d U N Qd i = — i — j j i È Î Í Í ˘ ˚ ˙ ˙ + Ú Ú m W W elem elem
FIGURE 44.1 Three-dimensional second-order Isoparametric finite elements. (a) Left, master tetrahedron, 10 nodes in al space(E, n, 4); right, actual tetrahedron, 10 nodes in global space (x, y, z).(b)Left, master hexahedron, 20 nodes in local space(S, n, S); right, actual hexahedron, 10 nodes in global space(x, y, z) 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, sen tangential field continuity is enforced [Bossavit, 1988]. The field equations [Eq.(44.1)] in terms of the ield intensity for the low frequency limit reduce to V×V×H+aH 0 (44.7) t c 2000 by CRC Press LLC
© 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 highfrequency 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