Nonlinear magnetostatic BEM formulation using one unknown double layer charge

Document Sample
Nonlinear magnetostatic BEM formulation using one unknown double layer charge Powered By Docstoc
					COMPEL: The International Journal for Computation and Mathematics in
Electrical and Electronic Engineering
Emerald Article: Nonlinear magnetostatic BEM formulation using one unknown
double layer charge
Kazuhisa Ishibashi, Zoran Andjelic

Article information:
To cite this document: Kazuhisa Ishibashi, Zoran Andjelic, (2011),"Nonlinear magnetostatic BEM formulation using one unknown
double layer charge", COMPEL: The International Journal for Computation and Mathematics in Electrical and Electronic Engineering,
Vol. 30 Iss: 6 pp. 1870 - 1884
Permanent link to this document:
Downloaded on: 15-05-2012
References: This document contains references to 15 other documents
To copy this document:
This document has been downloaded 165 times.

Access to this document was granted through an Emerald subscription provided by UNIVERSITY OF NEWCASTLE

For Authors:
If you would like to write for this, or any other Emerald publication, then please use our Emerald for Authors service.
Information about how to choose which publication to write for and submission guidelines are available for all. Additional help
for authors is available for Emerald subscribers. Please visit for more information.
About Emerald
With over forty years' experience, Emerald Group Publishing is a leading independent publisher of global research with impact in
business, society, public policy and education. In total, Emerald publishes over 275 journals and more than 130 book series, as
well as an extensive range of online products and services. Emerald is both COUNTER 3 and TRANSFER compliant. The organization is
a partner of the Committee on Publication Ethics (COPE) and also works with Portico and the LOCKSS initiative for digital archive
                                                                        *Related content and download information correct at time of download.
                                                   The current issue and full text archive of this journal is available at

30,6                                       Nonlinear magnetostatic BEM
                                          formulation using one unknown
                                                double layer charge
                                                             Kazuhisa Ishibashi and Zoran Andjelic
                                                 Corporate Research, ABB Switzerland Ltd, Baden, Switzerland

                                        Purpose – The purpose of this paper is to solve generic magnetostatic problems by BEM, by
                                        studying how to use a boundary integral equation (BIE) with the double layer charge as unknown
                                        derived from the scalar potential.
                                        Design/methodology/approach – Since the double layer charge produces only the potential gap
                                        without disturbing the normal magnetic flux density, the field is accurately formulated even by one
                                        BIE with one unknown. Once the double layer charge is determined, Biot-Savart’s law gives easily the
                                        magnetic flux density.
                                        Findings – The BIE using double layer charge is capable of treating robustly geometrical
                                        singularities at edges and corners. It is also capable of solving the problems with extremely high
                                        magnetic permeability.
                                        Originality/value – The proposed BIE contains only the double layer charge while the conventional
                                        equations derived from the scalar potential contain the single and double layer charges as unknowns.
                                        In the multiply connected problems, the excitation potential in the material is derived from the
                                        magnetomotive force to represent the circulating fields due to multiply connected exciting currents.
                                        Keywords Boundary integral equation, Double layer charge, Multiply connected problem,
                                        Nonlinear magnetostatic analysis, Scalar potential, Integral equations, Electric current
                                        Paper type Research paper

                                        1. Introduction
                                        The magnetostatic analysis by the boundary integral equation (BIE) derived from the
                                        scalar potential given in the book by Stratton (1941) is one of the basic numerical
                                        approaches and many papers have been reported (Rucker and Richter, 1988; Koizumi
                                        et al., 1990; Sawa and Hirano, 1990; Rucker et al., 1992, Minciunescu, 1998; Buchau et al.,
                                        2003, 2007). They utilize the scalar potential wH for the magnetic field H, and derive the
                                        BIEs with two unknowns of single and double layer charges, ss and sd as the state
                                                                                                           ^       ^
                                        variables. Hereafter, we call the BIEs “conventional BIE”. Tozoni and Mayergoyz
                                        (1974) have derived another BIE with the double layer charge sd from an integral
                                        representation of the scalar potential wB for the magnetic flux density B. It seems there
                                        are several advantages described below over the conventional BIEs, but only a few
                                        papers have been reported.
                                            In this paper, we derive a BIE with the help of the concept of wB and adopting the
                                        volume integral equation (VIE) approach (Wang, 1991). The double layer charge sd is
