Differences

This shows you the differences between two versions of the page.

--- users_manual [2020/10/19 00:16]
stanzurek [6.1 Finite Element Formulation]
+++ users_manual [2021/11/22 21:25] (current)
stanzurek [A.1 Modeling Permanent Magnets]
@@ Line 3133: / Line 3133: @@
 ==== 6.2 Linear Solvers ====
-For all problems, variations of the iterative Conjugate Gradient solver are used. This technique is
+For all problems, variations of the iterative Conjugate Gradient solver are used. This technique is appropriate for the sort of problem that FEMM solves, because the matrices are symmetric and very sparse. A row-based storage scheme is used in which only the nonzero elements of the diagonal and upper triangular part of the matrix are solved.
-appropriate for the sort of problem that FEMMsolves, because thematrices are symmetric and very
-sparse. A row-based storage scheme is used in which only the nonzero elements of the diagonal
+For magnetostatic and electrostatic problems, the preconditioned conjugate gradient (PCG) code is based on the discussion in [11]. Minor modifications are made to this algorithm to avoid computing certain quantities more than once per iteration. Although Silvester promotes the use of the Incomplete Cholesky preconditioner, it is not used in FEMM, because it nearly doubles the storage requirements–for each element of the matrix stored, a corresponding element of the preconditioner must also be stored. Instead, the Symmetric Successive Over-Relaxation (SSOR) preconditioner, as described in [13], is used. The advantage of this preconditioner is that it is built on the fly in a simple way using only the matrix elements that are already in storage. In general,
-and upper triangular part of the matrix are solved.
+the speed of PCG using SSOR is said to be comparable to the speed of PCG with Incomplete Cholesky.
-For magnetostatic and electrostatic problems, the preconditioned conjugate gradient (PCG)
-code is based on the discussion in [11]. Minor modifications are made to this algorithm to avoid
+For harmonic problems, the regular PCG algorithm cannot be used; the matrix that arises in the formulation of harmonic problems is Complex Symmetric (i.e. //A = A<sup>T</sup>// ), rather than Hermitian (i.e. //A = A∗//). Curiously, there is not much literature available on iterative solvers for complex symmetric problems, given the number of diverse applications in which these problems arise.
-computing certain quantities more than once per iteration. Although Silvester promotes the use
-of the Incomplete Cholesky preconditioner, it is not used in FEMM, because it nearly doubles
+However, there is a very good paper on the solution of linear problems with complex symmetric matrices via various flavors of Conjugate Gradient by Freund [14]. The techniques discussed by Freund allow one to operate directly on the complex symmetric matrix and take advantage of the symmetric structure to minimize the number of computations that must be performed per iteration. Although Freund supports Quasi-MinimumResidual approach, FEMMuses the complex symmetric version of biconjugate gradient also described in [14]. After coding and comparing the the speed of both BCG and QMR, it was found that BCG is somewhat faster due to a relatively smaller number of computations that must be performed per iteration (even though QMR has better convergence properties than BCG).
-the storage requirements–for each element of the matrix stored, a corresponding element of the
+However, using the algorithms as described by [14], solution times were unacceptably long. To decrease solution times, the complex symmetric BCG algorithm was modified to include the SSOR preconditioner (built in exactly the same way as for magnetostatic problems). Including the SSOR preconditioner in complex symmetric BCG problems usually yields an order of magnitude improvement in speed over no preconditioner.
-preconditioner must also be stored. Instead, the Symmetric Successive Over-Relaxation (SSOR)
-preconditioner, as described in [13], is used. The advantage of this preconditioner is that it is built
+In all problems, a node renumbering scheme is used. Although the conjugate gradient schemes work well without renumbering, the renumbering seems to roughly halve the solution time. There is an overall advantage to using the renumbering, because the of time required to perform the renumbering is small compared to the time required to run CG or BCG. Although there are many possible approaches to renumbering, FEMM uses the Cuthill-McKee method as described in [2]. Although there are newer schemes that yield a tighter profile, Cuthill-McKee does a relatively good job and requires very little to execute. The renumbering code is a hold-over from an early version of FEMM that employed a banded Gauss Elimination solver in which a good node numbering is essential to good performance. The renumbering speeds up CG and BCG by reducing the error between the SSOR approximation of //A<sup>-1</sup>// and the exact //A<sup>-1</sup>//. An interesting paper on the effect of the ordering of the unknowns on convergence in conjugate gradient methods is [15].
-on the fly in a simple way using only the matrix elements that are already in storage. In general,
-the speed of PCG using SSOR is said to be comparable to the speed of PCG with Incomplete
-Cholesky.
-For harmonic problems, the regular PCG algorithm cannot be used; the matrix that arises in
-the formulation of harmonic problems is Complex Symmetric (i.e. A = AT ), rather than Hermitian
-(i.e. A = A∗). Curiously, there is not much literature available on iterative solvers for complex
-symmetric problems, given the number of diverse applications in which these problems arise.
-However, there is a very good paper on the solution of linear problems with complex symmetric
-matrices via various flavors of Conjugate Gradient by Freund [14]. The techniques discussed
-by Freund allow one to operate directly on the complex symmetric matrix and take advantage
-of the symmetric structure to minimize the number of computations that must be performed per
-iteration. Although Freund supports Quasi-MinimumResidual approach, FEMMuses the complex
-symmetric version of biconjugate gradient also described in [14]. After coding and comparing the
-the speed of both BCG and QMR, it was found that BCG is somewhat faster due to a relatively
-smaller number of computations that must be performed per iteration (even though QMR has better
-convergence properties than BCG).
-However, using the algorithms as described by [14], solution times were unacceptably long.
-To decrease solution times, the complex symmetric BCG algorithm was modified to include the
-SSOR preconditioner (built in exactly the same way as for magnetostatic problems). Including the
-SSOR preconditioner in complex symmetric BCG problems usually yields an order of magnitude
-improvement in speed over no preconditioner.
-In all problems, a node renumbering scheme is used. Although the conjugate gradient schemes
-work well without renumbering, the renumbering seems to roughly halve the solution time. There
-is an overall advantage to using the renumbering, because the of time required to perform the
-renumbering is small compared to the time required to run CG or BCG. Although there are many
-possible approaches to renumbering, FEMM uses the Cuthill-McKee method as described in [2].
-Although there are newer schemes that yield a tighter profile, Cuthill-McKee does a relatively good
-job and requires very little to execute. The renumbering code is a hold-over from an early version
-of FEMM that employed a banded Gauss Elimination solver in which a good node numbering is
-essential to good performance. The renumbering speeds up CG and BCG by reducing the error
-between the SSOR approximation of A−1 and the exact A−1. An interesting paper on the effect of
-the ordering of the unknowns on convergence in conjugate gradient methods is [15].
 ==== 6.3 Field Smoothing ====
-Since first-order triangles are used by FEMM, the resulting solution for B and H obtained by
+Since first-order triangles are used by FEMM, the resulting solution for //B// and //H// obtained by differentiating //A// is constant over each element. If the raw //B// and //H// are used by the postprocessor, density plots of //B// and 2-D plots of field quantities along user-defined contours look terrible. The values of //B// and //H// also are not so accurate at points in an element away from the element’s centroid.
-differentiating A is constant over each element. If the raw B and H are used by the postprocessor,
-density plots of B and 2-D plots of field quantities along user-defined contours look terrible. The
+The use of smoothing to recover the accuracy lost by differentiating //A// is known as superconvergence. Of the greatest interest to FEMM are so-called “patch recovery” techniques. The basic idea is the the solutions for //B// are most accurate at the centroid of the triangular element (known as its //Gauss Point//). One desires a continuous profile of //B// that can be interpolated from nodal values, in the same way that vector potential //A// can be represented. The problem is, the “raw” solution of //B// is multivalued at any node point, those values being the different constant values of //B// in each element surrounding the node point. The general approach to estimating the “true” value of //B// at any node point is to fit a least-squares plane through the values of //B// at the Gauss points of all elements that surround a node of interest, and to take the value of the plane at the node point’s location as its smoothed value of //B// [16].
-values of B and H also are not so accurate at points in an element away from the element’s centroid.
-The use of smoothing to recover the accuracy lost by differentiating A is known as superconvergence.
+However, this approach to patch recovery has a lot of shortcomings. For the rather irregular meshes that can arise in finite elements, the least-squares fit problem can be ill-condition, or even singular, at some nodes in the finite element mesh. Furthermore, the superconvergence solution can actually be less accurate than the piece-wise constant solution in the neighborhood of boundaries and interfaces.
-Of the greatest interest to FEMM are so-called “patch recovery” techniques. The basic
-idea is the the solutions for B are most accurate at the centroid of the triangular element (known as
+One can note that the patch recovery method is merely a weighted average of the flux densities in all of the elements surrounding a given node. Instead of a least-squares fit, FEMM simply weights the values of flux density in each adjacent element’s Gauss point with a value inversely proportional to the distance from the Gauss point to the node point of interest. Away from boundaries, the results seem to be nearly as good as a least-squares fit. At boundaries and interfaces, the smoothed solution is no worse than the unsmoothed solution.
-its Gauss Point). One desires a continuous profile of B that can be interpolated from nodal values,
+A related approach is used for smoothing //D// and //E// in electrostatics problems. In electrostatics problems, however, nodal values of //D// are found by taking the gradient of a best-fit plane through the voltage of the neighboring nodes. A number of special cases must be caught so that reasonable results are obtained at various boundaries and surfaces.
-in the same way that vector potential A can be represented. The problem is, the “raw” solution of
-B is multivalued at any node point, those values being the different constant values of B in each
-element surrounding the node point. The general approach to estimating the “true” value of B at
-any node point is to fit a least-squares plane through the values of B at the Gauss points of all
-elements that surround a node of interest, and to take the value of the plane at the node point’s
-location as its smoothed value of B [16].
-However, this approach to patch recovery has a lot of shortcomings. For the rather irregular
-meshes that can arise in finite elements, the least-squares fit problem can be ill-condition, or even
-singular, at some nodes in the finite elementmesh. Furthermore, the superconvergence solution can
-actually be less accurate than the piece-wise constant solution in the neighborhood of boundaries
-and interfaces.
-One can note that the patch recovery method is merely a weighted average of the flux densities
-in all of the elements surrounding a given node. Instead of a least-squares fit, FEMM simply
-weights the values of flux density in each adjacent element’s Gauss point with a value inversely
-proportional to the distance from the Gauss point to the node point of interest. Away from boundaries,
-the results seem to be nearly as good as a least-squares fit. At boundaries and interfaces, the
-smoothed solution is no worse than the unsmoothed solution.
-A related approach is used for smoothing D and E in electrostatics problems. In electrostatics
-problems, however, nodal values of D are found by taking the gradient of a best-fit plane through
-the voltage of the neighboring nodes. A number of special cases must be caught so that reasonable
-results are obtained at various boundaries and surfaces.
 ===== Bibliography =====
@@ Line 3260: / Line 3206: @@
 ==== Appendix ====
 ==== A.1 Modeling Permanent Magnets ====
-FEMM accommodates permanent magnets, but there are some special rules associated with properly
+FEMM accommodates permanent magnets, but there are some special rules associated with properly modeling them. This appendix will explain how to distill enough information from a manufacturer’s literature to properly define the material in FEMM.
-modeling them. This appendix will explain how to distill enough information from a manufacturer’s
-literature to properly define the material in FEMM.
+The manufacturer provides information about their material in the form of a demagnetization curve. A sample curve for Alnico 5 is pictured in Figure A.1. The task is to get the appropriate information out of the curve put in a FEMM Block Property model.
-The manufacturer provides information about their material in the form of a demagnetization
-curve. A sample curve for Alnico 5 is pictured in Figure A.1. The task is to get the appropriate
+{{fig_a-1.png}}
-information out of the curve put in a FEMM Block Property model.
-Magnets can be modeled from several different, but equally valid, points of view. From the
-perspective finite element analysis, the most useful model is to think of the magnet as a volume of
-ferromagnetic material surrounded by a thin sheet of current, as shown in Figure A.2. From this
 Figure A.1: Sample demagnetization curve for Alnico 5
+Magnets can be modeled from several different, but equally valid, points of view. From the perspective finite element analysis, the most useful model is to think of the magnet as a volume of ferromagnetic material surrounded by a thin sheet of current, as shown in Figure A.2. From this point of view, the demagnetization curve is what occurs when different amounts of magnetomotive force are applied to a long magnet, acting in the direction opposing the field of the magnet. When enough MMF is applied so that the field is exactly cancelled out, the applied MMF must be exactly the same as the MMF that is driving the magnet. The B-H profile that is traversed on the way to the //B// = 0 point is just the B-H curve of the material inside the magnet.
+{{fig_a-2.png}}
 Figure A.2: Magnet as an equivalent current sheet.
-point of view, the demagnetization curve is what occurs when different amounts of magnetomotive
-force are applied to a long magnet, acting in the direction opposing the field of the magnet. When
-enoughMMF is applied so that the field is exactly cancelled out, the appliedMMF must be exactly
+Using these insights, the permanent magnet can be modeled. The //coercivity// (denoted //Hc//) of the magnet is the absolute value of the MMF that it takes to bring the the field in the magnet to zero. This value (in units of Amps/Meter) is entered in the ''H_c'' box in the Block Property dialog (see Figure 2.9). If the magnet material is nonlinear, the appropriate values to enter in the B-H data dialog can be obtained by shifting the curve to the right by exactly //Hc//, so that the //B// = 0 point lines up with the origin. For example, the shifted demagnetization curve corresponding to Alnico 5 is pictured in Figure A.3. If the demagnetization curve is straight enough to be considered linear, one can obtain the appropriate permeability by taking the slope of the demagnetization curve.
-the same as the MMF that is driving the magnet. The B-H profile that is traversed on the way to
-the B = 0 point is just the B-H curve of the material inside the magnet.
+Strong rare-earth materials at room temperature have a very linear demagnetization curve. Usually, a linear model is sufficient for these materials. In addition, these materials have a relative permeability very close to 1. The modeling of these materials can be simplified (while only incurring small errors) by assuming that the permeability is exactly 1. Then, if you know the energy product of the magnet material in units of MGOe (the unit in which the energy product is almost always given), the appropriate //Hc// can be determined via (A.1).
-Using these insights, the permanent magnet can be modeled. The coercivity (denoted Hc) of
-the magnet is the absolute value of the MMF that it takes to bring the the field in the magnet to
+(A.1)  $  Hc = \frac{5·10^5· \sqrt{E}}{π}  $
-zero. This value (in units of Amps/Meter) is entered in the H_c box in the Block Property dialog
-(see Figure 2.9). If the magnet material is nonlinear, the appropriate values to enter in the B-H data
+where //E// is the energy product in MGOe and the resulting //Hc// is in units of A/m (e.g. 40 MGOe ≈ 10<sup>6</sup> A/m).
-dialog can be obtained by shifting the curve to the right by exactly Hc, so that the B = 0 point lines
-up with the origin. For example, the shifted demagnetization curve corresponding to Alnico 5 is
+{{fig_a-3.png}}
-pictured in Figure A.3. If the demagnetization curve is straight enough to be considered linear, one
-can obtain the appropriate permeability by taking the slope of the demagnetization curve.
-Strong rare-earth materials at room temperature have a very linear demagnetization curve. Usually,
-a linear model is sufficient for these materials. In addition, these materials have a relative permeability
-very close to 1. The modeling of these materials can be simplified (while only incurring
-small errors) by assuming that the permeability is exactly 1. Then, if you know the energy product
-of the magnet material in units of MGOe (the unit in which the energy product is almost always
-given), the appropriate Hc can be determined via (A.1).
-Hc =
-(105)√E
-p
-(A.1)
-where E is the energy product in MGOe and the resulting Hc is in units of A/m (e.g. 40 MGOe
-≈ 106 A/m).
-With alnico magnets, great care must be taken in interpreting the finite element results. Unlike
-rare-earth magnets, these magnets exhibit a great degree of hysteresis when they are demagnetized.
-That is, when the flux density is pushed below the “knee” in the demagnetization curve, the
-flux level does not recover to the previous magnitude when the opposing MMF is removed. This
-hysteresis is illustrated in Figure A.4. This sort of demagntization and recoil can occur when the
 Figure A.3: Shifting the B-H curve of a permanent magnet
+With alnico magnets, great care must be taken in interpreting the finite element results. Unlike rare-earth magnets, these magnets exhibit a great degree of hysteresis when they are demagnetized. That is, when the flux density is pushed below the “knee” in the demagnetization curve, the flux level does not recover to the previous magnitude when the opposing MMF is removed. This hysteresis is illustrated in Figure A.4. This sort of demagntization and recoil can occur when the magnets are being handled prior to assembly into a device. In a motor, the magnets will demagnetize somewhat when the motor is first started. They will eventually end up running back and forth along a recoil line that is below the “virgin” demagnetization curve. The point is that the modeler cannot be sure exactly where the magnets are operating–an analysis that takes this sort of hysteresis into account is beyond the scope of FEMM. Note, however, that this caution applies only for nonlinear magnets; for practical purposes, rare-earth magnets generally do not exhibit this sort of hysteresis behavior.
+{{fig_a-4.png}}
 Figure A.4: Recoil in partially demagnetized Alnico 5.
-magnets are being handled prior to assembly into a device. In a motor, the magnets will demagnetize
-somewhat when the motor is first started. They will eventually end up running back and
-forth along a recoil line that is below the “virgin” demagnetization curve. The point is that the
-modeler cannot be sure exactly where the magnets are operating–an analysis that takes this sort
-of hysteresis into account is beyond the scope of FEMM. Note, however, that this caution applies
-only for nonlinear magnets; for practical purposes, rare-earth magnets generally do not exhibit this
-sort of hysteresis behavior.
-===== A.2 Bulk Lamination Modeling =====
+==== A.2 Bulk Lamination Modeling ====
-A great number ofmagnetic devices employ cores built up out of thin laminations for the purpose of
-reducing eddy current effects. One way to model these materials within a finite element framework
+A great number of magnetic devices employ cores built up out of thin laminations for the purpose of reducing eddy current effects. One way to model these materials within a finite element framework would be to model each discrete lamination (and the insulation between laminations) in the finite element geometry. An alternative is to treat the laminated material as a continuum and derive bulk properties that yield essentially the same results, while requiring a much less elaborate finite element mesh. FEMM has implemented this bulk approach to laminations.
-would be to model each discrete lamination (and the insulation between laminations) in the finite
-element geometry. An alternative is to treat the laminated material as a continuum and derive
+Consider that the flux can flow through the lamination in a combination of two ways: via the “easy” direction down the laminations, or the “hard” way, across the thickness of the laminations. The hard direction is difficult for flux for two reasons. First, the rolling process makes the iron somewhat less permeable than in the easy direction. Second, and most importantly, the flux must traverse the insulation between laminations, which typically has a unit permeability.
-bulk properties that yield essentially the same results, while requiring a much less elaborate finite
-element mesh. FEMM has implemented this bulk approach to laminations.
+The first assumption in deriving the bulk permeability model is that the permeability in the iron itself is isotropic. This isn’t quite true, but almost all of the reluctance in the hard direction results from crossing the gap between laminations. Having a significant error in the hard direction permeability in the iron itself only results in a trivial change in the bulk reluctance in the cross-lamination direction.
-Consider that the flux can flow through the lamination in a combination of two ways: via the
-“easy” direction down the laminations, or the “hard” way, across the thickness of the laminations.
-The hard direction is difficult for flux for two reasons. First, the rolling process makes the iron
-somewhat less permeable than in the easy direction. Second, and most importantly, the flux must
-traverse the insulation between laminations, which typically has a unit permeability.
-The first assumption in deriving the bulk permeability model is that the permeability in the
-iron itself is isotropic. This isn’t quite true, but almost all of the reluctance in the hard direction
-results from crossing the gap between laminations. Having a significant error in the hard direction
-permeability in the iron itself only results in a trivial change in the bulk reluctance in the crosslamination
-direction.
 Armed with this assumption, a circuit model can be produced for each direction of flux travel.
 For the easy direction, the circuit model is pictured in Figure A.5. There are two reluctances in
 parallel–one for flux that flows through the iron part of the laminations:
-Rez, f e =
-L
+(A.2) $ R_{ez,fe} = \frac{L}{μ_r μ_0 c W} $
-μrμocW
-(A.2)
 and another reluctance for flux that flows through the air between laminations:
-Rez,air =
-L
+(A.3) $ R_{ez,air} = \frac{L}{μ_r (1 - c) W} $
-μo(1−c)W
-(A.3)
+where //L// and //W// are the length and width of the path traversed, and //c// is the fraction of the path filled with iron. Adding these two reluctances in parallel yields:
-where L andW are the length and width of the path traversed, and c is the fraction of the path filled
-with iron. Adding these two reluctances in parallel yields:
+(A.4) $ R_{ez,fe} = \frac{L}{((1 - c) + c μ_r ) μ_0 W} $
-Rez =
-L
+Since //L// and //W// are arbitrarily chosen, the bulk permeability of the section is:
-((1−c)+cμr)μoW
-(A.4)
+(A.5)  $ μez = ((1−c) +c μ_r) μ_0 $
-Since L and W are arbitrarily chosen, the bulk permeability of the section is:
-μez = ((1−c)+cμr)μo (A.5)
+{{fig_a-5.png}}
 Figure A.5: Equivalent circuit for flux in the “easy” direction
-For the solution of nonlinear problems, the derivation of the changes to Newton’s method to accomodate
-the bulk lamination model are greatly simplified if it is assumed that (1−c) << cμr. In
+For the solution of nonlinear problems, the derivation of the changes to Newton’s method to accomodate the bulk lamination model are greatly simplified if it is assumed that $(1−c) << c μ_r$. In this case, $μ_{ez}$ can be approximated as:
-this case, μez can be approximated as:
-μez ≈ cμrμo (A.6)
+(A.6)  $ μ_{ez} ≈ c μ_r μ_0 $
-This approximation leads to only trivial errors until the fill factor approaches zero. For example, if
-μr = 1000 with a 90% fill, the difference between (A.5) and (A.6) is only about 0.01%.
+This approximation leads to only trivial errors until the fill factor approaches zero. For example, if //μr// = 1000 with a 90% fill, the difference between (A.5) and (A.6) is only about 0.01%.
-For the hard direction, a different equivalent circuit, pictured in Figure A.6 can be drawn. In
-this case, the circuit is two reluctances in series, as the flux has to cross the insulation and the
+For the hard direction, a different equivalent circuit, pictured in Figure A.6 can be drawn. In this case, the circuit is two reluctances in series, as the flux has to cross the insulation and the lamination in succession. These reluctances are:
-lamination in succession. These reluctances are:
-Rhard, f e =
+(A.7)  $ R_{hard,fe} = \frac{c L}{μ_r μ_0 W} $
-cL
-μrμoW
+(A.8) $ R_{hard,air} = \frac{(1−c) L}{μ_0 W} $
-(A.7)
-Rhard,air =
-(1−c)L
-μoW
-(A.8)
 Adding these two reluctances together in series yields:
-Rhard =
+(A.9) $ Rhard = \frac{c+(1−c)μ_r} {μ_r μ_0} \frac{L}{W} $
-c+(1−c)μr
-μrμo
+Since //L// and //W// are arbitrary, the bulk permeability in the hard direction is:
-L
+(A.10)  $ μ_{hard} = \frac{μ_r μ_0}{c+(1−c) μ_r} $
-W
-(A.9)
+{{fig_a-5.png}}
 Figure A.6: Equivalent circuit for flux in the “hard” direction.
-Since L and W are arbitrary, the bulk permeability in the hard direction is:
-μhard =
-μrμo
+If the material is laminated “in-plane,” all flux is flowing in the easy direction, and (A.6) is used as the permeability for each element. In problems that are laminated parallel to x or y, (A.6) and (A.10) are used as permeabilities in the standard fashion for elements with an anisotropic permeability.
-c+(1−c)μr
-(A.10)
+For harmonic problems, eddy currents flow in the laminations, and hysteresis causes additional loss. If the laminations are thin compared to the other dimensions of the geometry, the effects of eddy currents and hysteresis can be encapsulated in a frequency-dependent permeability [6]. In this case, the magnetostatic permeability, $μ_r$ is simply replaced by the frequency-dependent permeability $μ_{fd}$ in (A.6) and (A.10):
-If the material is laminated “in-plane,” all flux is flowing in the easy direction, and (A.6) is
-used as the permeability for each element. In problems that are laminated parallel to x or y,
+(A.11)  $ μ_{fd} = \frac { μ_{r} e^\frac{−jφ_h}{2} tanh \left[e^\frac{−jφ_h}{2} \sqrt{j ω σ μ_r μ_0} \frac{d}{2} \right] } {\sqrt{j ω σ μ_r μ_0} \frac{d}{2}} $
-(A.6) and (A.10) are used as permeabilities in the standard fashion for elements with an anisotropic
-permeability.
+In (A.11), $φ_h$ represents a constant phase lag between B and H due to hysteresis, $σ$ is the conductivity of the lamination material, $d$ is the thickness of the iron part of the lamination, and $ω$ is the frequency of excitation in rad/s. Note that the concept of hysteresis-induced lag can be applied to non-laminated materials as well, simply by multiplying the magnetostatic permeability by e^{−jφ_h} for harmonic problems.
-For harmonic problems, eddy currents flow in the laminations, and hysteresis causes additional
-loss. If the laminations are thin compared to the other dimensions of the geometry, the effects
-of eddy currents and hysteresis can be encapsulated in a frequency-dependent permeability [6].
-In this case, the magnetostatic permeability, μr is simply replaced by the frequency-dependent
-permeability μf d in (A.6) and (A.10):
-μf d =
-μre−jfh
-tanh
-h
-e−jfh
- √jwsμrμo
-d
-i
-√jwsμrμo
-d
-(A.11)
-In (A.11), fh represents a constant phase lag between B and H due to hysteresis, s is the conductivity
-of the lamination material, d is the thickness of the iron part of the lamination, and w is the
-frequency of excitation in rad/s. Note that the concept of hysteresis-induced lag can be applied to
-non-laminated materials as well, simply by multiplying the magnetostatic permeability by e−jfh
-for harmonic problems.
 ==== A.3 Open Boundary Problems ====
-Typically, finite element methods are best suited to problems with well-defined, closed solution
+Typically, finite element methods are best suited to problems with well-defined, closed solution regions. However, a large number of problems that one might like to address have no natural outer boundary. A prime example is a solenoid in air. The boundary condition that one would //like// to apply is $A = 0$ at $r = ∞$. However, finite element methods, by nature, imply a finite domain. Fortunately, there are methods that can be applied to get solutions that closely approximate the “open boundary” solution using finite element methods.
-regions. However, a large number of problems that one might like to address have no natural
-outer boundary. A prime example is a solenoid in air. The boundary condition that one would like
-to apply is A = 0 at r = ¥. However, finite element methods, by nature, imply a finite domain.
-Fortunately, there are methods that can be applied to get solutions that closely approximate the
-“open boundary” solution using finite element methods.
 === A.3.1 Truncation of Outer Boundaries ===
-The simplest, but least accurate, way to proceed is to pick an arbitrary boundary “far enough”
-away from the area of interest and declare either A = 0 or ¶A/¶n = 0 on this boundary. According
+The simplest, but least accurate, way to proceed is to pick an arbitrary boundary “far enough” away from the area of interest and declare either $A = 0$ or $ \partial A / \partial n = 0$ on this boundary. According to [17], a rule of thumb is that the distance from the center of the problem to the outer boundary should be at least five times the distance from the center to the outside of the objects of interest. Truncation is the method employed by most magnetics finite element programs, because it requires no additional effort to implement.
-to [17], a rule of thumb is that the distance from the center of the problem to the outer boundary
-should be at least five times the distance from the center to the outside of the objects of interest.
+The down side to truncation is that get an accurate solution in the region of interest, a volume of air much larger than the region of interest must also be modeled. Usually, this large region exterior to the area of interest can be modeled with a relatively coarse mesh to keep solution times to a minimum. However, some extra time and space is still required to solve for a region in which one has little interest.
-Truncation is the method employed by most magnetics finite element programs, because it requires
-no additional effort to implement.
-The down side to truncation is that get an accurate solution in the region of interest, a volume
-of air much larger than the region of interest must also be modeled. Usually, this large region
-exterior to the area of interest can be modeled with a relatively coarse mesh to keep solution times
-to a minimum. However, some extra time and space is still required to solve for a region in which
-one has little interest.
 === A.3.2 Asymptotic Boundary Conditions ===
-A thorough review of open boundary techniques is contained in [17]. Perhaps the simple way to
-approximate an “open” boundary (other than truncation) described in [17] is to use asymptotic
+A thorough review of open boundary techniques is contained in [17]. Perhaps the simple way to approximate an “open” boundary (other than truncation) described in [17] is to use asymptotic boundary conditions. The result is that by carefully specifying the parameters for the “mixed” boundary condition, and then applying this boundary condition to a circular outer boundary, the unbounded solution can be closely approximated. An example that employs an asymptotic boundary condition to obtain an unbounded field solution is the ''axi1.fem'' example included in the distribution.
-boundary conditions. The result is that by carefully specifying the parameters for the “mixed”
-boundary condition, and then applying this boundary condition to a circular outer boundary, the
+Consider a 2-D planar problem in polar coordinates. The domain is a circular shell of radius //ro// in an unbounded region. As //r → ∞//, vector potential //A// goes to zero. On the surface of the circle, the vector is a prescribed function of //θ//. This problem has an analytical solution, which is:
-unbounded solution can be closely approximated. An example that employs an asymptotic boundary
-condition to obtain an unbounded field solution is the axi1.fem example included in the distribution.
+(A.12)  $ A(r,θ) = \sum_{m=1}^{∞} \frac{a_m}{r^m} cos(mθ + α_m) $
-Consider a 2-D planar problem in polar coordinates. The domain is a circular shell of radius ro
-in an unbounded region. As r →¥, vector potential A goes to zero. On the surface of the circle,
+where the $a_m$ and $α_m$ parameters are chosen so that the solution matches the prescribed potential on the surface of the circle.
-the vector is a prescribed function of q. This problem has an analytical solution, which is:
-A(r,q) =
+One could think of this solution as describing the solution exterior to a finite element problem with a circular outer boundary. The solution is described inside the circle via a finite element solution. The trick is to knit together the analytical solution outside the circle to the finite element solution inside the circle.
-¥å
-m=1
+From inspecting (A.12), one can see that the higher-numbered harmonic, the faster the magnitude of the harmonic decays with respect to increasing r. After only a short distance, the higher numbered harmonics decay to the extent that almost all of the open-space solution is described by only the leading harmonic. If n is the number of the leading harmonic, the open-field solution for large, but not infinite, r is closely described by:
-am
-rm cos(mq+am) (A.12)
+(A.13) $ A(r,θ) ≈ \frac{a_n}{r^n} cos(nθ + α_n) $
-where the am and am parameters are chosen so that the solution matches the prescribed potential
-on the surface of the circle.
-One could think of this solution as describing the solution exterior to a finite element problem
-with a circular outer boundary. The solution is described inside the circle via a finite element
-solution. The trick is to knit together the analytical solution outside the circle to the finite element
-solution inside the circle.
-From inspecting (A.12), one can see that the higher-numbered harmonic, the faster the magnitude
-of the harmonic decays with respect to increasing r. After only a short distance, the highernumbered
-harmonics decay to the extent that almost all of the open-space solution is described by
-only the leading harmonic. If n is the number of the leading harmonic, the open-field solution for
-large, but not infinite, r is closely described by:
-A(r,q) ≈
-an
-rn cos(nq+an) (A.13)
 Differentiating with respect to r yields:
-¶A
-¶r
+(A.14) $ \frac{ \partial A }{\partial r}  = - \frac{n a_n}{r^{n+1}} cos(nθ + α_n) $
-= −
-nan
-rn+1 cos(nq+an) (A.14)
 If (A.14) is solved for an and substituted into (A.13), the result is:
-¶A
-¶r
+(A.15) $ \frac{ \partial A }{\partial r} + \left( \frac{n}{r} \right) A  = 0 $
-+
-n
+Now, (A.15) is a very useful result. This is the same form as the “mixed” boundary condition supported by FEMM. If the outer edge of the solution domain is circular, and the outer finite element boundary is somewhat removed from the area of primary interest, the open domain solution can be closely approximated by applying (A.15) the circular boundary.
-r
+To apply the Asymptotic Boundary Condition, define a new, mixed-type boundary condition. Then, pick the parameters so that:
-A = 0 (A.15)
-Now, (A.15) is a very useful result. This is the same form as the “mixed” boundary condition
+(A.16)  $ c_0 = \frac{n}{μ_0 r_0} $
-supported by FEMM. If the outer edge of the solution domain is circular, and the outer finite element
-boundary is somewhat removed from the area of primary interest, the open domain solution
+(A.17)  $ c_1 = 0 $
-can be closely approximated by applying (A.15) the circular boundary.
-To apply the Asymptotic Boundary Condition, define a new, mixed-type boundary condition.
+where $r_0$ is the outer radius of the region in meters (regardless of the working length units), and $μ_o = 4 π 10^{−7}$.
-Then, pick the parameters so that:
-c0 =
+Although the above derivation was specifically for 2-D problems, it turns out that when the same derivation is done for the axisymmetric case, the definition of the mixed boundary condition coefficients are exactly the same as (A.16).
-n
-μoro
+To apply the Asymptotic Boundary Condition to electrostatics problems, pick the parameters so that:
-(A.16)
-c1 = 0 (A.17)
+(A.18)  $ c_0 = \frac{ε_0 n}{r_0} $
-where ro is the outer radius of the region in meters (regardless of the working length units), and
-μo = 4p(10−7).
+$ c_1 = 0 $
-Although the above derivation was specifically for 2-D problems, it turns out that when the
-same derivation is done for the axisymmetric case, the definition of the mixed boundary condition
+where $r_0$ is the outer radius of the region in meters (regardless of the working length units), and ε0 = 8.85418781762e-012. Note that ε0 is defined in the Lua implementation in both the pre- and post-processors as the global variable ''eo'', which can be used in any script or edit box in the program.
-coefficients are exactly the same as (A.16).
-To apply the Asymptotic Boundary Condition to electrostatics problems, pick the parameters
+{{fig_a-7.png}}
-so that:
-c0 =
-eon
-ro
-(A.18)
-c1 = 0
-where ro is the outer radius of the region in meters (regardless of the working length units), and
-eo = 8.85418781762e-012. Note that eo is defined in the Lua implementation in both the preand
-post-processors as the global variable eo, which can be used in any script or edit box in the
-program.
 Figure A.7: Air-cored coil with “open” boundary condition
-Just like magnetics problems, it turns out that 2-D problems are also described by (A.18). One
-subtle difference, however, is that n =1 in the axisymmetric case corresponds to the case in which
+Just like magnetics problems, it turns out that 2-D problems are also described by (A.18). One subtle difference, however, is that n =1 in the axisymmetric case corresponds to the case in which there is a net charge with (i.e. the geometry looks like a point charge when viewed from a distance), whereas the n= 1 case corresponds to a dipole charge distribution in 2D planar problems. If charge is conserved in the geometry of interest in the axisymmetric case, one needs to use n =2 in Eq. (A.18). It should be noted that this is a departure from the magnetostatic case with the vector potential formulation in which n = 1 corresponds to a dipole arrangement in the axisymmetric case.
-there is a net charge with (i.e. the geometry looks like a point charge when viewed from a distance),
-whereas the n= 1 case corresponds to a dipole charge distribution in 2D planar problems. If charge
+Some care must be used in applying this boundary condition. Most of the time, it is sufficient to take n = 1 (i.e the objects in the solution region look like a dipole when viewed from a large distance). However, there are other cases (e.g. a 4-pole halbach permanent magnet array) in which the leading harmonic is something other than n= 1. You need to use your insight into your specific problem to pick the appropriate n for the leading harmonic. You also must put the objects of interest roughly in the center of the circular finite element domain to minimize the magnitude higher-order field components at the outer boundary.
-is conserved in the geometry of interest in the axisymmetric case, one needs to use n =2 in Eq.
-(A.18). It should be noted that this is a departure from the magnetostatic case with the vector
+Although the application of this boundary condition requires some thought on the part of the user, the results can be quite good. Figure A.7, corresponding to the axi1 example, represents the field produced by an air-cored coil in free space. The asymptotic boundary condition has been applied to the circular outer boundary. Inspecting the solution, flux lines appear to cross the circular boundary as if the solution domain were truly unbounded.
-potential formulation in which n = 1 corresponds to a dipole arrangement in the axisymmetric
-case.
+A quick note on computational efficiency: applying the absorbing boundary condition imposes no additional computing cost on the problem. The ABC is computationally no more time-consuming to apply than enforcing A = 0 at the outer boundary. Solution times for the PCG solver are equivalent in either case. It can also readily be derived that the ABC works exactly the same for harmonics problems. (To see this, just assume that the am in (A.12) can be complex valued, and follow the same derivation).
-Some care must be used in applying this boundary condition. Most of the time, it is sufficient
-to take n = 1 (i.e the objects in the solution region look like a dipole when viewed from a large
-distance). However, there are other cases (e.g. a 4-pole halbach permanent magnet array) in which
+=== A.3.3 Kelvin Transformation ===
-the leading harmonic is something other than n= 1. You need to use your insight into your specific
-problemto pick the appropriate n for the leading harmonic. You alsomust put the objects of interest
+A particularly good approach to “open boundary” problems is the Kelvin Transformation, a technique first discussed in the context of computational magnetics in [18] and [19]. The strengths of this technique are:
-roughly in the center of the circular finite element domain to minimize the magnitude higher-order
+  * the effects of the exterior region are, in theory, exactly modeled by this approach;
-field components at the outer boundary.
+  * a sparse matrix representation of the problem is retained (unlike FEM-BEM methods, which give the same “exact solution” but densely couples together the boundary nodes).
-Although the application of this boundary condition requires some thought on the part of the
+  * requires no “special” features in the finite element solver to implement the technique, other than the ability to apply periodic boundary conditions.
-user, the results can be quite good. Figure A.7, corresponding to the axi1 example, represents
-the field produced by an air-cored coil in free space. The asymptotic boundary condition has been
+The purposes of this note are to explain what the Kelvin transformation is derived and to show how it is implemented in the context of the FEMM finite element program.
-applied to the circular outer boundary. Inspecting the solution, flux lines appear to cross the circular
-boundary as if the solution domain were truly unbounded.
+** Derivation **
-A quick note on computational efficiency: applying the absorbing boundary condition imposes
+In the “far field” region, the material is typically homogeneous (e.g. air and free of sources. In this case, the differential equation that describes vector potential A is the Laplace equation:
-no additional computing cost on the problem. The ABC is computationally no more timeconsuming
-to apply than enforcing A = 0 at the outer boundary. Solution times for the PCG solver
+(A.19) $ \nabla^2 A = 0  $
-are equivalent in either case. It can also readily be derived that the ABC works exactly the same
-for harmonics problems. (To see this, just assume that the am in (A.12) can be complex valued,
+If we write (A.19) in polar notation, //A// is described by:
-and follow the same derivation).
-A.3.3 Kelvin Transformation
+(A.20) $ \frac{1}{r} \left( r \frac{\partial A}{\partial r} \right) + \frac{1}{r^2} \frac{\partial^2 A}{\partial θ^2} = 0  $
-A particularly good approach to “open boundary” problems is the Kelvin Transformation, a technique
-first discussed in the context of computational magnetics in [18] and [19]. The strengths of
+Assume that the “near field” region of the problem can be contained in a circle of radius //ro// centered at the origin. The far-field region is then everything outside the circle.
-this technique are:
-• the effects of the exterior region are, in theory, exactly modeled by this approach;
+One approach to unbounded problems is to attempt to map the unbounded region onto a bounded region, wherein problems can more easilby be solved. Specifically, we desire a way to transform the unbounded region outside the circle into a bounded region. One simple way to make such a mapping is to define another variable, //R//, that is related to //r// by:
-• a sparse matrix representation of the problem is retained (unlike FEM-BEM methods, which
-give the same “exact solution” but densely couples together the boundary nodes).
+(A.21) $ R = \frac{r_0^2}{r} $
-• requires no “special” features in the finite element solver to implement the technique, other
-than the ability to apply periodic boundary conditions.
+By inspecting (A.21), it can be seen that this relationship maps the exterior region onto a circle of radius //ro//.
-The purposes of this note are to explain what the Kelvin transformation is derived and to show how
-it is implemented in the context of the FEMM finite element program.
+The next step is to transform (A.19), the differential equation that the field must satisfy, into the mapped space. That is, (A.19) must be written in terms of //R// and //θ// rather than //r// and //θ//. We can evaluate derivatives in terms of //R// instead of //r// by employing the chain rule:
-Derivation
-In the “far field” region, the material is typically homogeneous (e.g. air and free of sources. In this
+(A.22) $ \frac{\partial}{\partial r} = \frac{\partial}{\partial R} \left( \frac{\partial R}{\partial r} \right) = - \frac{\partial}{\partial R} \left( \frac{R}{r_0} \right)^2  $
-case, the differential equation that describes vector potential A is the Laplace equation:
-Ñ2A = 0 (A.19)
+Now, we can note that at //r// = //R// = //ro//,
-If we write (A.19) in polar notation, A is described by:
+(A.23)  $ \frac{\partial A}{\partial r} = -\frac{\partial A}{\partial R} $
-r
-¶
-¶r
-r
-¶A
-¶r
-+
-r2
-¶2A
-¶q2 = 0 (A.20)
-Assume that the “near field” region of the problem can be contained in a circle of radius ro centered
-at the origin. The far-field region is then everything outside the circle.
-One approach to unbounded problems is to attempt to map the unbounded region onto a
-bounded region, wherein problems can more easilby be solved. Specifically, we desire a way
-to transform the unbounded region outside the circle into a bounded region. One simple way to
-make such a mapping is to define another variable, R, that is related to r by:
-R =
-r2
-o
-r
-(A.21)
-By inspecting (A.21), it can be seen that this relationship maps the exterior region onto a circle of
-radius ro.
-The next step is to transform (A.19), the differential equation that the field must satisfy, into
-the mapped space. That is, (A.19) must be written in terms of R and q rather than r and q. We can
-evaluate derivatives in terms of R instead of r by employing the chain rule:
-¶
-¶r
-=
-¶
-¶R
-dR
-dr
-= −
-¶
-¶R
-R
-ro
-(A.22)
-Now, we can note that at r = R = ro,
-¶A
-¶r
-= −
-¶A
-¶R
-(A.23)
 and we can substitute (A.22) into (A.19) to yield, after some algebraic manipulation:
-R
+(A.24) $  \frac{1}{R} \frac{\partial}{\partial R} \left(R \frac{\partial A}{\partial R} \right) + \frac{1}{R^2}  \frac{\partial^2 A}{\partial θ^2}  = 0 $
-¶
-¶R
+Eq. (A.24), the transformed equation for the outer region, has exactly the same form as inner region, only in terms of //R// rather than //r//. The implication is that for the 2-D planar problem, the exterior can be modeled simply by creating a problem domain consisting of two circular regions: on circular region containing the items of interest, and an additional circular region to represent the “far field.” Then, periodic boundary conditions must be applied to corresponding edges of the circle to enforce the continuity of //A// at the edges of the two regions. The is continuity of //A// at the boundary between the exterior and interior regions. For a finite element formulation consisting of first-order triangles, (A.23) is enforced automatically at the boundaries of the two regions. The second circular region exactly models the infinite space solution, but does it on a bounded domain - one could always back out the field for any point in space by applying the inverse of (A.21).
-R
+** Kelvin Transformation Example ** – ''open1.fem''
-¶A
-¶R
+As an example, consider an E-core lamination stack with a winding around it. Suppose that the objective is to determine the field around the E-core in the absence of any flux return path (i.e. when the magnetic circuit is open). In this case, the flux is not constrained to flow in a path that is a priori well defined, because the laminations that complete the flux path have been removed.
-+
+The geometry was chosen arbitrarily, the purpose here being more the procedure than the actual problem. The E-core was chosen to have a 0.5” thick center leg, 0.25” thick outer legs, and a slot depth of 0.75”. The material for the core is linear with a relative permeability of 2500. The coil carries a bulk current density of 2 MA/m2. The input geometry is picture in Figure A.8.
-R2
+In Figure A.8, the core is placed within a circular region, and a second circular region is drawn next to the region containing the core. Periodic boundary conditions are applied to the arcs that define the boundaries as shown in Figure A.8. The way that periodic boundary conditions are implemented in FEMM, each periodic boundary condition defined for the problem is to be applied to two and only two corresponding entities. In this case, each boundary circle is composed of two arcs, so two periodic boundary conditions must be defined to link together each arc with in the domain with the core to its corresponding arc in the domain representing the exterior region.
-¶2A
-¶q2 = 0 (A.24)
+Also notice that a point has been drawn in the center of the exterior region. A point property has been applied to this point that specifies that //A// = 0 at this reference point. The center of the circle maps to infinity in the analogous open problem, so it makes sense to define, in effect, //A// = 0 at infinity. If no reference point is defined, it is fairly easy to see that the solution is only unique to within a constant. The situation is analogous to a situation where Neumann boundary conditions have been defined on all boundaries, resulting in a non-unique solution for A. Due to the type of solver that FEMM employs, the problem can most likely be solved even if a reference point is not defined. However, defining a reference point eliminates the possibility of numerical difficulties due to uniqueness issues.
-Eq. (A.24), the transformed equation for the outer region, has exactly the same form as inner
-region, only in terms of R rather than r. The implication is that for the 2-D planar problem, the
+{{fig_a-8.png}}
-exterior can be modeled simply by creating a problem domain consisting of two circular regions:
-on circular region containing the items of interest, and an additional circular region to represent
-the “far field.” Then, periodic boundary conditions must be applied to corresponding edges of the
-circle to enforce the continuity of A at the edges of the two regions. The is continuity of A at the
-boundary between the exterior and interior regions. For a finite element formulation consisting
-of first-order triangles, (A.23) is enforced automatically at the boundaries of the two regions. The
-second circular region exactly models the infinite space solution, but does it on a bounded domain–
-one could always back out the field for any point in space by applying the inverse of (A.21).
-Kelvin Transformation Example – open1.fem
-As an example, consider an E-core lamination stack with a winding around it. Suppose that the
-objective is to determine the field around the E-core in the absence of any flux return path (i.e.
-when the magnetic circuit is open). In this case, the flux is not constrained to flow in a path that is
-a priori well defined, because the laminations that complete the flux path have been removed.
-The geometry was chosen arbitrarily, the purpose here being more the procedure than the actual
-problem. The E-core was chosen to have a 0.5” thick center leg, 0.25” thick outer legs, and a slot
-depth of 0.75”. The material for the core is linear with a relative permeability of 2500. The coil
-carries a bulk current density of 2 MA/m2. The input geometry is picture in Figure A.8.
-In Figure A.8, the core is placed within a circular region, and a second circular region is drawn
-next to the region containing the core. Periodic boundary conditions are applied to the arcs that
-define the boundaries as shown in Figure A.8. The way that periodic boundary conditions are
-implemented in FEMM, each periodic boundary condition defined for the problem is to be applied
-to two and only two corresponding entities. In this case, each boundary circle is composed of two
-arcs, so two periodic boundary conditions must be defined to link together each arc with in the
-domain with the core to its corresponding arc in the domain representing the exterior region.
-Also notice that a point has been drawn in the center of the exterior region. A point property
-has been applied to this point that specifies that A = 0 at this reference point. The center of the
-circle maps to infinity in the analogous open problem, so it makes sense to define, in effect, A = 0
 Figure A.8: Example input geometry.
-at infinity. If no reference point is defined, it is fairly easy to see that the solution is only unique to
-within a constant. The situation is analogous to a situation where Neumann boundary conditions
-have been defined on all boundaries, resulting in a non-unique solution for A. Due to the type of
+The resulting solution is shown in Figure A.9. As is the intention, the flux lines appear to cross out of the of the region containing the core as if unaffected by the presence of the boundary. The flux lines reappear in the domain representing the exterior region, completing their flux paths through the exterior region.
-solver that FEMM employs, the problem can most likely be solved even if a reference point is not
-defined. However, defining a reference point eliminates the possibility of numerical difficulties
-due to uniqueness issues.
+==== A.4 Nonlinear Time Harmonic Formulation ====
-The resulting solution is shown in Figure A.9. As is the intention, the flux lines appear to
-cross out of the of the region containing the core as if unaffected by the presence of the boundary.
+Starting with the the 3.3 version of FEMM, the program includes a “nonlinear time harmonic” solver. In general, the notion of a “nonlinear time harmonic” analysis is something of a kludge. To obtain a purely sinusoidal response when a system is driven with a sinusoidal input, the system must, by definition, be linear. The nonlinear time harmonic analysis seeks to include the effects of nonlinearities like saturation and hysteresis on the fundamental of the response, while ignoring higher harmonic content. This is a notion similar to “describing function analysis,” a widely used tool in the analysis of nonlinear control systems. There are several subtly different variations of the formulation that can yield slightly different results, so documentation of what has actually been implement is important to the correct interpretation of the results from this solver.
-The flux lines reappear in the domain representing the exterior region, completing their flux paths
-through the exterior region.
+An excellent description of this formulation is contained in [21]. FEMM formulates the nonlin-ear time harmonic problem as described in this paper. Similar to Jack and Mecrow, FEMM derives an apparent BH curve by taking H to be the sinusoidally varying quantity. The amplitude of B is obtained by taking the first coefficient in a Fourier series representation of the resulting B. For the purposes of this Fourier series computation, FEMM interpolates linearly between the user-defined points on the BH curve to get a set of points with the same H values as the input set, but with an adjusted B level. The rationale for choosing H to be the sinusoidal quantity (rather than B) is that choosing B to be sinusoidal shrinks the defined BH curve–the B values stay fixed while the H values become smaller. It then becomes hard to define a BH curve that does not get interpolated. In contrast, with H sinusoidal, the B points are typically larger than the DC flux density levels, creating a curve with an expanded range.
-A.4 Nonlinear Time Harmonic Formulation
-Starting with the the 3.3 version of FEMM, the program includes a “nonlinear time harmonic”
+{{fig_a-9.png}}
-solver. In general, the notion of a “nonlinear time harmonic” analysis is something of a kludge.
-To obtain a purely sinusoidal response when a system is driven with a sinusoidal input, the system
-must, by definition, be linear. The nonlinear time harmonic analysis seeks to include the effects
-of nonlinearities like saturation and hysteresis on the fundamental of the response, while ignoring
-higher harmonic content. This is a notion similar to “describing function analysis,” a widely used
-tool in the analysis of nonlinear control systems. There are several subtly different variations of
-the formulation that can yield slightly different results, so documentation of what has actually been
-implement is important to the correct interpretation of the results from this solver.
-An excellent description of this formulation is contained in [21]. FEMMformulates the nonlin-
 Figure A.9: Solved problem.
-ear time harmonic problem as described in this paper. Similar to Jack and Mecrow, FEMM derives
-an apparent BH curve by taking H to be the sinusoidally varying quantity. The amplitude of B is
+A “nonlinear hysteresis lag” parameter is then applied to the effective BH curve. The lag is assumed to be proportional to the permeability, which gives a hysteresis loss that is always proportional to |B|<sup>2</sup>. This form was suggested by O’Kelly [22]. It has been suggested that that the Steinmetz equation could be used to specify hysteresis lag, but the Steinmetz equation is badly behaved at low flux levels (i.e. one can’t solve for a hysteresis lag that produces the Steinmetz |B|<sup>1.6</sup> form for the loss as B goes to zero.)
-obtained by taking the first coefficient in a Fourier series representation of the resulting B. For the
-purposes of this Fourier series computation, FEMM interpolates linearly between the user-defined
+For nonlinear in-plane laminations, an additional step is taken to obtain an effective BH curve that also includes eddy current effects. At each H level on the user-defined BH curve, a 1D nonlinear time harmonic finite element problem is solved to obtain the total flux that flows in the lamination as a function of the H applied at the edge of the lamination. Then dividing by the lamination thickness and accounting for fill factor, and effective B that takes into account saturation, hysteresis, and eddy currents in the lamination is obtained for each H.
-points on the BH curve to get a set of points with the same H values as the input set, but with
-an adjusted B level. The rationale for choosing H to be the sinusoidal quantity (rather than B) is
-that choosing B to be sinusoidal shrinks the defined BH curve–the B values stay fixed while the H
-values become smaller. It then becomes hard to define a BH curve that does not get interpolated.
-In contrast, with H sinusoidal, the B points are typically larger than the DC flux density levels,
-creating a curve with an expanded range.
-A “nonlinear hysteresis lag” parameter is then applied to the effective BH curve. The lag
-is assumed to be proportional to the permeability, which gives a hysteresis loss that is always
-proportional to |B|2. This form was suggested by O’Kelly [22]. It has been suggested that that
-the Steinmetz equation could be used to specify hysteresis lag, but the Steinmetz equation is badly
-behaved at low flux levels (i.e. one can’t solve for a hysteresis lag that produces the Steinmetz
-|B|1.6 form for the loss as B goes to zero.)
-For nonlinear in-plane laminations, an additional step is taken to obtain an effective BH curve
-that also includes eddy current effects. At each H level on the user-defined BH curve, a 1D nonlinear
-time harmonic finite element problem is solved to obtain the total flux that flows in the
-lamination as a function of the H applied at the edge of the lamination. Then dividing by the lamination
-thickness and accounting for fill factor, and effective B that takes into account saturation,
-hysteresis, and eddy currents in the lamination is obtained for each H.

Unofficial FEMM documentation by Dr Stan Zurek

User Tools

Site Tools

Differences

Page Tools

Unofficial
FEMM
documentation
by Dr Stan Zurek