A new BEM technique for nonlinear 2D magnetostatics

Engineering Analysis with Boundary Elements 26 (2002) 795–801 www.elsevier.com/locate/enganabound

A new BEM technique for nonlinear 2D magnetostatics J. Lobry* Department of Electrical Engineering, Faculte´ Polytechnique de Mons, Bd Dolez, 31 B7000 Mons, Belgium Received 10 December 2001; revised 8 May 2002; accepted 12 May 2002

Abstract The paper presents a new boundary element technique for the solution of nonlinear 2D magnetostatic problems. It is based on a recurrent formulation of Poisson’s equation where the change in material property between the domains is modelled as a single layer along the interfaces. Then the nonhomogeneity appears as a source term. The nonlinear regions are decomposed into subdomains where the material is assumed homogeneous. The underlying principle lies on the transmission-line modelling (TLM) method that is commonly used for electric network analysis. A new numerical scheme is derived where the solution is obtained iteratively without any system of algebraic equations for open boundary problems. This is an important point with respect to the classical FEM and BEM treatment that always require the assembly and the solution of a set of simultaneous algebraic equations. Numerical results show significant reduction in memory requirements for our method due to the absence of any matrix storage. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: Nonlinear problems; Transmission-line modelling; Finite elements; Magnetostatics

1. Introduction Magnetostatic problems are usually solved by the finite element method due to its flexibility to deal with nonlinear materials. However, it is necessary to discretise unbounded regions or narrow air gaps that is time-consuming. These drawbacks have triggered research interests in nonlinear boundary element techniques since the early eighties. Seminal papers of Lean and Bloomberg [1], Shao and Zhou [2] are the first attempts but they require domain discretisation of the ferromagnetic regions so that the advantages of BEM are partially lost. Yamada and Wrobel [3] overcome this difficulty by using the dual reciprocity method. But the common bottleneck of all those numerical methods stands in the assembly and the solution of a system of algebraic equations so that a large amount of memory is required. In this paper, a new numerical scheme based on a BEM-type method and the transmission-line modelling (TLM) method is proposed in order to avoid this drawback. The present study follows previous work developed by Lobry et al. [4,5] by considering two-dimensional nonlinear problems. The use of the TLM method to solve nonlinear lumped networks has been introduced by Johns and O’Brien twenty years ago [6]. The general principle is to insert transmission * Tel.: þ32-65-374186; fax: þ32-65-374189. E-mail address: [email protected] (J. Lobry).

lines in order to separate nonlinear elements of a network from the rest of this network for some computational reasons. Surprisingly, this idea can be exploited here with interest for the BE solution of static field problems. Indeed, an electric analogue circuit may be assigned to 2D Poisson’s problem where infinitesimal resistors hold the material property of the problem. In piece-homogeneous problems the property evolves from one region to another. If the elementary resistors are put at the end of elementary links having all the same characteristic impedance, it is shown that the formulation of Poisson’s problem changes in such a way that it relates to a fictitious homogeneous medium. The real material property appears as a single layer density on the interfaces under the form of an infinite summation added to the effective excitation. The solution of the problem is then obtained in an iterative way due to the infinite sum. The boundary element treatment of an open boundary problem with this new scheme reveals that potential value at any point can be computed directly, without any system of equations to be solved. In the case of nonlinear materials, the domains are no longer homogeneous due to nonlinearity so that they must be decomposed in small internal cells that are subdomains where homogeneity may be assumed. Our technique can be applied as explained above but now the terminal resistors outlined above are nonlinear. The nonlinearity is treated separately for each resistor. The subdomains we consider are triangular elements just as in

0955-7997/02/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 9 5 5 - 7 9 9 7 ( 0 2 ) 0 0 0 4 3 - 7

796

J. Lobry / Engineering Analysis with Boundary Elements 26 (2002) 795–801

the finite element method. The difference scalar potential (DSP) formulation of Mayergoyz et al. [7] is used to describe the magnetic phenomena. Numerical simulation allows the comparison of the TLM scheme with the standard FEM technique.