COMPEL: The International Journal
for Computation and Mathematics in      defined as two gapless layer charges with opposite poles of the same values, and so sd
Electrical and Electronic Engineering   does not produce any gap between the normal components of B but it produces a
Vol. 30 No. 6, 2011
pp. 1870-1884                           potential gap between the layers. Since the continuity of B is automatically fulfilled by
q Emerald Group Publishing Limited      wB, it formulates easily the magnetostatic fields. With the help of the potential gap,
DOI 10.1108/03321641111168165           we derive a BIE with one unknown sd by enforcing the boundary condition of the
continuity of the tangential component of H, while the other boundary condition of the            Nonlinear
continuity of the normal component of B is automatically fulfilled because sd does not         magnetostatic
affect it. By virtue of these properties, it is expected to solve accurately wB, which is
definite and smoothly distributed even if B becomes infinite. Since sd given by solving       BEM formulation
the BIE is equivalent to the loop current, it is easy to evaluate B with the help of
Biot-Savart’s law even if B approaches infinity at the edge and corner.
   In the multiply connected problem (MCP), we introduce the excitation potential in                 1871
the material derived from the magnetomotive force (MMF) so as to meet easily the
requirement that the line integral of H along the toroidal core must equal the total
MMF due to the exciting currents, and then we formulate the magnetostatic field by
adopting the surface integral equation approach (Wang, 1991) with the help of the
scalar potential wH for H (Stratton, 1941). While the original BIEs derived from wH
contain ss and sd as unknowns, utilizing the boundary conditions, we derive a BIE
           ^      ^
only with sd .
   This paper also includes how to solve nonlinear magnetostatic problems.

2. Integral representation of scalar potential
The magnetostatic fields are represented as:
                                      7£H ¼ J;                                        ð1Þ
                                       7·B ¼ 0                                        ð2Þ
where H, B and J are the magnetic field, flux density and the electric current density. In
the VIE approach given in the book by Wang (1991), a material with the surface S and
volume V is replaced to the magnetization M as follows.
   With the magnetic permeability of free space m0, we rewrite equation (2) as:
                             m0 7 · H ¼ 27 · ðB 2 m0 H Þ:                             ð3Þ
Equation (3) is for determining the magnetic field H in the free space with the magnetic
charge distribution. The right side is regarded as fictitious magnetic charges due to M,
which is defined as:
                                   M ¼ B 2 m0 H :                                     ð4Þ
Since equation (1) is not affected by the material property, it remains as it is.
   We rewrite the integral representation of the potential wM due to M given in the
book by Bleaney and Bleaney (1976) as:
                                       M ·r
                       wMP ¼ mrio            dV
                                       4pr 3
                                   VZ            Z             
                                        M · dS        7·M
                            ¼ mrio             2           dV                     ð5Þ
                                         4pr        V 4pr
                                   ZS            Z           
                                        MS            mv
                            ¼ mrio           dS þ         dV
                                      S 4pr         V 4pr

where the subscript P denotes the value at an observation point Po, the relative
permeability mrio is defined as mrio ¼ m/m0 inside V, mrio ¼ 1 outside V, r is the
COMPEL   distance from the integral point Pi to Po, and mv and MS are called the volume and
         surface charges defined as:
                                   mv ¼ 27 · M ¼ 27 · ðB 2 m0 H Þ
                                  M S ¼ n · M ¼ n · ðB 2 m0 H Þ:
         with the unit outward normal n. The potential wM in equation (5) is composed of two
1872     potentials; the potential wMs due to MS and that wmv due to mv. The concept of magnetic
         shell (Stratton, 1941), which is composed of the double layer charges sd, suggests that
         wMs could be replaced by the potential due to sd. Taking this concept into account, we
         get the total potential wB as:
                                 wBP ¼ wBeP þ mrioP wmvP þ sd           dS                   ð7Þ
                                                              S   4pr 3
         where wBe and wmv are the potentials at Po produced, respectively, by the exciting
         source and mv, r is the distance from an integral point Pi on S to Po and n is the unit
         normal at Pi. The potential wmv is represented as:
                                            wmvP ¼         dV :                              ð8Þ
                                                     V 4pr

         3. Boundary integral equation
         The potentials on the material surfaces with the subscripts o and i denoting the outer
         and inner sides are given as:
                                                  VsP               n·r
                           wBoP ¼ wBeP þ wmvP þ        sdP þ sd           dS;               ð9Þ
                                                   4p         S     4pr 3
                                                       VsP                n·r
                      wBiP ¼ wBeP þ mrP wmvP 2 1 2          sdP þ sd            dS         ð10Þ
                                                       4p            S    4pr 3
         where Vs is the solid angle subtended at the singular point on the surface S (Bladel,
         1991). Applying the boundary condition of the continuous condition of the tangential
         magnetic field to equations (9) and (10), we derive a BIE. The boundary condition is
         given as:
                                          n P £ B osP n P £ B isP
                                                     ¼                                          ð11Þ
                                              m0          m
         where B isP ¼ mH isP with the subscript s denoting the value on S. Since B is given as
         B ¼ 2 grad(wB), equation (11) is rewritten as:
                                       wBos1 2 wBos2 wBis1 2 wBis2
                                                    ¼                                           ð12Þ
                                           m0 Dl          mDl
         where the subscripts 1 and 2 denote the values at the points P1 and P2 on the surface,
         and Dl is the distance from P1 to P2. Since equation (12) is identical to the condition that
         wBos/m0 ¼ wBis/m, we get the following BIE as:
                                    C sP sdP þ sd          dS ¼ 2wBeP                            ð13Þ
                                                 S   4pr 3
