The volume integral equation method in magnetostatic problem

Full Text

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 {P¯k , k = 1, . . . , L} is the set of all vertexes in all tetrahedrons {Gi}. Let us introduce the notation B¯k = B¯(P¯k ), H¯k = H¯ (P¯k ), M¯ k = M¯ (P¯k ). We denote fk (x¯) - as a node function, associated with vertex P¯k . The functions fk (x¯) on each tetrahedron are linear functions, equaled to 1 at the vertex P¯k 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¯)B¯k , k=1 L Hˆ (a¯) = fk (a¯)H¯k , 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ˆ = B¯1, B¯2, . . . , B¯L) , T Mˆ (Bˆ) = µ0M¯ (B¯1), µ0M¯ (B¯2), . . . , µ0M¯ (B¯L)) , 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¯| dvx¯dva¯. 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¯)(dS¯x , e¯k )(dS¯a , e¯l ) |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, e¯l) ((x¯ - a¯), dS¯x) |x¯ - a¯| - ∂Gm ∂Gn j ∂fn(a¯) ((a¯ - x¯), dS¯a) (5) - 0.5 ∂al ∂Gm ∂Gn fi(x¯)(dSx, e¯k ) , |a¯ - x¯| where e¯k , e¯l 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 Sˆ1 be a triangle ABtCt, obtained by stretching the triangle S1 in λ times with respect to point A. Similarly Sˆ2 is a triangle ADtF t, obtained by stretching triangle S2 in λ times with respect to point A. Let J (λ) be J (λ) = ( Sˆ1 ( Sˆ2 dSxdSy . |x¯ - y¯| Substituting variables x¯ = λx¯1, y¯ = λy¯1 integral J (λ) can be reduced to: ( ( J (λ) = λ3 dSx 1 dSy1 . (6) S1 S2 |x¯1 - y¯1| Let Tˆ1 be trapezium BtBCCt, Tˆ2 - 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  Tˆ1 S2 Sˆ1 T2 Tˆ1 Tˆ2  |x¯1 - y¯1|  ∂J (λ) ( ( dlx1 dSy1 ( ( dSx1 dly1 ∂λ λ=1 = h1 BC S2 |x¯1 - y¯1| + h2 S1 DF |x¯1 - y¯1| , (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 |x¯1 + h2 - y¯1| S1 DF |x¯1 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]Bˆk+1 = (Hˆ S + ([A] + [C]) Mˆ (Bˆk )), Bˆ0 = ¯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.

About the authors

Pavel G Akishin

Joint Institute for Nuclear Research

Author for correspondence.
Email: akishin@jinr.ru

Doctor of Physical and Mathematical Sciences, Deputy Head of Scientific Department of Computational Physics, Laboratory of Information Technologies

6 Joliot-Curie str, Dubna, Moscow region, 141980, Russian Federation

Andrey A Sapozhnikov

Joint Institute for Nuclear Research

Email: asap@jinr.ru

Junior Researcher of Scientific Department of Computational Physics, Laboratory of Information Technologies

6 Joliot-Curie str, Dubna, Moscow region, 141980, Russian Federation

References