Discrete and Continuous Models and Applied Computational ScienceDiscrete and Continuous Models and Applied Computational Science2658-46702658-7149Peoples' Friendship University of Russia2219710.22363/2658-4670-2019-27-1-60-69Research ArticleThe volume integral equation method in magnetostatic problemAkishinPavel G<p>Doctor of Physical and Mathematical Sciences, Deputy Head of Scientific Department of Computational Physics, Laboratory of Information Technologies</p>akishin@jinr.ruSapozhnikovAndrey A<p>Junior Researcher of Scientific Department of Computational Physics, Laboratory of Information Technologies</p>asap@jinr.ruJoint Institute for Nuclear Research15122019271606920112019Copyright © 2019, Akishin P.G., Sapozhnikov A.A.2019<p>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.</p>finite element methodmagnetostaticsvolume integral equationssystems of nonlinear equationscubature formulaeiterative processметод конечных элементов, магнитостатика, объемные интегральные уравнения, системы нелинейных уравнений, кубатурные формулы, итерационный процесс<p>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]. 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 the following nonlinear ratios [7]: are connected by H (a) = , B(a) |B(a)|) 0 M (a) = B(a) 0 - H(a), (1) where 0 is the absolute magnetic permeability of vacuum; (x) is the magnetic permeability. The following integral equation is true: ( H (a) = H S (a) + 1 1M (x), 1 \ a x dv , 4 a G |x - a| where H S (a) is the magnetic field from the current winding. G is the area filled with iron. The field H S (a) can be found according to the Biot-Savart law [8]: H S (a) = 1 Rot ( 4 a R3 J(x) |x - a| dvx, 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 param- eters 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 approx- imations is the finite element method (FEM) [9-12]. Let us divide area G into tetrahedrons {Gi}. 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 ), M k = 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 B(a) = fk (a)Bk , k=1 L H (a) = fk (a)Hk , k=1 L M (a) = fk (a)M k. k=1 We characterize a discretized formulation of magnetostatic problem, using the finite element linear approximation within division elements [13]: L ( f (a)f (a)H dv = ( f (a)H S (a)dv + i j j=1 G L j a i a G + ( f (a) a ( 1 f (x) M , a 1 \ dv dv , i = 1, L. (2) i j=1 G j 4 j G |x - a| x a Let matrix [C] be a matrix of [3L 3L]: [C11] [C1L] [C] = ... . . . . .. , [CL1] [CLL] where [Cij ] - is a diagonal matrix of [3 3] dimension such as 1 0 0 0 1 0 0 0 1 ( [Cij ] = G fi(a)fj (a)dva . Let [A] be a matrix of [3L 3L] dimension: [A11] [A1L] [A] = ... . . . . .. , [AL1] [ALL] where [Aij ] - is a matrix of [3 3] dimention, such as for any constant vector M the following ratio is true: ( ( a 1 1 \ 4 [Aij ]M = G fi(a)dva G fj (x) M, a |x - a| dvx . (3) Let the following be true: T B = B1, B2, . . . , BL) , T M (B) = 0M (B1), 0M (B2), . . . , 0M (BL)) , T ( H S = 0 G ( f1(a)H S (a)dva, 0 G ( f2(a)H S (a)dva, . . . , 0 G fL(a)H S (a)dva . Taking into account (1), the system (2) could be written in the following way: [C]B = H S + ([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. 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 require- ments 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 discretiza- tion 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 conver- gence 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 usu- ally 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 sec- ond 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 genera- tion based on electromagnetic fields modelling in large-scale electromagnetic machines have been discussed. A description of the three-dimensional adap- tive 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 individ- ually. 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]. 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 ( ( Г 2 1 l m,i,k = Gm Gn k fi(x)fj (a) x al |x - a| dvxdva. 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: Jn,j,l m,i,k = - fi (x)fj (a)(dSx , ek )(dSa , el ) |x - a| Gm Gn fm(x) fn(a) (dS , dS ) - 0.5 i j x a - xk al Gm Gn |x - a| - 0.5 i fm(x) xk fj (a)(dSa, el) ((x - a), dSx) |x - a| - Gm Gn j fn(a) ((a - x), dSa) (5) - 0.5 al Gm Gn fi(x)(dSx, ek ) , |a - x| where ek , el are unit coordinate system vectors. i The notation fm(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 {G} 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 0 f (x) = kf (x). Let us illustrate the method of integration of homogeneous functions from [22] taking as an example the integral: ( ( J0 = S1 S2 dSxdSy , |x - y| where S1, S2 are shown in Figure 1. Figure 1. Illustration of homogeneous function integration by two triangles Let S1 be a triangle ABtCt, obtained by stretching the triangle S1 in times with respect to point A. Similarly S2 is a triangle ADtF t, obtained by stretching triangle S2 in times with respect to point A. Let J () be J () = ( S1 ( S2 dSxdSy . |x - y| Substituting variables x = x1, y = y1 integral J () can be reduced to: ( ( J () = 3 dSx 1 dSy1 . (6) S1 S2 |x1 - y1| Let T1 be trapezium BtBCCt, T2 - trapezium DDtF tF . Let us calculate the limit of difference ratio: J () (J () - J (1)) =1 = lim 1 . ( - 1) Using the additivity of integrals as a set function over which they are taken we have: J () 1 ( ( ( ( ( ( dSx1 dSy1 = lim + + . (7) =1 From (7) we have: 1 - 1 T1 S2 S1 T2 T1 T2 |x1 - y1| J () ( ( dlx1 dSy1 ( ( dSx1 dly1 =1 = h1 BC S2 |x1 - y1| + h2 S1 DF |x1 - y1| , (8) where h1 - height LA, h2 - height AM . Differentiating (6) by and taking into account (8), we get: 1 ( ( dlx1 dSy1 ( ( dSx1 dly1 (9) J0 = J (1) = 3 h1 BC S2 |x1 + h2 - y1| S1 DF |x1 1 - y | . 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. 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+1 = (H S + ([A] + [C]) M (Bk )), B0 = 0, k = 1, 2, . . . . This process will be finished when equation residuals become less than the previously set value . For solving linear systems of equations [C]x = y the incomplete Kholetsky expansion method in combination with the conjugate gradient method are used [23]. 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. 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 Conclusion The issues of applying the volume integral equation method to the calcula- tions 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 sub- stitution 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 sug- gested. The results of magnet system simulations based on this approach are in agreement with the calculations by other programs.</p>[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.][A. G. Armstrong, GFUN3D User Guide. RL-76-029/A, CERN (1976).][A. A. Halacsy, Three-dimensional analysis of magnetic fields, in: Proceedings 3rd International Conference on Magnet Technology, Hamburg, 1970, pp. 113-128.][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.][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.][J. Simkin, C. W. Trowbridge, Three dimensional non-linear electro-magnetic field computations using scalar potentials, IEE Proceedings B - Electric Power Applications 127 (6) (1990) 368-374. doi:10.1049/ip-b.1980.0052.][J. A. Stratton, Electromagnetic theory, MCgraw-hill, 1941.][J. D. Jackson, Classical electrodynamics, 2nd Edition, John Wiley & Sons, 1975.][O. C. Zienkiewicz, The finite element method in engineering science, MCgraw-hill, 1971.][J. T. Oden, Finite elements of nonlinear continua, McGraw-Hill, New York, 1971.][W. G. Strang, G. J. Fix, An analysis of the finite element method, Prentice-Hall, 1973.][J. P. Aubin, Approximation of elliptic boundary-value problems, Wiley-Interscience, 1972.][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.][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.][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.][M. Berzins, Mesh quality: a function of geometry, error estimates or both?, Engineering with Computers 15 (1999) 236-247. doi:10.1007/s003660050019.][L. Durbeck, Evaporation: a technique for visualizing mesh quality, in: 8th International Meshing Roundtable, Sandia National Laboratories, South Lake Tahoe, 1999, pp. 259-265.][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][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][M. Abramowitz, I. Stegun, Handbook of mathematical functions with functions, graphs, and mathematical tables, 1964.][I. P. Mysovskikh, Interpolation cubature formulas [Interpolyatsionnyye kubaturnyye formuly], Nauka, Moscow, 1981, in Russian.][P. G. Akishin, The integral equation method in magnetostatic problems: abstract of a PhD thesis, Ph.D. thesis, JINR, Dubna, in Russian (1983).][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.]