where CP is the nonlinear coefficient of sdP defined as:                                            Nonlinear
                                      VsP mrs 2 VsP þ 4p                                      magnetostatic
                               CP ¼                                                 ð14Þ
                                          4pðmrs 2 1Þ                                       BEM formulation
with the relative permeability mrs defined as mrs ¼ jBisj/jm0Hisj on the material
surface. The other boundary condition of the continuity of normal magnetic flux
density is automatically fulfilled because sd does not affect the original normal                     1873
component of B.
   The exciting potential wBe due to the exciting coil current Ic is given as follows.
Since the potential wBe is defined as the line integral of the magnetic flux density Be
produced by Ic from the original point P0 of the potential to Po, it is given as:
                                         Z Po
                                 wBeP ¼       B e · dL                            ð15Þ

where Be is given with the help of Biot-Savart’s law as:
                                           I c £ rc
                             B e ¼ m0 Lc            dLc                             ð16Þ
                                             4pr 3

with the distance rc from an integral point along the coil length Lc to Po. The exciting
potential wBe in equation (15) is evaluated numerically, or with the help of the integral
formula as:
                                wBeP ¼ m0 I c ðVcP þ dc Þ                           ð17Þ
                                             nc · r c
                                 VcP ¼             3
                                                      dS c                         ð18Þ
                                          S c 4pr c
with a flat surface Sc surrounded by a loop line current Ic and the unit normal nc on Sc,
and the value of dc is the adjusting term given as: dc ¼ 2 0.5 when VcP . 0 and
dc ¼ 0.5 when VcP , 0. In the magnetostatic analysis of simply connected problems,
the exciting potential must be continuous and dc works for wBe to be continuous while
dc is unnecessary in MCP because the discontinuity of Vc works as the value of the cut
surface. When the exciting field is produced by the a permanent magnet, the exciting
potential wBe is given as:
                                             M e · re
                                wBeP ¼             3
                                                       dV e                        ð19Þ
                                         V e 4pr e
where Me is the magnetization of the permanent magnet, Ve is its volume and re is the
distance from the point at a segmental volume dVe to Po.

4. Flux density due to double layer charge
Here, we evaluate the magnetic flux density Bs due to the double layer charge at
any point utilizing the relation between the double layer charge and loop current
(Stratton, 1941).
   The concept of magnetic shell is shown in Figure 1. The magnetic shell is composed
of the double layer charges sd, which are replaced by the loop currents Jl when each
COMPEL                          segmental charge is regarded as constant. This assumption leads to Bs given with the
30,6                            help of Biot-Savart’s law as follows:

                                                                      Ns          I
                                                                                         uJ £ r L
                                                             B sP ¼         sdi               3
                                                                                                  dLi                ð20Þ
                                                                      i¼1             DLi 4pr L

                                where Ns is the total number of sd that equals to the number of Jl, uJ is the direction of
                                Jl, which circulates anticlockwise along the periphery of sd, DL is the length of Jl, and
                                rL is the distance from the integral point along Jl to Po.
                                    When we obtain sd by adopting the constant element in discretization of equation
                                (13), we can evaluate easily B at any point by employing equation (20) because only the
                                periphery of the constant element is taken into account as DL in equation (20). In this
                                case, Ns equals the number of the total elements. The value on the surface is not given
                                by equation (20) because it is singular on the surface. How to evaluate the B on the
                                surface is described below. Since the constant element is incapable of setting the
                                unknown sd at the edge and corner, B at these points cannot be evaluated.
                                    In order to obtain B at the edge and corner, we adopt the linear element in
                                discretization of equation (13). The double layer charge sd at each node of the
                                surface elements is obtained by solving equation (13), and sde with the subscript e
                                denoting the value within the surface element is interpolated with sd at the nodes.
                                Figure 2 shows one of the surface elements with sdi (i ¼ 1, 2, 3 and 4) shown by the
                                symbol W at each node. We interpolate sde at the calculating point Pe within the surface
                                element as:

Figure 1.
Magnetic flux density due
to magnetic shell, which is
equivalent to either sd or Jl

                                                             σd4                                   σd3


Figure 2.
Interpolation of double
layer charge within                                          σd1                                   σd2
surface element
                                                                       a                   a
                                    sde ¼         N i sdi                            ð21Þ
                                                                                             BEM formulation
where the shape function Ni is given as:
                                N1 ¼        ða 2 xÞðb 2 yÞ;
                                      4ab                                                                      1875
                                N2 ¼        ða þ xÞðb 2 yÞ;
                                N3 ¼        ða þ xÞðb þ yÞ;
                                N4 ¼        ða 2 xÞðb þ yÞ