2. DSP-TLM formulation for nonlinear 2D magnetostatics Two-dimensional nonlinear magnetostatics are described by the difference scalar potential (DSP) that was proposed in Ref. [7] as an alternative to the reduced and total potential formulations for 3D problems. The advantage of DSP is that it avoids numerical cancellation error in ferromagnetic regions and the potential is continuous across the problem domain. We do not use the magnetic vector potential for a bad behaviour in convergence was observed with the TLM technique. The DSP formulation involves a two-step procedure and is valid for single-connected iron regions. Suppose a 2D ferromagnetic (iron) domain Vfer of permeability m is field dependent. The region is plunged into an open boundary air region V0 of permeability m0. A current density flows normally into the latter domain. The first step consists of solving the problem in the case of infinite permeability in the iron region. A source field HS is primarily defined as the one generated in the free space by the excitation. A resulting magnetic field H0 is then computed taking into account the presence of the material with infinite permeability. This field vanishes in the iron part. For the second step, the actual permeability is used and a scalar potential c is introduced so that the effective magnetic field is: H ¼ H0 2 grad c:

ð1Þ

This potential is governed by the Poisson’s equation: div m grad c ¼ 2m0 H0n d›Vfer

in V ¼ Vfer < V0

ð2Þ

where H0n is the normal component of H0 along the air –iron interface ›Vfer and d›Vfer is a Dirac distribution defined along the interface. Eq. (1) is also the result of the application of Kirchhoff’s current law to an equivalent distributed-parameter electric circuit as shown in Fig. 1.

Fig. 1. Analogue circuit of Poisson’s problem.

Indeed, the x– y components of grad c may be seen as the voltage drops per unit length along the coordinate axes. Thence m grad c gives the components of the current flowing through the incremental resistors per unit width in V. Taking the divergence results in the locally injected current density m0H0n times d›Vfer. Boundary conditions appear as nodes with fixed voltage with respect to the ground for the Dirichlet condition or open-circuited nodes for the Neumann condition. Suppose that the distributed nonlinear resistors of the circuit of Fig. 1 are put at the end of lossless lines. Their respective surge impedance is either equal to m21 f dx=dy along the x-axis or m21 f dy=dx along the y-axis where mf is a fictitious permeability that is assumed to be constant over V. The transit time T is the same for all the lines. From elementary transmission-line theory [8], Eq. (2) becomes: " # ! N k X Y div mf grad cN þ 2 2 KN2l grad cN2k k¼1

l¼1

¼ 2m0 H0n d›Vfer or, otherwise stated: ! N k X Y m0 7 cN ¼ 2 H d 22 div 2 KN2l grad cN2k ; mf 0n ›Vfer k¼1 l¼1 2

ð3Þ where K is the local reflection coefficient: KN2l ¼

mf 2 mðlgrad cN2lþð1=2Þ lÞ : mf þ mðlgrad cN2lþð1=2Þ lÞ

The value of K depends on the one of lgrad cl computed from the x– y voltage drops at the end of the lines. Remark that from Eq. (1) lgrad cl is equal to the magnitude of the magnetic field H in nonlinear regions. Eq. (3) is the socalled TLM formulation, it can be seen as the substitution of a fictitious material of constant property mf to the effective one. The real property m is modelled as an excitation that can be deduced in an iterative way from preceding known potential values. For physical reasons, Eq. (3) converges towards the expected solution of Eq. (2). Moreover, convergence is unconditional since we are dealing with a circuit made of passive resistors and lines. It is evident that the formulation can be employed for modelling any problem described by Poisson equation. Suppose the nonlinear region Vfer is divided into internal ND cells that are sufficiently small to exhibit a constant reluctivity even though it is field dependent. Those cells may be taken as the elements of a triangular discretisation likewise in the finite element method. The problem is assumed to be unbounded. Let us denote ›Vj the respective boundary of each subdomain Vj. Using our transmissionline model, the material property of Vj is replaced by a fictitious one so that the global domain V is defined by the permeability mf. Henceforth the reflection coefficient K is

J. Lobry / Engineering Analysis with Boundary Elements 26 (2002) 795–801

