The volume integral equation method in magnetostatic problem Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

Discrete & Continuous Models 2019 27 (1) 60—69

. „ . http://journals.rudn.ru/miph

& Applied Computational Science

Research article

UDC 519.644

DOI: 10.22363/2658-4670-2019-27-1-60-69

The volume integral equation method in magnetostatic problem

Pavel G. Akishin, Andrey A. Sapozhnikov

Joint Institute for Nuclear Research 6 Joliot-Curie str, Dubna, Moscow region, 1^1980, Russian Federation

(received: February 22, 2019; accepted: October 21, 2019)

This article addresses the issues of volume integral equation method application to magnetic system calculations. The main advantage of this approach is that in this case finding the solution of equations is reduced to the area filled with ferromagnetic. The difficulty of applying the method is connected with kernel singularity of integral equations. For this reason in collocation method only piecewise constant approximation of unknown variables is used within the limits of fragmentation elements inside the famous package GFUN3D. As an alternative approach the points of observation can be replaced by integration over fragmentation element, which allows to use approximation of unknown variables of a higher order.

In the presented work the main aspects of applying this approach to magnetic systems modelling are discussed on the example of linear approximation of unknown variables: discretisation of initial equations, decomposition of the calculation area to elements, calculation of discretised system matrix elements, solving the resulting nonlinear equation system. In the framework of finite element method the calculation area is divided into a set of tetrahedrons. At the beginning the initial area is approximated by a combination of macro-blocks with a previously constructed two-dimensional mesh at their borders. After that for each macro-block separately the procedure of tetrahedron mesh construction is performed. While calculating matrix elements sixfold integrals over two tetrahedra are reduced to a combination of fourfold integrals over triangles, which are calculated using cubature formulas. Reduction of singular integrals to the combination of the regular integrals is proposed with the methods based on the concept of homogeneous functions. Simple iteration methods are used to solve non-linear discretized systems, allowing to avoid reversing large-scale matrixes. The results of the modelling are compared with the calculations obtained using other methods.

Key words and phrases: finite element method, magnetostatics, volume integral equations, systems of nonlinear equations, cubature formulae, iterative process

Introduction

The paper [1] gives an overview of existing methods and programs for the numerical simulation of magnetic systems. The issues concerning the volume

© Akishin P. G., Sapozhnikov A. A., 2019

This work is licensed under a Creative Commons Attribution 4.0 International License http://creativecommons.org/licenses/by/4.0/

integral equation method [2-5] for magnetic system calculations are discussed in this article. We consider the approach, which allows to construct high order accuracy approximation method for initial problem discretization. The main steps of the modelling process are being discussed: discretization of the initial equations for linear approximation of the magnetic field, algorithm of area division into elements, matrix element calculation procedures, a method of solving a corresponding system of nonlinear equations. The numerical results obtained for the modelling of quadrupole and dipole magnet agree with those obtained by the famous code TOSCA [6].

1. Volume integral equation method

Let B(a), H(a), M(a) be the induction, the intensity and the magnetization of the magnetic field at the point a. The values B, H, M are connected by the following nonlinear ratios [7]:

H(a) = B(al^ , M(a) = ^ - H(a), (1)

y {IB(a)l) ^o yo

where y0 is the absolute magnetic permeability of vacuum; y(x) is the magnetic permeability. The following integral equation is true:

H(a) = HS (a) + 4n V«/(M(x), Va ix—aj) dvx,

G

where HS(a) is the magnetic field from the current winding. G is the area filled with iron. The field HS(a) can be found according to the Biot-Savart law [8]:

1 f J(x) H S (a) = - Ro^ J-L dvx,

R3

where J(x) is current density. The difficulty of applying the integral approach is related to the singularity of the kernel of the integral equations.

This is the reason why only a piecewise approximation of unknown parameters within division element area is used in GFUN3D code [2]. Alternatively to collocation method, integration over dividing elements can be used.

It allows to use higher order approximation for unknown variables. The most convenient mathematical approach for constructing such type of approximations is the finite element method (FEM) [9-12].

Let us divide area G into tetrahedrons {G^}. We suppose that the frag-

N

mentation G = (J Gi satisfies the requirements of FEM. Let us assume

i= 1

{Pk, k = 1,... ,L} is the set of all vertexes in all tetrahedrons {Gi}. Let us introduce the notation Bk = B(Pk), Hk = H(Pk), Mk = M(Pk). We denote fk(x) — as a node function, associated with vertex Pk. The functions fk (x) on each tetrahedron are linear functions, equaled to 1 at the vertex Pk and to 0

at any other vertexes. Using these notations we define linear approximations for vectors B(a),H(a), M(a):

L L L

B(a) = Y fk(a)Bk, Hi(a) = Y fk(a)#k, M(a) = Y fk(a)Mk.

k=i

k=l

k=l

We characterize a discretized formulation of magnetostatic problem, using the finite element linear approximation within division elements [13]:

Y] fi(a)fj(a)Hjdv-a = j fi(a)HS(a)dv-a+

j=l G

L

G

+ Y / fi(a)

j=l G

Va

4n

G

fj (a) Mj , Va

1

|x — a|

dvx

dva, i = 1, L. (2)

Let matrix [C] be a matrix of [3L x 3L]:

"[Cii] ■■■ [Cil]

[C ] =

L[ClI] ■■■ [Cll]J

where [Cj] — is a diagonal matrix of [3 x 3] dimension such as

[Cij] = J fi(a)fj (a)dva

G

(1 o o^ 0 1 0 \o 0 1/

Let [A] be a matrix of [3L x 3L] dimension:

'[Aii] ■■■ [Ail]'

[A]

.[Ali] ■■■ [All],

where [Aij-] — is a matrix of [3 x 3] dimention, such as for any constant vector M the following ratio is true:

[Aij ]M

fi(a)dva

G

fj (X) Ma, Va

1

G

|x — a|

dvx

(3)

Let the following be true:

B= {Bi,B2 ,...,Bl)t ,

M(B) = (iiaM {Bi),VoM (B2),. . (BL)f,

T

HS = ( io J fi(a)HS (a)dva ,ioJ f2(a)HS (a)dva ,...,ioJ fl(o)Hs (a)dva

G G G

Taking into account (1), the system (2) could be written in the following way:

[C]B = HS + ([A] + [C]) M(B). (4)

Using node functions of higher order similarly to (2) it is possible to formulate a discretization with quadratic, cubical or higher approximation of variables within the element.

2. Finite element mesh generating

In order to build a discretization for the integral equations (2) the region of calculations should be divided into the tetrahedrons, satisfying the requirements of FEM. There are many articles dedicated to this problem [14-17]. Depending on the task, there are different requirements for mesh elements. In subregions, where the solution should change faster, more detailed discretization is needed, and as a result, the element size must be smaller. And, vice versa, within the regions of slow solution changing, detailed division leads to a big number of elements, thus complicating the solving of the final dis-cretized system of equations. That is why in such regions the element size must be larger. Moreover, the final mesh elements should not be degenerated. The degenerated elements affect the solution approximation and the convergence of iterative methods, used for solving the discretized problem. In the case of constructing mesh with heterogeneous materials there are additional restraints for a division related to medium edge borders. These borders usually look like curves on a plain or surfaces in space and should not be crossed by the edge of the mesh. In fact, these limitations show that every element of the mesh should consist of only one particular substance. The region edges can be approximated with lines or surfaces as well as with curves and second or further order surfaces. Such edges are approximated with resultant elements according to the size of the considered mesh sub-region.

In report [18] the problems of multidimensional finite element mesh generation based on electromagnetic fields modelling in large-scale electromagnetic machines have been discussed. A description of the three-dimensional adaptive mesh generator 3DFEMMesh, a review of mesh generation methods and a number of the existing criteria for a quality assessment of the constructed mesh are given in this work. The procedure of mesh construction is based on a representation of the problem domain as a combination of standard 3D macroblocks. After generating of a two-dimensional mesh on all macroblocks boundaries, three-dimensional mesh in each block may be constructed individually. The program has a graphical interface for the data entry and a visual assessment of the partition quality. The generator 3DFEMMesh is included into JINR programme library JINRLIB [19].

3. Matrix element calculation

The problem of defining matrix coefficients of the discretized equations can be reduced to calculating sixfold, singular in general case, integrals from (3) by two tetrahedrons. The simplest way to evaluate these integrals is to use a cubature. Given a big number of integrals to be calculated, requirements to cubature formula optimality are extremely important. Because of singularity of the function being integrated the necessity of using cubature formula of high accuracy arises [20]. In monograph [21] the issues of the general theory of cubature formula building are studied and many cubature formulas for different types of simplexes are listed. While calculating integrals, which are necessary for the discretized equations, there might be situations when the integrated function is singular. Moreover, there are situations when the function being integrated is singular in every point of the volume under integration. In those cases cubature formula application is not possible and we need to develop further methods for such kind of integration.

First of all, we consider the method which allows to decrease the coefficient calculation time of the discretizated systems. The matrix coefficients in [Aij] from (3) can be presented as sum of integrals such as:

jn,j,l -if

Gm Gn

Taking into account that {fk(x)} are linear functions and, as a consequence, the function gradient vectors are constant inside each tetrahedron, such volume integrals can be reduced to surface integrals thus:

fi(x)fj(a)

d2

1

dxk8al |x — a!

dvrdvn

jn,j,i _ f I fi(x)fj(a)(dSx,ek)(dSa,ei)

Jm'i'k —ff |x — a|

dGm dG n

_05dfm(x) ja) l l (dSx,dSa)

dxk dal J J |x — a|

dGm dG n

f & I f fj(a)(dSa,ei)((x ,— ■a)adSx)

dxk J J |x — a|

dGm dG n

n ç;d/n(a) / / _ ,((a — x),dSa)

— f f fi(x)(dSx,ek) — , (5)

dGm dG

where ek, el are unit coordinate system vectors.

The notation ôfim(x)/ôxk means that the derivative is calculated on the tetrahedron Gm. it is important to note that region G consists of a union of tetrahedrons. Then the borders {dG} are triangles. Thus, calculating the expressions (5) reduces a 6D integral to a sum of 4D integrals over two triangles. There are four types of the positional relationships of triangles in space: triangles do not cross, triangles have one mutual vertex, triangles have one mutual side and triangles are congruent. For the first type the cubature

formula application is possible. For the others it is not possible due to the singularity of the expressions being integrated. Let us note that the expression under integration in (5) can be represented as a sum of homogeneous functions. Let the function f (x) be a homogeneous function if for any A > 0

f (Ax) = Ak f (x).

Let us illustrate the method of integration of homogeneous functions from [22] taking as an example the integral:

j __dSx dSy

° |x - V1

S1 S2

where Si, S2 are shown in Figure 1.

Let S1 be a triangle AB'C', obtained by stretching the triangle S1 in A times with respect to point A. Similarly S2 is a triangle AD'F', obtained by stretching triangle S2 in A times with respect to point A. Let J (A) be

j (a) __dSx dSy

|x - v\'

Si S2

Substituting variables x = Axi = Avji integral J (A) can be reduced to:

J (A)=A3 // if-f. (6)

Si S2

Let rTi be trapezium B/BCCf, TT2 — trapezium DDfFfF. Let us calculate the limit of difference ratio:

dJ (A)

dA

_ lim (J(A\ - J(1))

x^i (A - 1)

Using the additivity of integrals as a set function over which they are taken we have:

dJ (X)

dX

1

A=1

lim

X - 1

[u+u+u]

\Ti S2 Si T2 Ti T2 J

dSXl dSy1 |xi - ?7i1

(7)

From (7) we have:

dJ (X)

dX

hi

dlX1 dSyi + h

J |xi - yih '"2 } J |xi - yi1

BC S2 Si DF

dSxi dlyi

A=i

(8)

where hi — height LA, h2 — height AM. Differentiating (6) by X and taking into account (8), we get:

Jo

J(1) = 3

hi

BC S2

dlXi dSyi T-Z-=7 + h2

- y i1

Si DF

dSXi dlyi

|xi - yi|

(9)

It must be mentioned that the expressions on the right hand of (9) are regular integrals. It is possible to use the cubature formula for calculations of these integrals.

After application of this method all singular integrals from (5) can be reduced to the superposition of regular integrals of lower order, for calculation of which a cubature formula can be used. In the case when the triangles have one common side or coincide, this approach should be applied successively twice and three times respectively. The method of integration illustrated above allows to reduce all singular integrals from (5) to superpositions of regular integrals of lower order, the calculation of which can be done by cubature.

4. Iterative methods for nonlinear systems solving

In practice to achieve the demanded approximation accuracy it is necessary to split region G into smaller elements, that leads to huge dimension rise of the nonlinear discretized systems of the equations. It is extremely difficult to apply methods which require inversion of high order matrixes. So, for solving discretizated systems of equations (4) simple iterative process is used:

[C]Bk+i = (HS + ([A] + [C]) M(Bk)),

B0 = 0, k =1, 2,....

This process will be finished when equation residuals become less than the previously set value e. For solving linear systems of equations [C]x = y the incomplete Kholetsky expansion method in combination with the conjugate gradient method are used [23].

5. Magnet system modelling

The method of volume integral equations with the linear approximation of magnetization has been used for modelling the dipole and quadrupole magnets. The model of a variant of the projected dipole magnet for CBM experiment (GSI, Darmstadt) is shown In Figure 2a. A splitting of the magnet into the tetrahedrons has been done with the help of the generator 3DFEMMesh. In the process of modelling the dipole symmetry of magnetic field has been taken into account, that allowed to reduce the number of unknown parameters by 8 times. One eighth of the magnet has been divided into 5264 tetrahedrons. There are 1363 vertexes in all tetrahedrons. In Figure 2b the distribution of the magnet field module inside the magnet is shown.The results given in Figure 2c show agreement of the present method with the famous code TOSCA [6] which is based on solving partial differential equations.

a. b. c

Figure 2. 3D modelling of dipole magnet CBM

The proposed approach has been also applied to the modelling of the BOOSTER quadrupole magnet of the projected accelerating complex NICA (JINR, Dubna).

In the process of modelling the quadrupole symmetry of magnetic field has been taken into account. This allowed to reduce the number of unknown parameters by 16 times. One sixteenth of the magnet was divided into 7040 tetrahedrons. The total number of vertexes was 1729. In Figure 3a there is a computer model of the magnet. In Figure 3b the allocation of magnet field module inside the magnet is shown. The resulting comparison obtained by the proposed method and TOSCA code is given in Figure 3c.

Figure 3. 3D modelling of quadrupole magnet BOOSTER

6. Conclusion

The issues of applying the volume integral equation method to the calculations of magnetic systems have been considered in this article. Within the framework of finite element method for discretization of continual equations an alternative approach has been used. This method is based on the substitution of observation points with integration by division elements. Thus, the main problem of applying the volume integral equation method which is related to kernel singularity is removed. The suggested approach allows to increase the order of approximation of the initial equations. Procedures for calculations of matrix coefficients for discretized equations and methods for solving a corresponding non-linear system of equations have been suggested. The results of magnet system simulations based on this approach are in agreement with the calculations by other programs.

References

1. C. W. Trowbridge, J. K. Sykulski, Some key developments in computational electromagnetics and their attribution, IEEE transactions on magnetics 42 (4) (2006) 503-508. doi:10.1109/TMAG.2006.872491.

2. A. G. Armstrong, GFUN3D User Guide. RL-76-029/A, CERN (1976).

3. A. A. Halacsy, Three-dimensional analysis of magnetic fields, in: Proceedings 3rd International Conference on Magnet Technology, Hamburg, 1970, pp. 113-128.

4. M. J. Friedman, Mathematical study of the nonlinear singular integral magnetic field equation I, SIAM Journal on Applied Mathematics 39 (1)

(1980) 14-20. doi:10.1137/0139003.

5. M. J. Friedman, Mathematical study of the nonlinear singular integral magnetic field equation II, SIAM Journal on Numerical Analysis 18 (4)

(1981) 644-653. doi:10.1137/0718042.

6. J. Simkin, C. W. Trowbridge, Three dimensional non-linear electromagnetic field computations using scalar potentials, IEE Proceedings B - Electric Power Applications 127 (6) (1990) 368-374. doi:10.1049/ip-b.1980.0052.

7. J. A. Stratton, Electromagnetic theory, MCgraw-hill, 1941.

8. J. D. Jackson, Classical electrodynamics, 2nd Edition, John Wiley & Sons, 1975.

9. O. C. Zienkiewicz, The finite element method in engineering science, MCgraw-hill, 1971.

10. J. T. Oden, Finite elements of nonlinear continua, McGraw-Hill, New York, 1971.

11. W. G. Strang, G. J. Fix, An analysis of the finite element method, Prentice-Hall, 1973.

12. J. P. Aubin, Approximation of elliptic boundary-value problems, Wiley-Interscience, 1972.

13. P. G. Akishin, A. A. Sapozhnikov, Linear approximation of volume integral equations for the problem of magnetostatics, in: EPJ Web Conferences. Mathematical Modeling and Computational Physics (MMCP 2017), Vol. 173, 2018, Article Number 03001. doi:10.1051/epjconf/201817303001.

14. Z. Cendes, Magnetic field computation using Delaunay triangulation and complementary finite element methods, IEEE Transactions on Magnetics 19 (1983) 2551-2554. doi:10.1109/TMAG.1983.1062841.

15. E. K. Buratynski, A three-dimensional unstructured mesh generator for arbitrary internal boundaries, in: Numerical Grid Generation in Computational Fluid Mechanics: Proceedings, Pineridge Press, Swansea, 1988, pp. 621-631.

16. M. Berzins, Mesh quality: a function of geometry, error estimates or both?, Engineering with Computers 15 (1999) 236-247. doi:10.1007/s003660050019.

17. L. Durbeck, Evaporation: a technique for visualizing mesh quality, in: 8th International Meshing Roundtable, Sandia National Laboratories, South Lake Tahoe, 1999, pp. 259-265.

18. P. G. Akishin, A. A. Sapozhnikov, Automatic generation of three-dimensional grids, JINR, Dubna, 2015.

URL http://www1.jinr.ru/Preprints/2015/058(P11-2015-58).pdf

19. P. G. Akishin, A. A. Sapozhnikov, 3DFEMMesh - program for automatic generation of three-dimensional Mesh.

URL http://wwwinfo.jinr.ru/programs/jinrlib/3dfemmesh/

indexe.html

20. M. Abramowitz, I. Stegun, Handbook of mathematical functions with functions, graphs, and mathematical tables, 1964.

21. I. P. Mysovskikh, Interpolation cubature formulas [Interpolyatsionnyye kubaturnyye formuly], Nauka, Moscow, 1981, in Russian.

22. P. G. Akishin, The integral equation method in magnetostatic problems: abstract of a PhD thesis, Ph.D. thesis, JINR, Dubna, in Russian (1983).

23. J. A. Meijerink, H. A. van der Vorst, An iterative solution method for linear systems of which the coefficient matrix is a symmetric M-matrix, Mathematics of Computation 31 (137) (1977) 148-162. doi:10.2307/2005786.

For citation:

P. G. Akishin, A. A. Sapozhnikov, The volume integral equation method in magnetostatic problem, Discrete and Continuous Models and Applied Computational Science 27 (1) (2019) 60-69. DOI: 10.22363/2658-4670-201927-1-60-69.

Information about the authors:

Pavel G. Akishin (Russian Federation) — Doctor of Physical and Mathematical Sciences, Deputy Head of Scientific Department of Computational Physics, Laboratory of Information Technologies, Joint Institute for Nuclear Research (e-mail: akishin@jinr.ru, phone: +7(49621)62630, ORCID: https://orcid.org/0000-0001-9542-046X, ResearcherID: B-8330-2016, Scopus Author ID: 56780566600)

Andrey A. Sapozhnikov (Russian Federation) — Junior Researcher of Scientific Department of Computational Physics, Laboratory of Information Technologies, Joint Institute for Nuclear Research (e-mail: asap@jinr.ru, phone: +7(49621)62630, ORCID: https://orcid.org/0000-0002-1270-043X, Scopus Author ID: 57200759121)