with the local coordinates, x and y, at Pe.
   Adopting sde interpolated as equation (21), we evaluate Bs as follows. The surface
element is divided further into the sub-elements as shown in Figure 3(a) for evaluating
Bs on the flat surface and Figure 3(b) for Bs at the edge and corner. The solid lines in
these fugures are for the original surface element and the broken lines are for
subdividing the surface element. The calculating points Po for obtaining sd1 2 sd4
have been set at the intersecting point of the solid lines shown by W. We set the point Pe
shown by V at the center of the sub-element surrounded by the solid and dotted lines
and evaluate sde by employing equation (21). The loop current is introduced so as to
circulate anticlockwise along the periphery of the sub-element as shown by the arrows
in the figure on the assumption that sde on the sub-elements are constant. When we
evaluate Bs at the vertex, we set one loop current Jl on the surface elements adjoining
the vertex as shown in Figure 3(b), where Jl, which is approximately equivalent to sd at
the vertex, circulates along each outer periphery of three adjoining sub-elements, and
others as shown in Figure 3(a).
   The total magnetic flux density B is given as:
                                B P ¼ B sP þ B eP þ B mvP                            ð22Þ
where Bs is the magnetic flux density due to Jl given as equation (20), Be due to the
exciting sources given in equation (16), and Bmv due to the equivalent magnetic charge
mv given in equation (6).
   We can apply equation (22) to evaluate B at any point except on the surface.
Even on the surface, equation (22) is capable of evaluating the normal component of B

               σd4                 σd3

                                                                                                              Figure 3.
                                                                                                    (a) Subdivisions for
                                                                                              evaluating magnetic flux
                                                                                             density; (b) subdivisions of
               σd1                 σd2                                                       surface elements adjoining
                          (a)                               (b)
COMPEL   by setting the calculating point Pe at the center of the sub-element, but incapable of
         doing the tangential component Bt. Since there is the gap DBs between the inner and
30,6     outer Bt as given in equation (11), we take DBs into account to evaluate Bt and get the
         following equations:
                                                            ðn P £ DB sP Þ
                                n P £ B oP ¼ n P £ B sP þ                  ;
1876                                                              2
                                                            ðn P £ DB sP Þ
                                n P £ B iP ¼ n P £ B sP   2
         where Bsp has been evaluated by employing equation (22). Substituting these
         equations into equation (11), DBs is determined as:
                                                    2ðm0 2 mÞ
                                   n p £ DB sP ¼              n p £ B sP :
                                                     m þ m0
         Finally, we get the tangential magnetic flux density as:
                                     n P £ B oP ¼          n P £ B sP ;
                                                    m þ m0
                                     n P £ B iP   ¼        n P £ B sP :
                                                    m þ m0
         Since the potential wB on the surface has already been given by solving equation (13),
         Bt also is evaluated as:
                                                 wB1 2 wB2
                                            Bti ¼          ;
                                                 m0 Bti                                     ð24Þ
                                           Bto ¼
         where wB1 and wB2 are the potentials at P1 and P2 on the surface, and Dl is the distance
         between these points. The magnetic field is evaluated as HP ¼ BP/m.

         5. Multiply connected problems
         Here, we introduce an integral representation for the exciting potential in the material
         so as to meet easily the circulating magnetic field produced by the exciting coil current
         Ic in the MCP. We derive the exciting potential wHem for H from the MMF produced in
         the material by Ic. A typical MCP is shown in Figure 4 where the segmental MMF dIm
         is defined as:
                                     dI m ¼ H m · dLm dS m ¼ H m dV                         ð25Þ
         where Hm is the magnetic field due to Ic, dSm is the surface of the segmental MMF, dLm
         is the thickness of the MMF with the direction of the same as that of Hm and dV is the
         segmental volume of the toroidal core. The potential dwHem at an observation point Po
         due to dIm is given as:
                                                            nm · r
                                          d wHemP ¼ dI m                                    ð26Þ
                                                            4pr 3
                             Po         r                    Hm                                    magnetostatic
                                                                                                 BEM formulation
                                    dϕmP                      dSm

                                                           dIm = Hm• dLmdSm
                                  Exciting                          Toroidal
                                  coil                              core
                                                                                                                Figure 4.
                                                                                                 Segmental MMF dIm due
                                                                                                   to exciting coil current

where, r is the distance from dIm to Po and nm is the unit vector perpendicular to dSm.
Taking the direction of nm is the same as that of Hm into consideration, the potential
wHem at Po due to the total dIm in the toroidal core is given as:
                                            r ·H m
                               wHemP ¼           3
                                                    dV :                            ð27Þ
                                           V 4pr

The magnetic field Hm due to Ic is evaluated with the help of Biot-Savart’s law.
   The concept of wHem based on the segmental MMF is similar to that based on the cut