piecewise defined. Let Kj be the related value in Vj, we have m 2 m0 K0 ¼ f ; mf þ m0 ð4Þ mf 2 mðlHj;N2kþð1=2Þ lÞ ðj ¼ 1; …; ND Þ; Kj;N2k ¼ mf þ mðlHj;N2kþð1=2Þ lÞ where Hj is the average value of the magnetic field on the element Vj. The coefficient K is discontinuous across the interfaces ›Vj so that the divergence terms in the right-hand side of Eq. (3) are singular. Then Dirac distributions appears and Eq. (3) can be rewritten as: "N ! N k D Y X X m0 2 7 cN ¼ 2 H d 22 2 Kj;N2l mf 0n ›Vfer k¼1 j¼0 l¼1  ›cN2k    d þ ›n ›V ›Vj

k Y

!

þ

ð5Þ

3. A numerical BEM –TLM scheme

Let us now derive an integral formulation in the same way as for the classical BEM technique [9]. Introducing the standard 2D Green function wi, some classical integration by parts over domain V leads to: ! þ ›wi ›c N m þ 2 wi dG þ 0 ci;N ¼ 2 cN wi H0n dG ›n ›N mf ›V

þ

›Vfer

k Y l¼1

j¼0

! 2 Kj;N2l

l¼1

! 2 KN2l

j¼1

þ

wi

›cN2k dG ›n

l¼1

wi

›V j

ð ›cN2k wi 72 cN2k dV dG þ ›n Vj

!# :

V0

þ

¼2

wi

 ›cN2k ›w 2 cN2k i dG; ›n ›n

ð8Þ

›V0

where ci is a geometric factor with its usual meaning. On the other hand, using the divergence theorem, the sum of boundary and domain integrals on the subdomains of Vfer reduces to: ð þ ›c wi N2k dG þ wi 72 cN2k dV ›n Vj ¼2

wi 72 cN2k dV

ð Vj

ð9Þ

grad wi grad cN2k dV:

Eqs. (8) and (9) are then introduced in Eq. (7) so that we obtain the final form of the TLM integral formulation: m ð ci;N ¼ 0 w H dG mf ›Vfer i 0n " ! N þ X ›wi k dG 2 ci ci;N2k ð2K0 Þ þ2 cN2k ›n k¼1 ›V0

# V

þ

In order to obtain an economical algorithm, the air and ferromagnetic subdomains must be treated separately. On the one hand, considering the air region, the use of the third Green identity leads to: ð wi 72 cN2k dV þ ci ci;N2k

›Vj

ð

! 2 Kj;N2l

›Vj

3.1. Mixed integral expression

k¼1

ND k Y X

ð7Þ

where ›/›n relates to the inward normal to ›Vj. The summation in the right-hand side represents a single layer density on ›Vj. Eq. (5) is the DSP-TLM formulation of our problem. It should be solved for N ¼ 1 until convergence. The convergence behaviour of this iterative scheme is the same as the transient on a mismatched lossless line supplied with a step function voltage. The choice of the fictitious permeability mf has a strong influence on the number of iterations for convergence within a given accuracy criterion. This problem is discussed in Section 5. Although the numerical treatment has not yet been discussed, a discretisation is implicit in the formulation due to the nonlinear material. This apparent anticipation is necessary for what follows.

þ2

›V0

l¼1

j

"N N k D X Y X

potential and its normal derivative vanish. Therefore, the associated boundary integral of the right-hand side is zero. Now, the domain integral of the Laplacians may be split into a summation over the subdomains Vj. Recall that the first term ðj ¼ 0Þ in this summation relates to the linear domain V0 so that the recurrent product of related reflection coefficients reduces to ð2K0 Þk : Under the considerations outlined, Eq. (6) can be written as: m ð ci;N ¼ 0 w H dG mf ›Vfer i 0n " ! N þ ð X ›cN2k k 2 dG þ ð2K0 Þ þ2 wi wi 7 cN2k dV ›n V0 k¼1

#

2 KN2l 72 cN2k

797

2 ð6Þ

