- Email: [email protected]
Effective magnetic field computation in Tokamaks in presence of magnetic materials
Fusion Engineering and Design 96–97 (2015) 641–644
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Effective magnetic field computation in Tokamaks in presence of magnetic materials Andrea G. Chiariello a , Alessandro Formisano a , Raffaele Fresa b,c , Francesco Ledda a,∗ , Raffaele Martone a , Francesco Pizzo a a
Dip. di Ing. Industriale e dell’Informazione, Seconda Università di Napoli, Via Roma 29, I-81031 Napoli, Italy Scuola di Ingegneria, Università della Basilicata, Potenza, Italy c Consorzio EURATOM/ENEA/CREATE, Italy b
a r t i c l e
i n f o
Article history: Received 3 October 2014 Received in revised form 15 June 2015 Accepted 16 June 2015 Available online 10 July 2015 Keywords: Magnetic field evaluation Ferromagnetic materials Fusion Technology Neutral Beam Injector
a b s t r a c t Magnetic materials play an important role in the magnetic fields distribution inside Tokamaks, and their contributions must be carefully computed. In some applications, when 3D geometries are involved and high accuracy must be achieved, finite elements may reveal too computationally demanding and different approaches must be considered. In this paper, two possible solutions are presented. The first one is based on the discretization of magnetic parts into rectangular prisms with constant magnetization density; analytical formulas are then used to compute the field contribution due to each brick. The second one is based on the representation of the effect of magnetic materials as a set of dipoles, by using again analytical formulas to evaluate the contribution to the magnetic field. Iterative procedures are applied to evaluate the magnetization density in prisms. The efficiency and accuracy of both models are validated against commercial codes on a test case, and then applied to the computation of magnetic field in presence of Neutral Beam Injector in an ITER-like geometry. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Accurate knowledge of magnetic fields in Tokamaks represents a key priority in the design of present devices, where size and power levels make the accurate control of fields mandatory to achieve desired performance. Discrepancies with respect to the design field map are originated not only by magnets fabrication and installation tolerances, by the joints and by the busbars, but also by the presence of ferromagnetic elements. Their combined effect can be very hard to be estimated using standard computational techniques, e.g. Finite Elements Method (FEM), since high accuracy is required, but no symmetry can be exploited to reduce the model complexity. A number of approaches were proposed to cope with deformed magnets, joints and busbars [1–3], but none of them, at authors’ best knowledge, is specifically tailored to cope with magnetic materials. A further example is the effect of ferromagnetic parts on magnetic diagnostics signals for plasma identification, where again the contribution from iron, eventually lacking any symmetry, must be accurately computed. Magnetic materials are present in Tokamaks, for example in the test blanket modules (TBMs) or in the TF ripple
∗ Corresponding author. E-mail address: [email protected] (F. Ledda). http://dx.doi.org/10.1016/j.fusengdes.2015.06.078 0920-3796/© 2015 Elsevier B.V. All rights reserved.
correction system, or in the passive shielding of Neutral Beam Injectors (NBI). Accurate field computation in large 3D structures in presence of magnetic materials has been widely discussed in the past. Different methods have been adopted, including integral methods [4,5], “smart” formulations implemented in FEM models to save memory and computational time (sparsification methods [6], Fast Multipole Methods [7]), or FEM-BEM formulations [8]. The development of FEM and BEM approaches did somehow put aside methods based on equivalent sources (see, e.g. [7,9,10]). Such methods may represent again a viable alternative thanks to the involved semianalytical formulas, relatively straightforward and well suited for implementation on architectures based on Graphic Processing Units [11]. In this paper, two approaches are described to reduce the computational effort in presence of fully 3D ferromagnetic parts. Both take advantage from discretization into elements possessing analytical expressions relating the field to the magnetic equivalent sources. Namely, a discretization into rectangular prisms with constant magnetization density and a representation using sets of magnetic dipoles are considered. The effectiveness of the methods is first tested against 3D FEM using a simple geometry, and then applied to the evaluation of the effect of the Neutral Beam Injector (NBI) passive shields in an ITER-like device on the magnetic
642
A.G. Chiariello et al. / Fusion Engineering and Design 96–97 (2015) 641–644
field at the first wall, where the diagnostics are placed. The paper is organized as follows. In Section 2 the basic mathematical models are discussed, in Section 3 a short description of the NBI passive shields is presented; Sections 4 and 5 show the results achieved in a couple of examples; finally, in Section 6 some conclusions are derived. 2. Mathematical model Although other cases could be treated in a similar way with minor modifications, for the sake of exposition in the following it is assumed that the time evolution of the magnetic field is much slower than the time constants of the passive structures. This is certainly valid during flat-top of the plasma discharges, but may be critical during disruptions or ramp-up and ramp-down. Under these assumptions, a non-dynamical model can be used to describe the flux density profile: div B = 0
(1)
rot H = J
(2)
If the field is required also when dynamical effects are relevant, then the model must be modified to include the full set of Maxwell equations. Work is in progress from authors to extend the computational codes to cope also with fast parts of the discharge. The magnetic constitutive relationship is assumed in the form B = 0 (H+M(H)), where M is the H-controlled non-linear function providing magnetization vector. It is assumed that the iron constitutive relationship is isotropic, and does not show hysteretic behavior. The whole analysis domain ˝ can be divided in four different regions: ˝plsm
J = Jpl M=0
Plasma region : the current is derived from
equilibrium equations
˝coil
J = Jext M=0
Active coils region : the current is imposed by
power supply system
˝iron
J=0 M = M(H)
˝iron
M=0
(3b)
Iron regions : the magnetization is derived from
the material characteristic
J=0
(3a)
Vacuum region
(3c)
(3d)
Note that, thanks to non-dynamicity assumption, passive structures can be neglected, since no eddy currents are present. Of course, a suitable set of boundary conditions have to be included on the external boundary of ˝0 . The computation of magnetostatic field in presence of magnetized bodies can be decomposed into three logical steps: 1. Evaluation of magnetic field Hsrc from impressed currents in free space; 2. Evaluation of magnetization distribution M in magnetic materials; 3. Evaluation of field contribution Hiron due to magnetization, in the free space.
Fig. 1. The elementary “brick” used for the iron parts discretization.
The first step can be done using Biot–Savart law [12], or any other method suited for impressed currents, and will not be discussed any more. The second step can be pursued by decomposing M on a suited representation basis, and solving the non-linear integral equation with a suitable tool, such as, e.g. the Picard approach. The best-known methods in this class represent the material with “equivalent” magnetic sources (magnetic dipoles or fictitious currents/charges [12]), and include the source parameters among the unknowns. This approach will be used here, discretizing the magnetic body into a suitable number of elementary subdomains, whose shape is chosen based on the availability of simple and accurate formulas to compute magnetization and magnetic field. Two possible analysis methodologies will be developed in the following: 2.1. Rectangular prisms (bricks) In the first approach, magnetic parts will be decomposed into rectangular prisms, “bricks”, and a uniform magnetization density is assumed in each brick as in Refs. [7,9]. The magnetic scalar potential due to an assigned magnetization map M(r ) can be written as: iron (r)
=−
1 4
∇ · M(r ) ˝iron
1 d˝iron |r − r |
(4)
where ˝iron is the region of space filled with magnetic material, r is the “field point” coordinate vector, and r is the “source point”. The contribution of the iron on the magnetic field is then computed as Hiron = iron . If M is constant over a “brick” aligned along one of the coordinate axes, (4) can be computed analytically in closed form, since the divergence operator is vanishing everywhere in the inner domain ˝iron of the iron domain ˝iron and just the contribution due to the discontinuity of the normal component of M must be considered. For the sake of exemplification, the magnetic potential z due ˆ as in the brick depicted in to a z-aligned magnetization M = Mz Z Fig. 1 is 1 (−1)k−1 4 k=1.8
z − z 2 k 2(z − zk ) atan + atan + 1 y − yk
z (r)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨
=
(x − xk )asinh
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬
y − yk
(x − xk ) + (z − zk ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ x − xk ⎪ ⎪ ⎩ (y − yk )asinh 2 2
2
(y − yk ) + (y − yk )
+
(5)
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ 2
where (xk , yk , zk ) represent the coordinates of brick vertices in Fig. 1, and
1 = (x − xk )2 + (z − zk )2 + (x − xk ) 2 = (y − yk )(z − zk )
(x − xk )2 + (y − yk )2 + (z − zk )2
A.G. Chiariello et al. / Fusion Engineering and Design 96–97 (2015) 641–644
643
Eq. (5) was derived starting from (4) by using the symbolic solver of Matlab© R14 to achieve best performance when using HPC algorithms and it has been compared to similar expression from Ref. [9] to assess correctness of resulting values. The magnetic field Hiron of a uniform magnetization along z in the brick will then be given by: Hiron (r) = Mz (hxz xˆ + hyz yˆ + hzz zˆ ) where hkz = ∂ z /∂k (k = x,y,z). Similar results of the scalar potential can be obtained also when the magnetization M is x or y-aligned. Since computation in this stage assume that magnetization map is assigned, it is possible to decompose any distribution of magnetization into the three components M = Mx xˆ + My yˆ + Mz zˆ , and obtain:
⎡
⎤
hxx
hyx
hzx
Hiron (r) = hM(r ) h = ⎣ hxy
hyy
hzy ⎦
hzx
hzy
hzz
⎢
⎥
(6)
where the tensor h links the effect of the uniform 3D magnetization of the brick to the magnetic field. It should be noticed that (6) could be easily generalized for a brick not aligned with the coordinates as in Ref. [13]. When more than one brick or field points are considered, each pair of them will be treated as in (6) and a suitable matrix of tensors h should be considered including three rows for each field point and three columns for each brick. Limiting in a Cartesian coordinate system, the tensors matrix can be suitably treated as 3 × 3 partitioned matrix. An iterative scheme based on collocation method is used to evaluate the magnetization of the bricks. First (6) is used at each iteration to obtain the magnetic field in the center of each brick from the magnetization map at previous iteration. Using the non-linear constitutive relationship M = M(H), a new estimate of M (constant over each brick) is then obtained, until convergence in Cauchy sense is achieved: Hnk+1 = Hnsrc +
n,m
m=1:Nbricks
h
Mm k
Mnk+1 = ˇkn M(Hnk+1 ) + (1 − ˇkn )Mnk
(7) (8)
n=1,2,...Nbricks
where Nbricks is the total number of bricks, k is the Picard iteration index, hn,m is the 3 × 3 submatrix of h characterizing the contribution of the m-th brick on the n-th one whose elements are computed according to (6), and ˇkn = 1/(1 − nk min ) is an underrelaxation factor, required to achieve convergence [7], k is the k-th iteration estimation of the relative permeability of the n-th brick, and min is the smallest eigenvalue of h. A Newton schemes may speed up the convergence of the nonlinear iterative scheme but they would require a matrix inversions at each iteration step. Using (7) and (8), very large magnetic structures may be treated by assembling “on the fly” the required matrix elements, saving memory space and possibly relying on GPU speed to limit the computing burden. 2.2. Equivalent dipoles The bricks model allows self-consistently computing magnetization map M and field due to magnetic parts Hiron . Anyway, once the magnetization is known, either using bricks or using FEM, the effect of magnetic materials can be taken into account in remote regions (e.g. in the plasma chamber) by using faster formulas for the “long distance” effect. In fact, the magnetic material can be partitioned into smaller volumes, e.g. spheres, and each sphere
Fig. 2. A schematic view of the NBI and vessel considered geometry. Vessel sections are used to define test points.
of radius rsp can be associated to a magnetic dipole of assigned momentum: m=
4 3 r M 3 sp
(9)
The value at the generic point r of the magnetic field associated to a magnetic dipole is: Hiron (r) =
−m + 3(m · l rr )lrr
(10)
4|r − r |3
where r is the point where the dipole is collocated and ˆirr is the versor oriented along the direction r − r . 3. Neutral beam magnetic field reduction system Two Heating (HNBI) and one Diagnostic (DNBI) Neutral Beam Injectors are foreseen in ITER and a third heating injector could be installed in a successive operation phase. The residual magnetic field inside the NBI region is required to be low enough to avoid deflection of the ion beam and to achieve the design performance of these components. Both passive and active shielding systems have been foreseen, to reduce the high stray field inside the NBI due to the plasma and to the Poloidal Field Coils to an acceptable level. The effect of NBI passive shielding anyway is not negligible also outside the injector: in this paper, the interest is on the field contribution in the vessel surface, where magnetic diagnostics are placed (Fig. 2). 4. Test case In order to assess the accuracy and the computational effectiveness of the two models, a simple test case is considered: a hollow rectangular prism of a magnetic material, with similar dimensions of the external NBI iron shield is immersed into a uniform Hsrc field, vertically directed along z-axis. The resultant magnetization map M has been compared among the results given by the two models and those provided by a commercial FEM. Details about geometry, magnetic field and discretization levels are reported in Table 1. The magnetic field H along the symmetry axis z has been compared Table 1 Details of the test case. Geometry of hollow prism External dimensions Wall thickness Magn. permeability
x = 13.8 m; y = 5.11 m; z = 5.03 m 0.15 m rel = 1000 (constant)
External field
Hsrc = 1000 AT/m along z direction
Details of models Bricks Dipoles FEM
363 on long sides, 121 on short, 1694 total Same as bricks 3840662 DoF
644
A.G. Chiariello et al. / Fusion Engineering and Design 96–97 (2015) 641–644
Fig. 3. The comparison between the FEM, dipoles and bricks.
ferromagnetic shield. The shield is supposed to be constructed using SS400 steel. The source field distribution Hsrc is due to Toroidal field coils, Central solenoids, poloidal field coils and plasma. Geometries and currents are similar to one of the ITER reference scenario, namely Start of Flat Top equilibrium in the 15MA inductive scenario [14]. To assess the effectiveness of the approach, the logarithmic variation of field log10 (||Biron ||/||Bsrc ||) on the vessel is computed. In Fig. 4 the norm of the magnetization M on the NBI shield walls is reported; the effect of the shield on the vessel first-wall is reported in Fig. 5 both as average and as maximum value in each poloidal cross section along the toroidal direction. 6. Conclusions Elementary sources with known analytical expression can provide useful tools, alternative to FEM codes, for the evaluation of the effect of ferromagnetic materials. In this paper, two methods, based on analytical expressions for the magnetic field due to ferromagnetic parts immersed in a field due to external sources have been proposed. Advantages in terms of accuracy and computational times have been assessed using NBI passive shields as examples of ferromagnetic parts. Work is in progress by authors to extend the model to cope with “fast” current dynamics, and possibly with hysteretic materials. Acknowledgments
Fig. 4. NBI magnetization map.
This work was partly supported by the EU under the Contract of Association between EURATOM and ENEA/CREATE (F4E-2008-OPE06-04-03(ES-AC) Task Order n .3 Analysis of the error field induced in the magnet system of ITER), and partly by the Italian MIUR under PRIN grant 2010SPS9B3. References
Fig. 5. Effect of magnetization on the vessel.
using a FEM model, bricks model and dipoles model, with M provided by bricks. The results are reported in Fig. 3 and they allow to state that the three different models are all effective in computing the required field at long distance from the source while the dipoles model can provide inaccurate result at the proximity of the source. On the other hand, the dipole model on serial C code is about 82 times faster than the brick model to compute the magnetic field due to an assigned magnetization map. 5. Analysis of NBI impact In order to test the methods in a realistic case, the effect of a single NBI passive shield was computed using the brick model. Since the interest in this work is to assess the effectiveness of the proposed methods, the geometric NBI details are not considered and a simple hollow rectangular prism is considered to model the NBI
[1] J. Knaster, et al., ITER non-axisymmetric error fields induced by its magnet system, Fusion Eng. Des. 86 (6) (2011) 1053–1056. [2] V. Amoskov, et al., Assessment of error field from solitary ferromagnetic elements located outside of ITER tokamak, Plasma Devices Oper. 16 (3) (2008) 171–179. [3] R. Fresa, et al., Sensitivity of the diamagnetic sensor measurements of ITER to error sources and their compensation (in press on Fusion Engineering and Design). [4] R. Albanese, et al., Electromechanical analysis of end windings in turbo generators, Int. J. Comput. Math. Electr. Electron. Eng. 30 (6) (2011) 1885–1898. [5] D. Bloomberg, V. Castelli, Reformulation of nonlinear integral magnetostatic equations for rapid iterative convergence, IEEE Trans. Magn. 21 (2) (1985) 1174–1180. [6] R. Fresa, G. Rubinacci, S. Ventre, An eddy current integral formulation on parallel computer systems, Int. J. Numer. Methods Eng. 62 (2005) 1127–1147. [7] I. Mayergoyz, Nonlinear magnetostatic calculations based on fast multipole method, IEEE Trans. Magn. 39 (3) (2002) 1103–1106. [8] O. Chadebec, J.L. Coulomb, F. Janet, A review of magnetostatic moment method, IEEE Trans. Magn. 42 (4) (2006) 515–520. [9] Z. Haznadar, Z. Stih, S. Berberovic, The least square deviations criterion and inverse-source magnetostatic problems, Int. J. Appl. Electromagn. Mech. 15 (2002) 207–212. [10] L. Urankar, High accuracy field computation of magnetized bodies, IEEE Trans. Magn. 21 (6) (1985) 2169–2172. [11] A.G. Chiariello, A. Formisano, R. Martone, Fast magnetic field computation in fusion technology using GPU technology, Fusion Eng. Des. 88 (9–10) (2013) 1635–1639. [12] E. Durand, Magnetostatique, Masson, Paris, 1968. [13] H. Lei, W. Lian-Ze, W. Zi-Niu, Integral analysis of a magnetic field for an arbitrary geometry coil with rectangular cross section, IEEE Trans. Magn. 38 (6) (2002) 3589–3593. [14] N. Mitchell, Scenarios for Coil, Power Supply and Cryoplant Analysis, ITER, Rep. ITER D 2FTVKV, 2009.
Contents lists available at ScienceDirect
Fusion Engineering and Design journal homepage: www.elsevier.com/locate/fusengdes
Effective magnetic field computation in Tokamaks in presence of magnetic materials Andrea G. Chiariello a , Alessandro Formisano a , Raffaele Fresa b,c , Francesco Ledda a,∗ , Raffaele Martone a , Francesco Pizzo a a
Dip. di Ing. Industriale e dell’Informazione, Seconda Università di Napoli, Via Roma 29, I-81031 Napoli, Italy Scuola di Ingegneria, Università della Basilicata, Potenza, Italy c Consorzio EURATOM/ENEA/CREATE, Italy b
a r t i c l e
i n f o
Article history: Received 3 October 2014 Received in revised form 15 June 2015 Accepted 16 June 2015 Available online 10 July 2015 Keywords: Magnetic field evaluation Ferromagnetic materials Fusion Technology Neutral Beam Injector
a b s t r a c t Magnetic materials play an important role in the magnetic fields distribution inside Tokamaks, and their contributions must be carefully computed. In some applications, when 3D geometries are involved and high accuracy must be achieved, finite elements may reveal too computationally demanding and different approaches must be considered. In this paper, two possible solutions are presented. The first one is based on the discretization of magnetic parts into rectangular prisms with constant magnetization density; analytical formulas are then used to compute the field contribution due to each brick. The second one is based on the representation of the effect of magnetic materials as a set of dipoles, by using again analytical formulas to evaluate the contribution to the magnetic field. Iterative procedures are applied to evaluate the magnetization density in prisms. The efficiency and accuracy of both models are validated against commercial codes on a test case, and then applied to the computation of magnetic field in presence of Neutral Beam Injector in an ITER-like geometry. © 2015 Elsevier B.V. All rights reserved.
1. Introduction Accurate knowledge of magnetic fields in Tokamaks represents a key priority in the design of present devices, where size and power levels make the accurate control of fields mandatory to achieve desired performance. Discrepancies with respect to the design field map are originated not only by magnets fabrication and installation tolerances, by the joints and by the busbars, but also by the presence of ferromagnetic elements. Their combined effect can be very hard to be estimated using standard computational techniques, e.g. Finite Elements Method (FEM), since high accuracy is required, but no symmetry can be exploited to reduce the model complexity. A number of approaches were proposed to cope with deformed magnets, joints and busbars [1–3], but none of them, at authors’ best knowledge, is specifically tailored to cope with magnetic materials. A further example is the effect of ferromagnetic parts on magnetic diagnostics signals for plasma identification, where again the contribution from iron, eventually lacking any symmetry, must be accurately computed. Magnetic materials are present in Tokamaks, for example in the test blanket modules (TBMs) or in the TF ripple
∗ Corresponding author. E-mail address: [email protected] (F. Ledda). http://dx.doi.org/10.1016/j.fusengdes.2015.06.078 0920-3796/© 2015 Elsevier B.V. All rights reserved.
correction system, or in the passive shielding of Neutral Beam Injectors (NBI). Accurate field computation in large 3D structures in presence of magnetic materials has been widely discussed in the past. Different methods have been adopted, including integral methods [4,5], “smart” formulations implemented in FEM models to save memory and computational time (sparsification methods [6], Fast Multipole Methods [7]), or FEM-BEM formulations [8]. The development of FEM and BEM approaches did somehow put aside methods based on equivalent sources (see, e.g. [7,9,10]). Such methods may represent again a viable alternative thanks to the involved semianalytical formulas, relatively straightforward and well suited for implementation on architectures based on Graphic Processing Units [11]. In this paper, two approaches are described to reduce the computational effort in presence of fully 3D ferromagnetic parts. Both take advantage from discretization into elements possessing analytical expressions relating the field to the magnetic equivalent sources. Namely, a discretization into rectangular prisms with constant magnetization density and a representation using sets of magnetic dipoles are considered. The effectiveness of the methods is first tested against 3D FEM using a simple geometry, and then applied to the evaluation of the effect of the Neutral Beam Injector (NBI) passive shields in an ITER-like device on the magnetic
642
A.G. Chiariello et al. / Fusion Engineering and Design 96–97 (2015) 641–644
field at the first wall, where the diagnostics are placed. The paper is organized as follows. In Section 2 the basic mathematical models are discussed, in Section 3 a short description of the NBI passive shields is presented; Sections 4 and 5 show the results achieved in a couple of examples; finally, in Section 6 some conclusions are derived. 2. Mathematical model Although other cases could be treated in a similar way with minor modifications, for the sake of exposition in the following it is assumed that the time evolution of the magnetic field is much slower than the time constants of the passive structures. This is certainly valid during flat-top of the plasma discharges, but may be critical during disruptions or ramp-up and ramp-down. Under these assumptions, a non-dynamical model can be used to describe the flux density profile: div B = 0
(1)
rot H = J
(2)
If the field is required also when dynamical effects are relevant, then the model must be modified to include the full set of Maxwell equations. Work is in progress from authors to extend the computational codes to cope also with fast parts of the discharge. The magnetic constitutive relationship is assumed in the form B = 0 (H+M(H)), where M is the H-controlled non-linear function providing magnetization vector. It is assumed that the iron constitutive relationship is isotropic, and does not show hysteretic behavior. The whole analysis domain ˝ can be divided in four different regions: ˝plsm
J = Jpl M=0
Plasma region : the current is derived from
equilibrium equations
˝coil
J = Jext M=0
Active coils region : the current is imposed by
power supply system
˝iron
J=0 M = M(H)
˝iron
M=0
(3b)
Iron regions : the magnetization is derived from
the material characteristic
J=0
(3a)
Vacuum region
(3c)
(3d)
Note that, thanks to non-dynamicity assumption, passive structures can be neglected, since no eddy currents are present. Of course, a suitable set of boundary conditions have to be included on the external boundary of ˝0 . The computation of magnetostatic field in presence of magnetized bodies can be decomposed into three logical steps: 1. Evaluation of magnetic field Hsrc from impressed currents in free space; 2. Evaluation of magnetization distribution M in magnetic materials; 3. Evaluation of field contribution Hiron due to magnetization, in the free space.
Fig. 1. The elementary “brick” used for the iron parts discretization.
The first step can be done using Biot–Savart law [12], or any other method suited for impressed currents, and will not be discussed any more. The second step can be pursued by decomposing M on a suited representation basis, and solving the non-linear integral equation with a suitable tool, such as, e.g. the Picard approach. The best-known methods in this class represent the material with “equivalent” magnetic sources (magnetic dipoles or fictitious currents/charges [12]), and include the source parameters among the unknowns. This approach will be used here, discretizing the magnetic body into a suitable number of elementary subdomains, whose shape is chosen based on the availability of simple and accurate formulas to compute magnetization and magnetic field. Two possible analysis methodologies will be developed in the following: 2.1. Rectangular prisms (bricks) In the first approach, magnetic parts will be decomposed into rectangular prisms, “bricks”, and a uniform magnetization density is assumed in each brick as in Refs. [7,9]. The magnetic scalar potential due to an assigned magnetization map M(r ) can be written as: iron (r)
=−
1 4
∇ · M(r ) ˝iron
1 d˝iron |r − r |
(4)
where ˝iron is the region of space filled with magnetic material, r is the “field point” coordinate vector, and r is the “source point”. The contribution of the iron on the magnetic field is then computed as Hiron = iron . If M is constant over a “brick” aligned along one of the coordinate axes, (4) can be computed analytically in closed form, since the divergence operator is vanishing everywhere in the inner domain ˝iron of the iron domain ˝iron and just the contribution due to the discontinuity of the normal component of M must be considered. For the sake of exemplification, the magnetic potential z due ˆ as in the brick depicted in to a z-aligned magnetization M = Mz Z Fig. 1 is 1 (−1)k−1 4 k=1.8
z − z 2 k 2(z − zk ) atan + atan + 1 y − yk
z (r)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨
=
(x − xk )asinh
⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬
y − yk
(x − xk ) + (z − zk ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ x − xk ⎪ ⎪ ⎩ (y − yk )asinh 2 2
2
(y − yk ) + (y − yk )
+
(5)
⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ 2
where (xk , yk , zk ) represent the coordinates of brick vertices in Fig. 1, and
1 = (x − xk )2 + (z − zk )2 + (x − xk ) 2 = (y − yk )(z − zk )
(x − xk )2 + (y − yk )2 + (z − zk )2
A.G. Chiariello et al. / Fusion Engineering and Design 96–97 (2015) 641–644
643
Eq. (5) was derived starting from (4) by using the symbolic solver of Matlab© R14 to achieve best performance when using HPC algorithms and it has been compared to similar expression from Ref. [9] to assess correctness of resulting values. The magnetic field Hiron of a uniform magnetization along z in the brick will then be given by: Hiron (r) = Mz (hxz xˆ + hyz yˆ + hzz zˆ ) where hkz = ∂ z /∂k (k = x,y,z). Similar results of the scalar potential can be obtained also when the magnetization M is x or y-aligned. Since computation in this stage assume that magnetization map is assigned, it is possible to decompose any distribution of magnetization into the three components M = Mx xˆ + My yˆ + Mz zˆ , and obtain:
⎡
⎤
hxx
hyx
hzx
Hiron (r) = hM(r ) h = ⎣ hxy
hyy
hzy ⎦
hzx
hzy
hzz
⎢
⎥
(6)
where the tensor h links the effect of the uniform 3D magnetization of the brick to the magnetic field. It should be noticed that (6) could be easily generalized for a brick not aligned with the coordinates as in Ref. [13]. When more than one brick or field points are considered, each pair of them will be treated as in (6) and a suitable matrix of tensors h should be considered including three rows for each field point and three columns for each brick. Limiting in a Cartesian coordinate system, the tensors matrix can be suitably treated as 3 × 3 partitioned matrix. An iterative scheme based on collocation method is used to evaluate the magnetization of the bricks. First (6) is used at each iteration to obtain the magnetic field in the center of each brick from the magnetization map at previous iteration. Using the non-linear constitutive relationship M = M(H), a new estimate of M (constant over each brick) is then obtained, until convergence in Cauchy sense is achieved: Hnk+1 = Hnsrc +
n,m
m=1:Nbricks
h
Mm k
Mnk+1 = ˇkn M(Hnk+1 ) + (1 − ˇkn )Mnk
(7) (8)
n=1,2,...Nbricks
where Nbricks is the total number of bricks, k is the Picard iteration index, hn,m is the 3 × 3 submatrix of h characterizing the contribution of the m-th brick on the n-th one whose elements are computed according to (6), and ˇkn = 1/(1 − nk min ) is an underrelaxation factor, required to achieve convergence [7], k is the k-th iteration estimation of the relative permeability of the n-th brick, and min is the smallest eigenvalue of h. A Newton schemes may speed up the convergence of the nonlinear iterative scheme but they would require a matrix inversions at each iteration step. Using (7) and (8), very large magnetic structures may be treated by assembling “on the fly” the required matrix elements, saving memory space and possibly relying on GPU speed to limit the computing burden. 2.2. Equivalent dipoles The bricks model allows self-consistently computing magnetization map M and field due to magnetic parts Hiron . Anyway, once the magnetization is known, either using bricks or using FEM, the effect of magnetic materials can be taken into account in remote regions (e.g. in the plasma chamber) by using faster formulas for the “long distance” effect. In fact, the magnetic material can be partitioned into smaller volumes, e.g. spheres, and each sphere
Fig. 2. A schematic view of the NBI and vessel considered geometry. Vessel sections are used to define test points.
of radius rsp can be associated to a magnetic dipole of assigned momentum: m=
4 3 r M 3 sp
(9)
The value at the generic point r of the magnetic field associated to a magnetic dipole is: Hiron (r) =
−m + 3(m · l rr )lrr
(10)
4|r − r |3
where r is the point where the dipole is collocated and ˆirr is the versor oriented along the direction r − r . 3. Neutral beam magnetic field reduction system Two Heating (HNBI) and one Diagnostic (DNBI) Neutral Beam Injectors are foreseen in ITER and a third heating injector could be installed in a successive operation phase. The residual magnetic field inside the NBI region is required to be low enough to avoid deflection of the ion beam and to achieve the design performance of these components. Both passive and active shielding systems have been foreseen, to reduce the high stray field inside the NBI due to the plasma and to the Poloidal Field Coils to an acceptable level. The effect of NBI passive shielding anyway is not negligible also outside the injector: in this paper, the interest is on the field contribution in the vessel surface, where magnetic diagnostics are placed (Fig. 2). 4. Test case In order to assess the accuracy and the computational effectiveness of the two models, a simple test case is considered: a hollow rectangular prism of a magnetic material, with similar dimensions of the external NBI iron shield is immersed into a uniform Hsrc field, vertically directed along z-axis. The resultant magnetization map M has been compared among the results given by the two models and those provided by a commercial FEM. Details about geometry, magnetic field and discretization levels are reported in Table 1. The magnetic field H along the symmetry axis z has been compared Table 1 Details of the test case. Geometry of hollow prism External dimensions Wall thickness Magn. permeability
x = 13.8 m; y = 5.11 m; z = 5.03 m 0.15 m rel = 1000 (constant)
External field
Hsrc = 1000 AT/m along z direction
Details of models Bricks Dipoles FEM
363 on long sides, 121 on short, 1694 total Same as bricks 3840662 DoF
644
A.G. Chiariello et al. / Fusion Engineering and Design 96–97 (2015) 641–644
Fig. 3. The comparison between the FEM, dipoles and bricks.
ferromagnetic shield. The shield is supposed to be constructed using SS400 steel. The source field distribution Hsrc is due to Toroidal field coils, Central solenoids, poloidal field coils and plasma. Geometries and currents are similar to one of the ITER reference scenario, namely Start of Flat Top equilibrium in the 15MA inductive scenario [14]. To assess the effectiveness of the approach, the logarithmic variation of field log10 (||Biron ||/||Bsrc ||) on the vessel is computed. In Fig. 4 the norm of the magnetization M on the NBI shield walls is reported; the effect of the shield on the vessel first-wall is reported in Fig. 5 both as average and as maximum value in each poloidal cross section along the toroidal direction. 6. Conclusions Elementary sources with known analytical expression can provide useful tools, alternative to FEM codes, for the evaluation of the effect of ferromagnetic materials. In this paper, two methods, based on analytical expressions for the magnetic field due to ferromagnetic parts immersed in a field due to external sources have been proposed. Advantages in terms of accuracy and computational times have been assessed using NBI passive shields as examples of ferromagnetic parts. Work is in progress by authors to extend the model to cope with “fast” current dynamics, and possibly with hysteretic materials. Acknowledgments
Fig. 4. NBI magnetization map.
This work was partly supported by the EU under the Contract of Association between EURATOM and ENEA/CREATE (F4E-2008-OPE06-04-03(ES-AC) Task Order n .3 Analysis of the error field induced in the magnet system of ITER), and partly by the Italian MIUR under PRIN grant 2010SPS9B3. References
Fig. 5. Effect of magnetization on the vessel.
using a FEM model, bricks model and dipoles model, with M provided by bricks. The results are reported in Fig. 3 and they allow to state that the three different models are all effective in computing the required field at long distance from the source while the dipoles model can provide inaccurate result at the proximity of the source. On the other hand, the dipole model on serial C code is about 82 times faster than the brick model to compute the magnetic field due to an assigned magnetization map. 5. Analysis of NBI impact In order to test the methods in a realistic case, the effect of a single NBI passive shield was computed using the brick model. Since the interest in this work is to assess the effectiveness of the proposed methods, the geometric NBI details are not considered and a simple hollow rectangular prism is considered to model the NBI
[1] J. Knaster, et al., ITER non-axisymmetric error fields induced by its magnet system, Fusion Eng. Des. 86 (6) (2011) 1053–1056. [2] V. Amoskov, et al., Assessment of error field from solitary ferromagnetic elements located outside of ITER tokamak, Plasma Devices Oper. 16 (3) (2008) 171–179. [3] R. Fresa, et al., Sensitivity of the diamagnetic sensor measurements of ITER to error sources and their compensation (in press on Fusion Engineering and Design). [4] R. Albanese, et al., Electromechanical analysis of end windings in turbo generators, Int. J. Comput. Math. Electr. Electron. Eng. 30 (6) (2011) 1885–1898. [5] D. Bloomberg, V. Castelli, Reformulation of nonlinear integral magnetostatic equations for rapid iterative convergence, IEEE Trans. Magn. 21 (2) (1985) 1174–1180. [6] R. Fresa, G. Rubinacci, S. Ventre, An eddy current integral formulation on parallel computer systems, Int. J. Numer. Methods Eng. 62 (2005) 1127–1147. [7] I. Mayergoyz, Nonlinear magnetostatic calculations based on fast multipole method, IEEE Trans. Magn. 39 (3) (2002) 1103–1106. [8] O. Chadebec, J.L. Coulomb, F. Janet, A review of magnetostatic moment method, IEEE Trans. Magn. 42 (4) (2006) 515–520. [9] Z. Haznadar, Z. Stih, S. Berberovic, The least square deviations criterion and inverse-source magnetostatic problems, Int. J. Appl. Electromagn. Mech. 15 (2002) 207–212. [10] L. Urankar, High accuracy field computation of magnetized bodies, IEEE Trans. Magn. 21 (6) (1985) 2169–2172. [11] A.G. Chiariello, A. Formisano, R. Martone, Fast magnetic field computation in fusion technology using GPU technology, Fusion Eng. Des. 88 (9–10) (2013) 1635–1639. [12] E. Durand, Magnetostatique, Masson, Paris, 1968. [13] H. Lei, W. Lian-Ze, W. Zi-Niu, Integral analysis of a magnetic field for an arbitrary geometry coil with rectangular cross section, IEEE Trans. Magn. 38 (6) (2002) 3589–3593. [14] N. Mitchell, Scenarios for Coil, Power Supply and Cryoplant Analysis, ITER, Rep. ITER D 2FTVKV, 2009.