surface. The difference is as follows. The potential wHec due to the coil current is given
as the right side of equation (17) divided by m0 on condition that dc ¼ 0, that is,
wHecP ¼ I c VcP , which is regarded as due to the total MMF concentrated on the cut
surface while equation (27) is due to the segmental MMF distributed in the entire core.
   Next, adopting wHem in equation (27) as the exciting potential in the material, we
shall derive BIEs for the MCP with the help of the SIE approach (Wang, 1991) by using
the potential wH for H (Stratton, 1941). The potentials in the outer and inner sides of
material, wHo and wHi, are given as:
                                   Z                  Z
                                       ›wHo 1                  n·r
                  wHoP ¼ wHecP þ                dS 2 wHo             dS;              ð28Þ
                                     S ›n 4pr           S     4pr 3
                                    Z                  Z
                                        ›wHi 1                 n·r
                   wHiP ¼ wHemP þ                dS 2 wHi            dS               ð29Þ
                                      S ›n 4pr           S    4pr 3
where m0(›wHoP/›n) ¼ m(›wHiP/›n) for the normal B and wHoP ¼ wHiP for the tangential
H on the surface. We call wH “double layer charge sd ” and call ›wHiP/›n “single layer
charge ss ”. Utilizing m0(›wHoP/›n) ¼ m(›wHiP/›n) on the surface, which implies that:
                           Z                   Z
                             ›wHo 1               ›wHi 1
                        m0             dS ¼ m                dS;
                            S ›n 4pr            S ›n 4pr

we get BIE for solving sd as:
                      ^                                   n·r
             ðmr þ 1Þ     þ ðmr 2 1Þ                sd
                                                    ^           dS ¼ mr wHemP 2 wHecP :   ð30Þ
                       2                        S         4pr 3
COMPEL   Once we have determined sd by solving equation (30), we obtain ss from BIE derived by
                                     ^                                    ^
         equating wHoP in equation (28) and wHiP in equation (29) on the surface as:
30,6                        Z                                      Z
                           1           1                                    n·r
                      1þ           ss
                                    ^     dS ¼ 2wHemP 2 wHecP þ 2 sd     ^        dS:      ð31Þ
                           mr S 4pr                                    S    4pr 3
         It is important to notice that we must adopt the constant element in discritization of
1878     equations (30) and (31) because each potential adjacent to MMF is intrinsically
         discontinuous even if the segmental MMF is extremely thin.

         6. Iterative solution of nonlinear BIE
         If the magnetization property is linear, the unknowns of the BIE are sd of the surface
         element, but if it is nonlinear, the unknowns are sd and mv of the volume element.
         Since the number of volume elements are much bigger than that of surface elements,
         it requires large amount of computer memory to solve sd and mv. In order to solve with
         small amount of computer memory, we propose an iterative approach adopting
         Newton-Raphson method combined with a simple iterative approach (called “combined
         iterative approach”).

         A. Magnetic flux density due to volume charge
         At first, we examine how to evaluate mv in equation (6) and Bmv due to mv in equation
         (22). When we adopt the constant volume element in discritization employed for the
         magnetic moment method (Takahashi et al., 2007), the surface magnetic charges mvs
         appear on the surface of the volume element. We shall divide the inside of the magnetic
         material into Nv volume elements and approximate mv as:
                                        mvs ¼ 1 2       n vs · B i                          ð32Þ
         where nvs is the unit normal on the element surface, and mr is the relative permeability
         defined as mr ¼ jBij/jm0Hij. The potential wmvs in equation (8) due to mvs is given as:
                                                         1 n vs · B i
                                 wmvsP ¼             12               dS                    ð33Þ
                                          Nv    DS v    mr     4pr

         where Nv is the number of the volume elements, and DSv is the area of the surface on
         the volume element.
            Since mvs is not equivalent to mv itself but to mv including sd, BmvP in equation (22)
         is given as:

             B mvP ¼ B mvsP 2 B MsP
                    ¼ ðmrP 2 1ÞB eP
                                        1 ðn vs · B i Þr      XZ       uJ £ r
                      þ mrP          12            3
                                                         dS 2       sd        dL
                            Nv  DS v    mr   4pr              Ns DL    4pr 3

         where Bmvs is the magnetic flux density due to mvs, BMs is that due to the surface
         magnetic charge Ms defined in equation (6), Ns is the number of the surface elements
         and DL is the periphery of the surface element.
B. Iterative solution of nonlinear equation                                                      Nonlinear
In Newton-Raphson method, the initial value of sd obtained as the solution of equation       magnetostatic
(13) is improved by the correction dsd, which is given by solving the following BIE as:
                                                                                           BEM formulation
                       ›C P sdP               n·r
                                dsdP þ dsd          dS ¼ 2dwBP                     ð35Þ
                         ›sdP            S    4pr 3