for any internal point in V. As we deal with open problems, the exterior boundary ›V is put at the infinity where both

ND k X Y j¼1

l¼1

! 2 Kj;N2l

ð Vj

# grad wi grad cN2k dV

ð10Þ

Expression (10) is an iterative prescription where the potential value at any point i and at TLM step N is explicitly calculated from the values on the interfaces at

798

J. Lobry / Engineering Analysis with Boundary Elements 26 (2002) 795–801

preceding steps. Therefore it can be easily predicted that no system of equations is to be expected after the discretisation as demonstrated in Section 3.2. 3.2. BEM discretisation As outlined in Section 2, the geometric discretisation is implicit in Eq. (10) since the subdomains Vj are chosen to manage the nonlinearity. The discretisation of ›V0 is then automatically obtained so that it is decomposed into NG boundary elements Gp. Next, an approximation must be adopted for the scalar potential c on the triangular elements. In order to be consistent with the objective of the discretisation of Vfer, a linear variation is considered so that the flux density is constant over the triangles Vj. The piecewise homogeneity is then naturally satisfied. Thus, after discretisation, Eq. (10) leads to

ci;N ¼ bsi þ 2

N X

(

" ð2K0 Þk

4. Performances of the method NG X

ð2Þ ðhð1Þ ip cp1;N2k þ hip cp2;N2k Þ

p¼1

k¼1

# 2 ci ci;N2k 2

ND k Y X j¼1

!

)

2 Kj;N21 Gij ðgrad cN2k Þj

l¼1

ð11Þ where bsi is the computed contribution of the source term and the hip and Gij factors results from the evaluation of the boundary and domain integrals. These last integrals are generally computed in a numerical manner by using Gaussian quadrature rules but analytical computation is preferable when the point i is one of the vertices of the element. Finally, a better form of Eq. (11) for computer implementation is the recurrence:

ci;N ¼ bsi þ

NG X

ð2Þ ðhð1Þ ip bp1;N þ hip bp2;N Þ 2 ci bi;N 2

p¼1

