Linear Approximation of Volume Integral Equations for the Problem of Magnetostatics

The volume integral equation method is considered for magnetic systems. New modeling results are reported.


Introduction
Issues concerning the volume integral equation method for the calculations of magnetic systems are considered in the present report.We focus the following discussion on three new points: 1. Short description of the generators of a finite element based 3D grid enabling numerical solution.

Numerical results obtained for quadrupole magnet agree with those obtained by the famous code
TOSCA [1].
Other issues will be discussed in a forthcoming regular paper.

Volume integral equation method
Let B(ā), H(ā), M(ā) be the induction, the intensity and the magnetization of the magnetic field at the point ā.The values B, H, M are connected with the following nonlinear ratios: where µ 0 is the absolute magnetic permeability of the vacuum, µ(x) is the magnetic permeability.
The following integral equation takes place: where HS (ā) is the magnetic field from current winding, G is the area filled by iron.The field HS (ā) can be found according to the Biot-Savart law: e-mail: akishin@jinr.ruwhere J( x) is the current density at x.The difficulty of applying the integral approach is connected with the singularity of the kernel of the integral equations.This is the reason why only a piecewise approximation of unknown parameters within element division area is used in the famous code GFUN3D [2].The alternative for the collocation method is an integration over dividing elements.It allows to use higher order approximation for the unknown variables.The most convenient mathematical approach for constructing such type of approximations is the finite elements method (FEM).
Let us divide the area G into tetrahedrons {G i }.We suppose that the fragmentation G = N i=1 G i satisfies the requirements of FEM.Let us assume { Pk , k = 1, . . ., L} is the set of all vertices in all tetrahedrons {G i }.Let us introduce the notation Hk = H( Pk ), Mk = M( Pk ), Bk = B( Pk ).We denote f k ( x) as a node function, associated with the vertex Pk .The functions { f k ( x)} on each tetrahedron are linear functions, equal to 1 at the vertex Pk and to 0 at any other vertices.Using these notations we define the linear approximations for the vectors B(ā), H(ā), M(ā): We define a discretized formulation of the magnetostatics problem, using the finite element linear approximation within the elements of the division: (

5) Let the matrices [C] and [A] be matrices of [3L × 3L] dimension, consisting of diagonal matrices [C i j ] and non-diagonal matrices [A
The matrices [C i j ] are the following diagonal matrices and the matrices [A i j ] satisfy the following relation for any constant vector M: We use the following notations: Taking into account (1) the system (5) can be written as: Using node functions of higher order similarly to (5) it is possible to formulate a discretization with quadratic, cubic or higher approximation of variables within the element.

Generation of the finite element mesh
To build a discretization for the integral equations ( 5), the region of calculations should be divided into tetrahedrons satisfying the requirements of FEM.Depending on the task, there are different requirements for the mesh elements.In subregions where the solution changes 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 huge number of elements, thus complicating solving the final discretized system of equations.Details of the multidimensional finite element mesh generating are discussed in [3], where the 3D mesh generator 3dfemmesh oriented on modeling electromagnetic fields in large-scale electrophysical machines has been proposed.The algorithm used is based on the representation of the problem domain as a combination of standard macroblocks with initial generation of a two-dimensional mesh on their boundary followed by generation of a threedimensional mesh in each block individually.The program has a graphical interface for the data entry and a visual assessment of the quality of the partition, it calculates a number of criteria for evaluating the quality of the resulting mesh.

Matrix elements calculation
The problem of defining matrix coefficients of the discretized equations can be reduced to calculating sixfold, singular in general, integrals from ( 6) by two tetrahedrons.The simplest way to evaluate these integrals done by a cubature.It might happen that the integrand function shows an isolated singularity or, even worse, the function being integrated is singular at every point of the volume under integration.In these cases the use of a cubature formula is not possible and we need suitable procedures for evaluating the matrix coefficients.The method which allows to decrease the time necessary for the calculation of the coefficients of the discretizated systems is described below.We start from the remark that the matrix coefficients entering (6) can be represented generically as integrals: Taking into account that { f p ( x)} are linear functions and, as a consequence, the function gradient vector is constant inside each tetrahedron, such volume integrals can be reduced to surface integrals: The notation ∂ f m i ( x)/∂x k means that the derivative is calculated on the tetrahedron G m .It is important to note that the region G consists of a union of tetrahedrons.Then the borders {∂G} are triangles.Thus, calculating the expressions (8) reduces a 6D integral to the sum of 4D integrals over two triangles.There are four types of positional relationships of the triangles in space: triangles do not cross, triangles have one mutual vertex, triangles have two mutual vertices, and triangles are congruent.For the first type use of a cubature formula is possible.Meanwhile, in the other cases, the cubature is unreliable due to the singularity of the expressions being integrated.Let us note that the expression under integration in (8) can be represented as a sum of homogeneous functions.The method of integration from [4] allows to reduce the singular integrals of homogeneous functions to superpositions of regular integrals of lower order, the calculation of which can be reliably done by cubature.

Iterative method for solving the nonlinear discretized system
To achieve the requested approximation accuracy in practice it is necessary to split the region G into smaller elements, which leads to huge dimension rise of the nonlinear discretized systems of the equations.The use of methods for handling high order matrices is extremely difficult.So, for solving discretizated systems of equations (7) a simple iterative process is used: This process ends when the equations residuals become less than a previously defined value ε.For solving linear systems of equations [C] x = ȳ , use is made of the incomplete Cholesky expansion method in combination with the conjugate gradient method [5].

Magnet systems modeling
The above mentioned method of volume integral equations was tested for different magnet systems modeling.Here we give a single instance: the computer model of the BOOSTER quadrupole magnet of the NICA project with the help of the generator 3dfemmesh is given Fig. 1a.The yoke field saturation is presented in Fig. 1b.The results shown in Fig. 1c point to the very good agreement of the present method and of the famous code TOSCA [1] based on solving partial differential equations.

Figure 1 .
Figure 1.a. -FEM discretization of BOOSTER quadrupole magnet volume; b. -yoke field saturation; c. -comparison of the outputs of the present method and TOSCA