where dwBP is the residual defined as the difference between the right and left sides of
equation (13) and given as:
                      dwBP ¼ wBeP þ C P sdP þ sd            dS:                    ð36Þ
                                                 S    4pr 3
The coefficients of the term of Jacobian matrix in equation (35) is given as:

                        ›C P sdP             ›C P ›mrs ›BsP
                                 ¼ C P þ sdP                                       ð37Þ
                         ›sdP                ›mrs ›BsP ›sdP
with the absolute value Bs of Bs because Cp is not the explicit function of sd but
of mrs, which is the function of Bs. The solution sold obtained previously is improved
by dsd as:

                                  snew ¼ sold þ dsd :
                                   d      d                                        ð38Þ

In the simple iterative approach, we evaluate mvs in equation (32), which works for
improving sd. The evaluation of mvs requires Bi obtained according to equation (22)
with the help of sd newly obtained and also mvs previously obtained. In this way, mvs is
successively improved at the end of each iteration in the combined iterative approach.
   The nonlinear BIE (equation (13)) is solved iteratively with the nonlinear
permeability m, which are approximated according to the B-H curve as:

             m ¼ m0 ð1; 400 2 1; 900ðBi 2 0:8Þ2 Þ; when 0 # Bi , 1:5T
               m ¼ m0 ð1 þ 468 expð26ðBi 2 1:5ÞÞÞ; when 1:5T # Bi :

The flowchart to solve the nonlinear magnetostatic problems is shown in Figure 5:
  (1) We evaluate the source potential wBe.
  (2) We set the initial value of the surface relative permeability mrs.
  (3) We solve equation (13) to obtain the initial values of sd.
  (4) Employing equation (22), we evaluate Bi and obtain mr.
  (5) With the help of equation (23), we evaluate the magnetic flux density Bs on the
      surface in order to obtain mrs.
  (6) We evaluate the residual dwBP in equation (36).
  (7) We solve equation (35) to obtain the correction dsd.
  (8) With dsd, we improve sd in equation (38).
  (9) We get back to (4) and repeat the procedures from (4) to (9) until the value sd
COMPEL                                                                   START
                                                                      Evaluate jBe.

                                                                         Set mrs.

1880                                                             Solve (13) to obtain sd.

                                                            Evaluate Bi in (22) to obtain mr.
                                                        Evaluate Bs in (22) and (23) to obtain mrs.

                                                            Evaluate Jacobian Matrix in (37).
                                                                 Evaluate djB in (36).

                                                                 Solve (35) to obtain dsd

                                                                   Improve sd in (38).
Figure 5.
Flowchart for nonlinear
magnetostatic analysis                                             CONVERGENCE
using combined iterative

                           The convergence in Newton-Raphson approach is good, however in the simple iterative
                           routine it is difficult to get good convergence and adopting the deceleration coefficient
                           h we improve the evaluated Bi as follows:

                                                          Bi ¼ hBOld þ ð1 2 hÞBNew
                                                                 i             i                                  ð40Þ
                           where the superscripts Old and New denote the value Bi previously and newly
                           evaluated by equation (22). When the value of h is set larger, the convergence becomes
                           more certain but the number of iterations increases.

                           7. Numerical validation of proposed approach
                           A. Cylindrical core
                           In order to check the adequacy and effectiveness of the proposed approach, first we
                           solve an axisymmetric problem shown in Figure 6 (Magele et al., 1988). A cylindrical
                           core with radius 5 cm and height 20 cm is excited by a coil with inner and outer radii
                           6 and 8 cm, and height 15 cm. The cylindrical axis is parallel to the Z-axis. It is confirmed
                           that the computed results of B are almost equal to those given in the paper by Magele
                           et al. (1988). The coil current is increased to 40 kAT so as to evaluate B in more saturated
                           case. We show the computed results of the magnetic flux density B in Figure 7(a) and (b)
                           comparing with those by using the magnetic moment method (Takahashi et al., 2007).
                           The magnetic flux density Bz2 z along Z-axis is shown in Figure 7(a), and Bz2 r along the
                           R-axis shown in Figure 7(b), where the symbol W denotes the results by the proposed
                           method and x denotes those by the magnetic moment method.
                                                                                                                                      2 cm                                            Nonlinear
                                                                                        6 cm
                                                                                     5 cm                                    1 cm                                               BEM formulation

                                                                                                                                                  R                                              1881
                                                                         20 cm
                                                                                                                                          15 cm

                                                                Coil                                                            Coil                                                        Figure 6.
                                                            I = 40 kAT                                                      I = 40 kAT                                              Axisymmetric model

                                   2                                                                                        2
 Magnetic flux density Bz-z (T)

                                                                                          Magnetic flux density Bz-r (T)

                                                            Cylinder Air
                                  1.5                                                                                      1.5

                                                                                                                                       Cylinder         Coil Air
                                   1                                                                                        1

                                                                                                                                                  Air                                           Figure 7.
                                  0.5                                                                                      0.5
                                                                                                                                                                                  (a) Computed results of
                                                                                                                                                                                    magnetic flux density
                                   0                                                                                        0                                                       Bz2 z along Z-axis; (b)
                                        0               5            10            15                                            0               5            10           15         computed results of
                                                                                                                                                                                    magnetic flux density
                                            Distance along Z-axis from center [cm]                                                   Distance along R-axis from center (cm)
                                                                                                                                                                                        Bz2 r along R-axis
                                                             (a)                                                                                      (b)