( with

ND X

Gij Hj;N

j¼1

bi;0 ¼ 0

bi;N ¼ 2K0 ðbi;N21 þ 2ci;N21 Þ

Hj;0 ¼ 0

Hj;N ¼ 2Kj;N21 ½Hj;N21 þ ðgrad cN21 Þj 

j ¼ 1; …; 1:

receiving end field values normally converge towards each other. Remark that the nonlinearity is solved for each element, separately.Eq. (12) requires the computation of the potential values at the discretisation nodes wherever point i may be. Thence those nodes must always be included in the set of the computation points. The computation at any external node (in V0) is carried out from the potential values obtained in Vfer, after the convergence of the BEM – TLM algorithm. Finally, the choice of the fictitious permeability mf is important in the TLM method. It has been observed experimentally that there exists an optimal value mf,opt that allows a minimum of iteration steps. This value can be found from experience or dynamically by adapting mf during the TLM iterations.

ð12Þ

As pointed out previously, no system of equations is involved for the calculation of the potential value ci at a given point i in V with our numerical scheme. This important result was already obtained in previous papers for linear homogeneous domains [5]. The reflection coefficient Kj;N21 in the nonlinear region is field dependent and should be computed iteratively from Eq. (4) with k ¼ 1: This leads to an implicit scheme since the permeability m relates to the value of the magnetic field at the receiving end of the transmission lines. This can be avoided if we take the field Hj;N21 at the sending end that is a known value. This approximation does not alter the result since the sending and

As it has already been emphasised, the main advantage of the BEM –TLM scheme (12) is that no system of algebraic equations is involved for open boundary situations so that some gain in memory requirements are to be expected compared to standard methods. Classically, a nonlinear 2D magnetostatic problem can be solved by using a FEM code. Then all the regions (air, copper, iron) have to be discretised. The infinite air region is usually replaced by a bounded box whose boundary is generally assumed to be at zero magnetic vector potential. The relative size of the meshes of the regions depends on the geometry of the problem. Let N0 and Nfer be the number of nodes in the air and ferromagnetic regions, respectively, the ratio a ¼ N0 =Nfer is generally greater than one. A FEM code requires at least the vector of nonzero elements of the stiffness matrix. Any node in a Delaunay triangular decomposition is connected to about six others so that, taking the symmetry of the matrix into account, the size of the vector mentioned above is of the order 4ðN0 þ Nfer Þ ¼ 4ða þ 1ÞNfer : The memory requirements are still greater than that due to the workspace needed by e.g. an ICCG solver. The implementation of the TLM algorithm as given by Eq. (12) only requires up to 2ND þ Nfer for the storage of the bi terms and of the x –y components of the Hj terms in Eq. (12). This nearly equals 5Nfer since the number of triangles is about twice the number of nodes in a two-dimensional triangular mesh. Thus a gain of more than 4(a þ 1)/5 is to be expected. In the case of a FEM – BEM algorithm, the air region is not meshed but a nonsymmetric dense block appears in the stiffness matrix due to the boundary element discretisation. It can be shown that our TLM algorithm remains interesting in this situation. A drawback of our technique is that CPU time grows as N 2 whereas it is N log N with FEM-type algorithm due to conjugate gradients.

J. Lobry / Engineering Analysis with Boundary Elements 26 (2002) 795–801

799

5. Numerical results The practical performances of our method are presented by comparing with the standard FEM. Both the algorithms have been implemented in Fortran 77 on a Pentium III1 GHz. The FEM code is classical and uses the Newton – Raphson iteration process to solve the nonlinearity with a convergence criterion based on the magnetic energy. An error less than 1026 is considered in both the FEM and TLM algorithms. The ICCG technique solves the Jacobian system in the FEM code. When comparing the numerical results for a given example, the discretisation of the ferromagnetic region is supposed to be the same for both the simulations. Three examples have been treated. The first one consists of a ferromagnetic cylinder in a uniform magnetic field where an analytical solution is available. The second case treats the problem of a shielded transmission line where magnetic energy is confined inside a close domain. Finally, a more practical example is presented. 5.1. Ferromagnetic cylinder in a uniform magnetic field The situation is depicted in Fig. 2a together with an optimal mesh (Fig. 2b). The saturation curve is the same for all the examples presented after and is plotted in Fig. 3.

Fig. 2. Ferromagnetic cylinder in a uniform magnetic field: (a) situation and (b) optimised mesh.

Fig. 3. Magnetic flux distribution of the problem. B is constant inside the cylinder.

The flux lines are displayed in Fig. 4 as an illustration. In the FEM discretisation, the air boundary consists of a 1 m £ 1 m square box with Dirichlet and Neumann conditions chosen so that the flux density B0 is indeed of 1 T in the x-direction. In the BEM –TLM algorithm, a source field HS is forced to 1/m0 A/m. We intend to compare the respective accuracies of the FEM and TLM algorithms for an increasing refinement of the discretisation. The magnetic energy inside the cylinder is taken as the representative value. The numerical results are shown in Fig. 5 with the analytical value. It is noticed that the two TLM and FEM treatments give upper and lower bounds for energy, respectively, with a clearly better convergence for the TLM algorithm. This is not surprising since BEM is intrinsically better suited to open problems than FEM. Moreover, the necessity of meshing a large free space domain in the FEM code makes the BEM – TLM method attractive. The ratio a discussed in Section 4 is about 30 if an optimised mesh is considered. This rather high value is due to a constant value of magnetic field inside the cylinder so that a fine mesh in not required except along its edge (Fig. 2b). The expected gain in memory requirements is at least equal to 25.

Fig. 4. B –H curve of the iron.

800

J. Lobry / Engineering Analysis with Boundary Elements 26 (2002) 795–801

Fig. 5. Convergence of the magnetic energy inside the cylinder vs. the size of the problem.

Fig. 7. Convergence of the magnetic energy inside the shield vs. the size of the problem.

5.2. Shielded transmission-line This example consists of the computation of the magnetic field around a parallel bar transmission-line enclosed in a ferromagnetic hollow cylinder so that the magnetic field is mainly confined inside the shield due to the high value of the permeability (Fig. 6). In a classical FEM

situation this case can be treated as a closed problem with the A ¼ 0 Dirichlet condition on the outer boundary of the cylinder. Fig. 7 shows a similar accuracy of the two algorithms. Although the air domain is small, a fine discretisation is required due to a rapidly varying field in this region. Then the a value is high and a gain of about 20 is

Fig. 6. Shielded parallel bar transmission-line: (a) situation and (b) field pattern.

Fig. 8. C-shaped yoke inductor: (a) situation and (b) field pattern.

J. Lobry / Engineering Analysis with Boundary Elements 26 (2002) 795–801

801

size of the problem. This feature was already observed in previous papers, e.g. Ref. [5].

6. Conclusions

Fig. 9. Convergence of the magnetic energy inside the yoke vs. the size of the problem.

reported for an optimised discretisation obtained from an adaptive meshing algorithm.

In this paper, a new BEM iterative scheme has been developed for nonlinear 2D magnetostatics. It is based on the substitution of fictitious excitations to the change in material property across the interface between different media. This approach leads to the absence of system of algebraic equations to be solved that is a new feature in finite element-type methods. Therefore, there is much to be gained by way of economy in memory requirements. This may be interesting for the treatment of large problems or when the available random access memory in a computer is not sufficient. The main drawback of the technique is the computation time growth. Future work should be oriented to the validation of the technique with other examples and its extension to three-dimensional problems where magnetic scalar potential models are normally used.

5.3. Ferromagnetic yoke inductor with air-gap Finally, consider the C-shaped steel yoke inductor presented in Fig. 8. In the FEM discretisation, the air boundary is placed at about 15 device diameters away. Here again, the accuracy for both methods is similar as shown in Fig. 9. The ratio a (about 6) is lower here than in the previous examples because of a more complicate shape of the yoke. Then a gain of at least 5 is to be expected. 5.4. Discussion of the numerical results The examples presented above show that the main advantage of our BEM – TLM method lies in a generally high gain in memory amount for solving a given problem with respect to the common FEM algorithm. The accuracy is also slightly better since the open region is considered as it is due to the BEM though it is not the case with FEM where a bounding box is required. However, the reported execution times for TLM technique are much higher than the one of the classical algorithm. This drawback was to be expected as the computation time in the BEM assembly generally grows with N 2. Finally, we note that the number of TLM iterations is generally about a few hundred depending on the problem and is quasi-independent of the

References [1] Lean MH, Bloomberg DS. Nonlinear boundary element method for two-dimensional magnetostatics. J Appl Phys 1984;55(6):2195 –7. [2] Shao KR, Zhou KD. The iterative boundary element method for nonlinear electromagnetic field calculations. IEEE Trans Magn 1988; 24(1):150– 3. [3] Yamada T, Wrobel LC. A new approach to magnetic field analysis by the dual reciprocity boundary element method. Int J Numer Meth Engng 1993;36:2073– 85. [4] Lobry J, Broche C, Tre´cat J. Use of the transmission-line modelling in BEM for solution of piece-homogeneous static field problems. IEE Proc Sci Meas Technol 1996;143(3):157 –62. [5] Lobry J, Broche C, Tre´cat J. Direct BEM solution of the open boundary Poisson’s problem with the TLM method. In: Brebbia CA, editor. Boundary elements XVIII. Southampton: Computational Mechanics Publications; 1996. p. 183 –92. [6] Johns PB, O’Brien M. Use of the transmission-line modelling (TLM) method to solve nonlinear lumped networks. Radio Electron Engng 1980;50(1/2):59 –70. [7] Mayergoyz ID, Chari MVK, D’Angelo J. A new scalar potential formulation for three-dimensional magnetostatic problems. IEEE Trans Magn 1987;23(6):3889 –94. [8] King RWP. Transmission-line theory. New York: McGraw-Hill; 1955. [9] Brebbia CA, Dominguez J. Boundary elements—an introductory course. Southampton/New York: Computational Mechanics Publications/McGraw-Hill; 1992.