CC BY
7y^K: 519.8
MATHEMATICAL MODELLING OF THE MAGNETIC SYSTEM BY A.N. TIKHONOV
REGULARIZATION METHOD
Polyakova R. V. Yudin,I.P.
Joint Institute for Nuclear Research, Dubna, Russia
Abstract. A nonlinear magnetostatic inverse problem is investigated for a case when it is needed to create a required magnetic field using conductors which coordinates vary on condition that current value is equal in all the conductors. It is known that the same problems fall under the category of noncorrect problems. Mathematical statement of this type of nonlinear magnetostatic inverse problems is given. The proposed numerical algorithm using the regularization method by A.N.Tikhonov permits to overcome comparatively easily the difficulties connected with badly conditioned equation systems to which usually the magnetostatic inverse problems reduce. This algorithm allows one to calculate the existing winding geometry of superconducting dipole ironless megnet
which ensures a magnetic field homogeneity up to 10 5 within a rectangular aperture.
Key words: magnet systems, inverse problem, Fredholm equation, regularized iterative processes
Introduction
When designing magnetic system it's necessary to solve the inverse problem, that is, via a given magnetic field to define current parameters or, its geometrical characteristics, or all that simultaneously. The definition of the beam density distribution in the magnetic system, in which the geometry is known, is a linear inverse problem for the given field. When the required field must be created with the help of conductors, the value of the current which varies similar to the coordinates of their position providing the current in all the conductors is the same, we come to the solving of the inverse problem. In this paper we consider the construction of a mathematical model of the magnetic system for this kind of the problem and the methods and numerical algorithms for their solution by using the Tikhonov regularization methods. Because a magnetic field is supposed to be given by one of its components (Hx, Hy, Hz) depending on a specific problem, so further H
sample will be used for notation.
Mathematical Model of the Magnetic System
Let in a region U with the help of the sources of current distributed in the region S a field H should be created with the given characteristics (for example, the whole homogeneous field in the region U .) It is known [1], that the field in any point z of set U is defined by the expression
H(z) = £ J(s)G(z, s)ds, z eU, s e S
where J(s) is a distributed density function of the current in the system, G(z, s) is a Green function ,
that analytically depends both on the geometry of the source of the magnetic system and on the point z eU. The inverse problem, namely, a definition over the given density of distribution of current in the magnetic system with the known geometry is a linear inverse problem (model 1). Then the mathematical problem reduces to the solution of the Fredholm linear integral equation of the first order with unknown function J (s).
If the composition of the magnetic field includes not only variant density of current and arrangement source of current, then we must solve the nonlinear inverse problem (model 2) with unknown J(s) and s e S .
Method for Solving the Inverse Problem (Model 1)
It is known that the problem of solving the first order Fredholm integral equation (1) belongs to the incorrectly defined class of problems, because large changes in the solution of J(s) can correspond to small changes in the input data H(z). To obtain a stable solution of incorrectly defined problems, A.N. Tikhonov developed regularization algorithms [1-4]. Here we use the second order regularization method for solving the problem. To do this, we construct a smooth parametric functional
Fa[J (s), H (z)] = 0[J (s), H (z)] + a Q[ J (s)], m
where
0[J(s),H(z)] = £ [H(z)-¡s J(s)G(z,s)dsf dz
is the quadratic deviation of the operator A[z, J(s)] = I J(s)G(z, s)ds of the function H(z),
J S
Ц J (s)] = £ J 2(s)ds
is the regularization functional, or stabilizer, and a is the numerical parameter of the regularization (a > 0).
Theorem1. For any function H(z) e L2 and any a > 0 there exists one and only one 2(n +1) differential function J a (s) that realizes the minimum of the smooth functional Fa[ J (s), H (z)] of the form (2).
Theorem 2. If H(z) = A[z, J(s)], J(s) e C^n+1), then for any S> 0 and auxiliary values 0 < yx <Y there exists S(s, Y\,Yi, J) , so that if
i • и нsi?) - H(z) ii2 < 8, Hs(z)e l2 ;where
2. a = a(S) has the orderS2;
3. Y <Yi,
then Jas that realize the minimum F" [J(s), f/ fV(z)] are
□ seS,i = 0,l,-,n 5<5Q{£,yvy2,J).
It follows that there exists some function Ja (s) that realizes the minimum of the functional
F^[J(s),H(z)] in the form (2), which weakly reduces to the solution of the equation (1) J(s) . When
applying the regularization method, the choice of the parameter a is one of the main problems. The point is that
the discrepancy principle is not always fulfilled for the obtained smooth solution, i.e. the inequality holds
□ (5)
where 5 is the accuracy of the approximate input data ff ri, H is the precise value of the input data. In
practice, when solving correct and incorrect problems, it is necessary to find a solution that satisfies the required accuracy. V.A. Morozov [2] proposed the discrepancy principle as the main criterion for choosing the regularization parameter.
Discrepancy Principle. Set any 0 < S <S0 and any 0 < h < h0 with the condition
X(h, S,J) = (U AhJ -AJ q, +5) x (1 + J3(S, h))m <UAhJ-Hs 0,
where f(S, h) is the positive function, such that limS f(S, h) = 0, Ja is the solution that gives the
minimum of the functional Fa (J, H). In this case there exists at least one value of the regularization parameter a = a(S, h) > 0 so that
Psh(a(S,h)) = x\h,s,J) and lim Jash= Ja where Jash= Jashsh).
S,h^0
Numerical Algorithm for Solving the Problem of Model 1
In the expression (2), if presenting integrals in the form of sums, we obtain
N M M
Fa=£ [Hj (z}) - £ J1 (s, )Kj (z}, ^ )]2 + a£ J (s, )As,
j=x i=x i=x (6)
where N is the number of points from the set U , M is the number of points from the set S , M < N, K. = [ G(z , 5)ds. Suppose As — As = const, Az. = Az = const.
' J As j
The condition of the minimum of the functional Fa :
dFa _Q dFa _Q dFa _Q
a^i ' dJi ' (7)
Taking into account (7), we obtain
p pa N N M
= "£ HKAz + £ £ JKyK, Az + aJ As = 0, I = i M.
J j=i j=i i=i
(8)
Thus, we have a system M of linear algebraic equations with the unknowns
M of the form J :
£ Ji £ Kj + aJA = £ H]KyAz, I = 1+ M.
'=1 j=i j=i (9)
Supposing aAs = a'Az, we obtain
M N N
X 4 X KK +a'Jr = X HJK'J' / = 1 -M■
1=1 j =1 j =1 (10)
Obviously, a' retains the meaning of the arbitrary coefficient a , therefore, the system of equations for J is finally written in the form:
M N N
X J X KjKj +aJ/ =X HjKj, / = 1 - M ■
1=1 j=1 j=1 (11)
If the magnetic system is a discrete set of coils, the field H(z) in any point z e U is defined in the following way:
m r
H(z) = X J1 f G(z,5>fc, i=1 1 (12)
where M is the number of coils, J is current density in the i-th coil, A,- is the cross section of the i-th
coil.
Having solved the system of equations (11), we obtain a discrete set J, / = 1 — M, which is a solution of the problem (1). Similar to that, we define the distribution of current density in the magnetic system for creating the field H(z;), j = 1 — N, z; e U.
Particular Case of the Mathematical Model 2
Suppose that in some region S , using the arrangement M of conductors with the same current /0 , it is necessary to create the field H. In the system H(z) , z e U we have
H(z) = Io £ Gis,, z),
'=1 (13)
where G(s,z) is the Green's function for the i-th conductor. It is required to define both the current I0 and the coordinates s of the conductors, which would best provide the given field H(z), z e U. The function G(s, z) is usually nonlinear with respect to the coordinates of the conductors s ; therefore, the
problem under consideration is nonlinear inverse.
Additional difficulties in solving the inverse problem are caused by the restriction on the parameters [3]. However, in some particular cases one can effectively find a solution that satisfies the conditions of the problem. Let us consider this case. Let the parameter of the conductors arrangement in the region S be only one coordinate, for example x; the admissible region of the conductors arrangement along the axis x is known,
xi Xj X2.
In this case, the equation (13) has the form:
M
H(z) = I0 £ G(X, z), z e U, X < X < X2.
'= (14)
H(z),
It is necessary to determine I , x to create the field in the magnetic system.
z eU
The problem (14) is a nonlinear inverse problem. Numerical Algorithm for Solving the Problem in Model 2
To avoid the difficulties associated with solving such systems of equations [4], problem (14) is solved in two stages. At the first stage, current density in the turns is distributed continuously within the given region. The
rx2
equation (13) has the form H (z) = I J( x)G( X, z)dx.
Jx!
This problem and the algorithm for its solution are considered in sections 2, 3. When choosing the solution Ja (x), (a is the regularization parameter), we take into account the following conditions of the problem:
1. The accuracy of the calculation of Н (z) cannot be worse than the required accuracy of the magnetic field in the magnetic system being created;
2. For the entire admissible interval [ x, x2 ], the function Ja (x) must keep the sign;
3. I Ja (x) |< Jdon is admissible current density.
Suppose that there exists some continuous solution Ja (x), that satisfies all three conditions. At the second stage, we divide the interval [x, x2 ] into М subintervals [x[, x'2 ], i = 1 ^ М. In this case
H(z) = £ Jj2 Ja (x)G(x, z)dx, z e U.
'= x (15)
For each subinterval [x\, x\ ], the conditions of the mean value theorem are satisfied (this is how the function Ja (x) is chosen), therefore
H(z,) = £ G(xj, z,) Jx Ja (x)dx, j = 1 * N,
'= xi (16)
where N is the number of points in the region U , in which the field Н, xj , i.e. the point in the i-th interval, is analyzed. The limits x', x\ are chosen so that
J4 Ja (x)dx = Jx^ Ja (x)dx = I0, т.е.
•M + (17)
fX2 J a( x)dx
t _ Jx1
1 A '
0 M
Ml
H ( z7 ) = /0 £ G( xj, Zj ) (18)
Obviously, for different z. ther exists its point xj , but based on the mean value theorem, it is always in
the interval [x\, x\ ]. It means that the coordinate xi is also in the i-th interval and is defined from the condition of the minimum of the functional
P(X ) = £[ JJ2 J" (x)G(Xzj)dx " /oG(xi, zj)]
7 a
j-i x1
- £ [J X2 Ja ( X)G( X, Z; )dx - /oG( Xi, Z; )] X - 0
CXi j-1 X1 °Xi XÎ < X< X^ (19)
In the Cartesian coordinate system s — (X, y }, z — {x, y}
G(s , z) - y-y^ ln(X - Xi + a)2 + (y - yi + b)2 +
2 (x - x, - a)2 + (y - y, + b)2
y - - bin( x - x-- a)2 +( y - y - b)2 +
2 (x - x, + a)2 + (y - yt - b)2
f \ f \ x - x. + a x - x. + a . . x - x. - a x - x. - as
+ ( x - x - a) x arctg-i—- - arctg -
(x - x + a) x
arctg----arctg ■ / _____^ _____^
y - y, - b y - y, + b ) y y - y, + b y - y, - b
(20)
where a is the half-size of the tire along x, b is the half-size of the tire along y , G(si,z) is the Green's
function for the rectangular tire in the Cartesian coordinate system. Thus, the solution of the problem reduces to solving M sequential nonlinear equations of the form (19) with one unknown, moreover, limits of the existence of the solution are known. It is noteworthy that we have analyzed the algorithms for creating a magnetic system using infinite thin conductors. It is easy to show that for the finite size of the conductor, the algorithms are fully preserved, but in this case the Green's function is under the sign of integration by the cross section of the conductor. The case when the admissible geometric region of arranging conductors is defined not by a line, but by a plane, also does not involve particular difficulties and can be described by a similar algorithm. An example of the mathematical modeling of a real superconducting (SC) dipole magnet, comprising rectangular current excitation windings, is described below.
Modeling a real magnetic system
This paper considers the solution of the inverse magnetostatics problem for an iron-free superconducting (SC) dipole magnet consisting of rectangular excitation current windings. The required magnetic field is created by selecting the position of the windings. In this case, the current windings must be powered from only one power supply and the current in all conductors (turns) forming the windings is the same.
Specific conditions, arising from the possibilities of manufacturing and using a magnet at the accelerator of LHE JINR, are imposed on the geometry of the magnetic system (Fig. 1), the magnitude of the current in the SC turn, as well as on the nature and magnitude of the magnetic field. Mathematically, the problem is reduced to the class of ill-posed problems and solved by the regularization method according to A.N. Tikhonov (described in detail above).
2
y
Fig. 1 Fig. 2
Li
Let it be necessary to create a uniform magnetic field H(x, y) with a single non-zero component Hy (x, y) in the aperture of the magnet (Fig. 1). We assume that the longitudinal dimensions of the magnet are much larger than the transverse ones and the calculations do not depend on the choice of the section. Let us write down the well-known expression (Bio-Savart law) to calculate Hy (x, y):
un jCXfMx -
Hy (X, y) = TT
S fx - X)2 + (y - y) *
I0 G (/? x, y),
(21)
~ 3/2 2 *
where (x, y)0SH, S = (Ox,S7, I0 is the value of the current in the section S7, J(x,ft) is the current density, G(/,x, y) is the Green function for the rectangular turn given above. The current winding
(region Si ) consists of m — N + N turns, where N is the number of winding turns along the axis x, N is the number of winding turns along the axis y . The transverse dimensions of the turn are 2a = 0.40 cm in width and 2b = 0.15 cm in height. As mentioned above, the current density in the turn is constant, i.e. J(x, y)= Jo . In this case, for the region Si , consisting of m identical turns with the current density Jo ,
Eq. (21) can be presented in the form
m x + ayt + b
Jo tfe T T
'=1 x - ay - b
(x - i
- x) + (y - y) |
x
,372 = Hy (x y)=
OA2 *
(22)
where (xt, yt) are the coordinates of the center of the i-turn; Ho = const in the region Sh . Thus, to
solve the problem, it is required to select Jo and the centers (xf, yf) so that Hy (x, y)= Ho in the region
SH . Such a problem belongs to the class of nonlinear inverse problems, which, as is known, for sufficiently large m, belong to the class of ill-posed problems.
The solution of our ill-posed problem is reduced to solving an ill-conditioned system of 2m +1 nonlinear equations [1, 4]. To avoid the difficulties associated with solving such systems of equations, problem (22) is solved in two stages. At the first stage, let us consider that the turns are fixed and uniformly distributed over the
winding, and the desired function is the current density J(x,y). Thus, problem (21) becomes linear. Let us select two parts from its solution (function J(x,y)): J(x£)= Jp= const), continuous along the axis
and Jy (y)= J(X^ const, y), continuous along the axis y. At the second stage, the obtained
x,
continuous distributions Jx (x) and Jy (y) are modeled using piecewise constant functions, i.e. "blocks"
under the condition of the constant current in one turn. ("Blocks" are formed from coils of standard sizes: 0.4 cm x 0.15 sm). It is assumed that the "block" functions composed in this way give a solution to problem (22) with all the restrictions imposed.
Stage I. It is noteworthy that equation (21) at this stage takes the form of a linear Fredholm integral
equation of the first kind. The problem of determining the distribution J(X,y) from the known Hy (x, y) is
unstable and belongs to the class of ill-posed problems. This is illustrated in Fig. 2, which shows the solution of this problem obtained by the authors by one of the usual methods without introducing the regularization. Later, the problem is solved by the Tikhonov regularization method. The solution is reduced to defining the function
J (X,y) minimizing the smoothing functional Ma :
Ma (J, H)= 0 Jy), Hy (x, y)*+ a
(23)
where
0
$), H y ( x, y)P= TT [H y ( x, y) - TT J (X?°y)K (x, y, X>y)dx&y]2 dxdy
S
H
S
(24)
is the quadratic deviation of the operator tt J ( XyyK (x, y, x,y)axSy from the given function
SI
Hy (x, y) = H0 in the region SH, K (x, y, X0.
x -
, 3/2
0A2 ^
in the region Sj
(x - X) +(y - y) I
o TT j2 (Y,fyd xd y is the stabilizing functional, a is the regularization parameter, which is
selected from the condition
tt [Hy (x, y) - tt Ja (XX,(x, y, X,y)dXay]2 dxdy J d (25) SH SI
where Ja is the element at which, for the given a , Ma reaches a minimum for the given 5 . To perform numerical calculations, let us discretize. We write the previous expressions in the following
form:
Ma=Y
j=i
m
Hj ( xj, yj Jr (, y. ) Kj ( xj, yj, xi, yi ) i=l
Jt 2 ( x., yt )AxAy
i=1
(26)
where (xi, yi) e Si , ^x^, yj ^ G SH , n is the number of points in the region Sh , m is the number of points in the region SI ,
xi +ayi +a ( xj - x ) dxdy y - yj + ô ( x. - xj + a ) 2+( yi - yj + b)
2
Ki, = i I
xj -ayi -a
(xj -x)2+(yi-y)2
2/3
(27)
(xi -xj -a)2+(yi -yj -b)
The elements ,/ . i — 1, 2,..., m, minimizing Ma , are found as a solution to the system of equations
dM
a
dJ
dMa
■ = 0,
dJL
= 0,
m
(28)
where
CM
a
cJ,
n n HI / \
: -2 S H K • AxAy + 2 S £ J ( X, y ) KuKaAx Ay + 2a- J -Ax - Ay j=l J j j=ii=i iv 1 l! j j 1
n m
7 J=1 - " J=u=-
Assuming a ■ AxAy = a-A x ■ A y , we get
(29)
S J £ KaKj- J =S HJ Hj
i=1 j=1 j=1 (30)
The parameter a is chosen so that the quadratic deviation
'=£
j=i
H -
m t \
S Jia(x, y ) Ka (XJ, yj, xi, yi )
AxAy
i=i
was minimal, such that F <S< 10"6. Here Ja is solution (30) for the given a . In (30), Ax = Ay = 0,25 sm, A x = A y = 0,055 sm and a = a ■ 21. Table 1 shows the dependence F on a in the interval a = 10"8 ^10"16 . a = 10"11 (a = 21-10"11 j is selected as the best. Figs. 3a and 3b show the dependences Jx on x (x = 0 + 3.6 sm; y = 4 sm) and Jy on y (x = 4 sm; y = 0 + 3.6 sm) obtained as a
solution to system (28) at a = 10 11.
Table 1.
Dependence of the functional F, the square of the residual on the regularization parameter a
m
n
n
2
a 10-8 10-9 10-10 10-11 10-12
F 271-10-10 79 -10-10 5,5-10-10 0,14-10-10 2,6 -10-10
a 10-13 10-14 10-15 10-16
F 55 -10-1 0 876 -10-10 2,8-10-10 0,354
Stage II. At the second stage of solving the problem, the continuous distributions Jx and J are modeled
by piecewise constant functions, namely, "blocks". In this case, it is required to fulfill the following conditions arising from the practical possibilities of manufacturing the magnets under study:
1) the current density in all blocks is constant and the same: = J{). i — 1, 2,..., k, k is the number of blocks;
2) each block must correspond to an integer number of turns, i.e. the length of each block is a multiple of the height of one turn 2b (in our case 2b = 0.15 cm);
3) in the transition from the continuous distributions Jx and Jy to the piecewise constant distribution, it is necessary to lose as few turns as possible.
a b
Fig. 3 shows the dependences " x on x (x
J,
0 ■ 3.6 cm;
■ y
4 cm),
J,
on
y (X = 4
cm
• y = n +
0 ■ 3.6 cm), obtained as a solution to system (30) at
a = 10
-it
In the first approximation, the blocks are placed so that the gaps between them correspond to the minima on the curves Jx, Jy. The block areas St are determined by the numerical integration Jx, Jy using the
trapezoidal method from the left minimum of the resulting continuous distributions Jx, Jto the right
minimum. The current densities in the i-block Ji are determined as Ji = si / , where lt is the distance
between the adjacent minima. Then, from the condition J{) — maxJ^ i —1,2, the current density,
which is the same in all blocks, is chosen. In this case, the block sizes do not go beyond li . Further, since there
must be an integer number of turns in each block, the block lengths l; are determined as l- = 0,15 x N., where
N is the integer part of the ratio s^ /10,15 x Jq j. In this first approximation, placing the blocks between all
the minima results in a large loss of turns. The smallest losses of turns occur when three blocks are placed on each of the continuous distributions Jx and Jy (Fig. 3).
The dimensions of the magnet blocks with different configurations are shown in Table 2. After such a preliminary placement of the blocks, the problem is further solved using the FUMILI program, i.e. the minimization of the quadratic functional /5/. The coordinates of the block centers are taken as variable parameters. The results of such minimization are given in Table 3. Fig. 4 demonstrates one of the possible
—4 .
configurations of a real magnet, in which the field uniformity AB / B < 1-10 'is 80% of the aperture size.
In this paper, the solution of a nonlinear inverse magnetostatics problem is studied in the case when the required magnetic field must be created using conductors whose coordinates vary, provided that the current in all conductors is the same. It is known that such problems belong to the class of ill-posed problems. A mathematical formulation of this class of nonlinear inverse magnetostatics problems is given. The proposed numerical algorithm, using the regularization method according to A.N. Tikhonov, makes it relatively easy to overcome the difficulties related to solving ill-conditioned systems of equations, to which inverse magnetostatics problems are usually reduced, and to calculate the geometry of the excitation windings of an iron-free superconducting dipole
magnet, which ensures the uniformity of the magnetic field up to 10 5 inside the rectangular aperture.
Fig. 4a, b: a - SC dipole configuration with a flat winding; b - own of the possible real configurations of the SC dipole.
Conclusion
1. A method of solving nonlinear inverse problems, which are necessary to describe the mathematical model of magnetic systems of some class, is considered.
2. The developed numerical algorithm, based on the method of regularization of solving incorrect problems with restrictions on the required parameters, leads the nonlinear problem of type (14) to solving M sequential nonlinear equations with one unknown. It allows avoiding the difficulties related to solving a system of nonlinear equations, the solution of which is usually given by inverse problems.
3. To implement the proposed method on a computer, a numerical algorithm is developed, and a Fortran package is written.
4. Using the given package, some practical problems [5-7] are solved, an example of one of which is given.
5. The results of numerical modeling of some real systems were used in the design and construction of the superconducting accelerator of JINR LHE.
The work is dedicated to the memory of Professor E.P. Zhidkov and Professor I.A. Shelaev.
References
[1] V.Y. Arsenin, A.N. Tikhonov. Methods for Solving Ill-Posed Problems. Moscow: Nauka, 1979, in Russian.
[2] V.A. Morozov. Numerical Methods and Programming. Issue 8, MSU Publishing House, 1967, pp. 6395, in Russian.
[5] B.T. Polyak. Iteration Methods for Solving Some Ill-Posed Variation Problems. Numerical Methods and Programming. Issue 12, MSU Publishing House, 1969, pp. 38-52, in Russian.[4]
[5] E.P. Zhidkov, R.V.Polyakova, I.P. Yudin. A.N. Tikhonov Regularization in a
Magnitostatic Problem. MCE: Proceedings of the 13th International Conference, 2006, pp. 11-177, in Russian.
[6] E.P. Zhidkov, L.L. Zinovieva, R.V. Polyakova, I.A. Shelaev. Solution of a certain class of nonlinear inverse problems of magnetostatics. 11-10845, Dubna, JINR, 1977, in Russian.
[7] E.P. Zhidkov, I.V. Kuts, R.V. Polyakova, I.A. Shelaev, I. P. Yudin. Solution of a nonlinear inverse problem of magnetostatics by the regularization method. P11-88-335, Dubna, JINR, 1988, in Russian.