B. Hollow sphere (shielding problem)
Next, we shall solve a problem with a magnetic hollow sphere in the uniform magnetic
field He (2,400 A/cm) shown in Figure 8, where the direction of He is parallel to the
Z-axis. We analyze linear and nonlinear cases. In the linear case, mr is supposed as
mr ¼ 1,000 and computed results are compared with analytical values given in the
book by Bleaney and Bleaney (1976). The computed results of B along X- and Z-axis,
together with analytical ones, are shown in Figure 9(a) and (b), where the symbols W
and x denote those by the proposed approach and FEM, and the solid line the



                                                                                        3 cm                                 5 cm                                                             Figure 8.
                                                             He = 2,400A/cm                                                                                                     Hollow sphere in uniform
                                                                                                                                                                                          magnetic field
COMPEL                       analytical ones. In the linear analysis, the results are almost the same as analytical
30,6                         values and omitted.

                             C. Toroidal core (MCP)
                             The last analysis model is shown in Figure 10. In this problem, the magnetization
                             property of toroidal core is supposed to be linear with the relative magnetic
1882                         permeability 500. The dimensions of the model is as follows: core radius Rt ¼ 10 cm,
                             core cross-section dimensions Lt ¼ 5 cm and Ht ¼ 5 cm, coil dimensions Lc ¼ 10 cm
                             and Hc ¼ 10 cm, coil angle uc ¼ 608. The coil current is 500 AT. The computed results
                             of the magnetic field H in the core at the midpoint of cross-section circulating along the
                             broken line from P are shown in Figure 11(a) and those along from Ps to Pe shown in
                             Figure 11(b). The symbols W and x show H computed by the proposed BIE and the
                             BIEs of the surface charge-current formulation (Ishibashi et al., 2011), respectively.
                             Both the results are almost the same inside the core but somewhat different outside the
                             core. It is too sensitive to obtain ss from sd by employing equation (31). And so, once
                                                                 ^       ^
                             sd and ss have been determined, ss is smoothed according to the following equation:
                              ^        ^                           ^

                                                                                                            Ne          Z
                                                                                                                                           nP · r
                                                                                                  ssP ¼ 2
                                                                                                                  ^                             3
                                                                                                            i¼1                         DS 4pr

                             where ss is the smoothed value of ss , Ne is the number of the total constant elements.

                                                             100                                                                                       100
                              Magnetic flux density Bz (T)

                                                                                                                        Magnetic flux density Bz (T)

                                                                    o, x : Nonlinear B-H                                                                      o, x : Nonlinear B-H
                                                             10–1        : Linear mr = 1,000                                                           10–1        : Linear mr = 1,000

                                                                                               Shielding                                                                                 Shielding

Figure 9.                                                    10–2                                                                                      10–2
(a) Computed results of
magnetic flux density B
along X-axis; (b) computed                                   10–3                                                                                      10–3
                                                                0      1     2      3      4     5      6                                                 0      1     2      3      4     5      6
results of magnetic flux
                                                              Distance from sphere center along X-axis (cm)                                             Distance from sphere center along Z-axis (cm)
density B along Z-axis
                                                                                   (a)                                                                                       (b)

                                                                                                                                                       Toroidal core
                                                                      Exciting coil
                                                                       I = 500AT
                                                                                                       Core radius Rt                                  P                         Lt
                                                                                                                                                                          Ps             Pe
                                                                                                  Coil angle qc                                                                               Hc
Figure 10.                                                                                                                                                                       Ht
Analysis model with
magnetic core excited by
coil currents                                                                                                                                                               Cross-section of
                                                                                                                                                                             core and coil
                                                   10                                                        Nonlinear

                         Magnetic field H (A/cm)
                                                                                                       BEM formulation

                                                        0   10     20    30     40      50    60
                                                            Distance along specified line (cm)

                         Magnetic field H (A/cm)

                                                                  Toroidal core
                                                   20                                                                 Figure 11.
                                                   15                                                     (a) Computed results of
                                                   10                                                     magnetic field H along
                                                    5                                                          circulating line; (b)
                                                    0                                                         computed results of
                                                     –5 –4 –3 –2 –1 0 1 2 3                  4     5      magnetic field H along
                                                           Position from Ps to Pe (cm)                 specified line from Ps to Pe
                                                                                                              shown in Figure 10

8. Conclusions
We have proposed a BIE with one unknown of the double layer charge sd. The main
difference between the conventional and proposed BIEs is as follows. The conventional
BIE is derived from scalar potential wH for the magnetic field H while the proposed BIE
from wB for the magnetic flux density B. This slight difference makes significant great
works. Even if the unknown of the BIE is one while the conventional BIEs contain two,
the boundary conditions are fulfilled completely and besides B is easily and directly
evaluated by Biot-Savart’s law because sd is equivalent to the loop current. Since sd at
the edge/corner is easily determined by the BIE and sd gives B, we can treat robustly
the geometrical singularities at the edge/corner. In order to solve the MCP, we have
derived a generic integral representation to give the exciting potential produced by the
MMF so as to meet arbitrarily shaped exciting coils. This makes the magnetostatic
analysis generic from the CAD-modeling point of view.

Bladel, J.V. (1991), Singular Electromagnetic Fields and Sources, Oxford University Press, Oxford.
Bleaney, B.I. and Bleaney, B. (1976), Electricity and Magnetism, 3rd ed., Oxford University Press,
COMPEL   Buchau, A., Hafla, W. and Rucker, W.M. (2007), “Accuracy investigations of boundary element
                methods for the solution of Laplace equations”, IEEE Trans. Magn., Vol. 43 No. 4, pp. 1225-8.
30,6     Buchau, A., Rucker, W.M., Rain, O., Rishmuller, V., Kurz, S. and Rjasanow, S. (2003),
                “Comparison between different approaches for fast and efficient 3-D BEM computations”,
                IEEE Trans. Magn., Vol. 39 No. 3, pp. 1107-10.
         Ishibashi, K., Andjelic, Z. and Pusch, D. (2011), “Magnetostatic analysis of multiple-connected
1884            problems by boundary integral equations”, Material Science Forum, Vol. 670, pp. 299-310,
                available at: (accessed 2011).
         Koizumi, M., Onisawa, M. and Utamura, M. (1990), “Three-dimensional magnetic field analysis
                using scalar potential formulated by boundary element method”, IEEE Trans. Magn.,
                Vol. 26 No. 2, pp. 360-3.
         Magele, Ch., Stoengner, H. and Preis, K. (1988), “Comparison of different finite element
                formulation for 3D magnetostatic problems”, IEEE. Trans. Magn., Vol. 24 No. 1, pp. 31-4.
         Minciunescu, P. (1998), “Contributions to integral equation method for 3D magnetostatic
                problems”, IEEE Trans. Magn., Vol. 34 No. 5, pp. 2461-4.
         Rucker, W.M. and Richter, K.R. (1988), “Three-dimensional magnetostatic field calculation using
                boundary element method”, IEEE Trans. Magn., Vol. 24 No. 1, pp. 23-6.
         Rucker, W.M., Magele, C., Schlemmer, E. and Richter, K.R. (1992), “Boundary element analysis of
                3-D magnetostatic problems using scalar potentials”, IEEE Trans. Magn., Vol. 28 No. 2,
                pp. 1099-102.
         Sawa, K. and Hirano, T. (1990), “An evaluation of the computational error near the boundary
                with magnetostatic field calculation by BEM”, IEEE Trans. Magn., Vol. 26 No. 2, pp. 403-6.
         Stratton, J.A. (1941), Electromagnetic Theory, McGraw-Hill, New York, NY.
         Takahashi, Y., Matsumoto, C. and Wakao, S. (2007), “Large-scale and fast nonlinear magnetostatic
                field analysis by magnetic moment method with adoptive cross approximation”, IEEE
                Trans. Magn., Vol. 43 No. 4, pp. 1277-80.
         Tozoni, O.B. and Mayergoyz, I.D. (1974), Computation of the Three-dimensional Electromagnetic
                Fields, Tehnika, Kiev.
         Wang, J.J.H. (1991), Generalized Moment Method in Electromagnetics, Wiley, New York, NY.

         About the authors
         Kazuhisa Ishibashi graduated from Electrical Engineering Department of Tokyo Metropolitan
         University and received PhD degree in 1986. From 1967 to 1992, he was with R&D of
         Toyo Seikan. Since 1992, he had been with Mechanical Engineering of Tokai University as a
         Professor and retired in 2010. Now, he is currently visiting ABB Corporate Research. His main
         interest is computational electromagnetics. Kazuhisa Ishibashi is the corresponding author and
         can be contacted at:
             Zoran Andjelic is Senior Principal Scientist at ABB Corporate Research in Baden,
         Switzerland. There he supervises research and software development projects related to
         electromagnetic and coupled problems in power engineering devices and apparatus. One of his
         major contributions is the synthesis of 3D simulation and automatic optimization in the design
         process of power transformers and switchgears.

         To purchase reprints of this article please e-mail:
         Or visit our web site for further details: