Fusion Engineering and Design EM zooming procedure in ANSYS Maxwell 3D

The severity of electromagnetic (EM) loads produced by plasma disruptions is one of the most concerning issues for the ITER in-vessel components design. To investigate the e ﬀ ects of fast EM transients on plasma surrounding structures during a disruption the Secondary Excitations (SE) method is used. This is an interface procedure to couple 2D plasma equilibrium codes with Finite Elements (FE) software. The Zooming Approach (ZA) used for the analyses presented here is a particular implementation of the SE method. The aim of this work is the de- monstration that the ZA can be e ﬀ ectively applied in case of ANSYS Maxwell 3D analyses combining the ease of use of the Maxwell code with the computational e ﬃ ciency of the ZA. The work has been carried out evaluating the EM loads acting on the ITER Diagnostic Equatorial Port Plug (EPP) during major disruptions scenario and comparing these loads with those obtained in previous analyses. Additional analyses have been performed to study the e ﬀ ect of ferromagnetic materials on EM loads in order to investigate ANSYS Maxwell capabilities in simulating non-linear magnetic properties.


Introduction
Electromagnetic (EM) loads produced by plasma disruptions have to be considered for a safe design of the ITER in-vessel components. The development of numerical methods providing the EM loads evaluation particularly during fast transients, is thus of meaningful importance to verify the robustness of ITER in-vessel components and to support their design.
According to the Faraday's law, when a volume of conductive material is immersed in a variable magnetic field, eddy currents will be induced in it. They will produce a magnetic field which opposes to the magnetic flux variation. The interaction of the eddy currents with the external field produces Lorentz's forces on the conductive volume. The EM loads analyses are performed with FE software and face essentially two main problems: the reproduction of the plasma behavior during disruption events; the modelling of the empty regions because of the field propagation through the vacuum.
The plasma behavior is analyzed in 2D equilibrium codes [1] in which the plasma is simulated by means of current filaments varying in intensity and position during time. Their direct use, in a FE analysis, would require a moving and detailed meshing of the plasma region. To overcome this problem, the Secondary Excitations (SE) interface method (described in [2]) is adopted to interface the 2D plasma equilibrium code outputs with the FE model.
The main drawback of modelling all the vacuum regions is instead represented by the high number of elements necessary for an accurate modelling of the components to be analyzed. This implies indeed high computational costs and problems in mesh connectivity between the conductive structures and the surrounding environment. In order to reduce the number of elements without losing results accuracy, a specific application of the SE method has been used: the Zooming Approach (ZA) [2].
In what follows the ZA is described through its implementation on the ITER Diagnostic Equatorial Port Plug (EPP) during major disruptions.

Zooming approach
To perform reliable EM analyses, the interface procedure (to adapt plasma equilibrium code outputs to FE fixed grids) must ensure a correct reproduction of the EM transient in terms of magnetic field and time derivative of the magnetic field. In doing that, the SE method replaces the axisymmetric filamentary conductors of the plasma equilibrium code, the so called Primary Excitations (PE), with an array of fixed axisymmetric conductors, the SE. The field and field time derivative of the original PE are reproduced in the FE model assigning specific currents (calculated through EMAG code [2]) to the fixed conductors assembly.
The SE method can be used with two different approaches (Fig. 1): a) the Global Approach (GA), where a full tokamak sector is modelled and the plasma region is completely included in the SE; b) the ZA, where the local model is enclosed in the SE array and the plasma region is outside.
In EM analyses of single components, the effects of the environment have also to be taken into account. In these cases the GA requires a hard modelling work of the surrounding environment as well as a big computational effort. The ZA instead allows a local analysis (from a whole domain to a region): only the sample and the main surrounding structures enclosed in the SE array need to be modelled.
The evaluation of the SE for the ZA is performed by EMAG code (described in [2] and used in ITER IO). The starting point is the assessment of the eddy currents induced by the plasma currents (obtained by DINA [1]) in the main toroidal continuous structures outside the sample. Subsequently, these eddy currents are added to the plasma source and used as a new PE set to evaluate the current in the SE. The calculated currents reverted in sign and imposed to the SE array will reproduce the EM transient inside the considered region. The details of the mathematical framework used for the SE evaluation can be found in Roccella et al. [2].

ZA in ANSYS Maxwell 3D
The aim of this work is the implementation of the ZA in ANSYS Maxwell 3D software [3] to investigate the Diagnostic EPP during a major disruption. The basic idea is to take advantage of the ZA computational efficiency in the user-friendly environment of ANSYS Maxwell 3D. ANSYS Maxwell is an interactive software, based on Maxwell's equations applied to a finite region of space which uses FE analysis (under appropriate boundary conditions and user-specified initial conditions) to solve electric or magnetic problems.
An accurate EM analysis requires often a very time-consuming meshing procedure because the mesh has to be finer as we get closer to the sample and with proper connectivity at limiting surfaces. ANSYS Maxwell relies on automatic meshing tools with limited control over the general meshing parameters (i.e. element maximum size or maximum number on specified geometries). In order to enhance (and fulfil) the analysis accuracy it is also possible to import from the Magnetostatic Solver a finer mesh refined by the software until the convergence criteria set by the user are fulfilled. Fig. 2 shows the model implemented with tetrahedron elements and reshaped to fit the SE zooming region. It represents a 20°sector of ITER Vacuum Vessel (VV) completed of Diagnostic EPP, Blanket Modules, Poloidal Field (PF) Coils and Central Solenoid (CS).
The Master/Slave boundary condition is imposed on the two lateral surfaces of the external vacuum volume to allow the modelling of only one sector of the periotic structure. The toroidal field component is in fact correctly reproduced by this boundary condition since it matches the field at the slave boundary to the field at the master boundary.  Because of the complexity of the model, computational schemes to satisfy (normal or tangential) continuity conditions of vector fields across the material interface are used. Accordingly, the components of a magnetic field that are tangential to the edges of an element are explicitly stored at the vertices; those tangential to the face of an element and normal to an edge are instead explicitly stored at the midpoint of selected edges. The value of a vector field at an interior point is interpolated from the nodal values; while the desired field in each element is approximated with a quadratic polynomial equation of 2nd order [3].
In the present work the ANSYS Maxwell Transient solver is used. The components modelled inside the zooming region are illustrated in Fig. 3. The VV and the Diagnostic EPP design have been simplified removing details reckoned to have relevant impact on the total element number count while negligible impact on the EM analysis results. The components dimensions and total volume remain almost unchanged; in this way, we can guarantee the correct evaluation of the EM loads (volume dependent quantities).
The Diagnostic EPP is the specific component on which the EM analysis is focused on. As shown in Fig. 4(a), it consists of the EPP structure, the Diagnostic Shield Modules (DSMs), and the Diagnostic First Walls (DFWs) [4]. The Diagnostic EPP model is shown in Fig. 4(b). Bolts, keys and features not relevant for the EM analysis are neglected. The electric contacts are reproduced by means of touching surfaces. The DSMs thickness is incremented compared to the design value in order to use the same conservative approach adopted in the Princeton Plasma Physics Laboratory (PPPL) EM analysis [5].

Current input data
The Diagnostic EPP EM loads have been evaluated for two major disruption events: Upward Major Disruption with 36 ms linear current quench (MD UP LIN 36 ms) and Upward Major Disruption with 16 ms exponential current quench (MD UP EXP 16 ms). Both the examined disruptions terminate with an upward movement of the plasma column.
The zooming SE and the PF and CS coils currents have been calculated elaborating the DINA output in the EMAG code. The current in the TF coil is imposed as constant and equal to the operating value of about 9.1 MA) corresponding to the reference field of 5.3T at the nominal major radius of the ITER plasma (6.2 m).
The SE array is composed of 72 toroidal volumes (in green in Fig. 2) enclosing the sample and replacing the plasma current and the vacuum vessel induced currents.
In order to guarantee that non-meaningful field variations appear at the beginning of the FE simulation, the EMAG current data have been modified for each SE and coil by imposing a constant current for the first three seconds of the EM transient.

Materials properties
The electric resistivity of the materials used in the EM analyses is reported in Table 1. The Diagnostic EPP components, the port and the VV are modelled assuming a 100% stainless steel structure.
Different electrical resistivity have been used for the blanket First Wall (FW) and the blanket Shielding Block (SB). An anisotropic electric resistivity is used for the blanket FW: the toroidal component is lower compared to the radial and poloidal ones, which are instead bigger due to the need to recreate the discontinuity caused by the fingers in these directions. For the definition of the anisotropic resistivity, a reference coordinate system has been created for the first wall of each blanket module.
For SE and coils, copper has been used.

Mesh sensitivity analysis
The meshing procedures discussed above have been taken into account to investigate their effects on the solution.
Firstly the Length-based Inside-selection refinement has been selected to limit the edge length of all tetrahedrons generated automatically by the Transient Solver. The maximum tetrahedron edge length in the Diagnostic EPP has been defined in order to be less than the skin depth. The resulting mesh is composed of about 1,350,000 elements.
Secondly the FE mesh has been obtained importing the adaptive mesh from the Magnetostatic solver. The mesh refinement is based on  the user required energy error that, for the present analyses, has been chosen equal to 0.01%. Because of the complex geometry, an initial mesh has been created through the definition of some tetrahedrons edge length limits aiming to increase the mesh density in the areas of interest. The mesh scheme obtained with the last procedure is composed of about 750,000 elements.
Comparing the Diagnostic EPP moments calculated by the two different mesh approaches in Fig. 5, it is possible to note that the difference in values interest mainly the radial moment and it is always below the 10%. For this reason, the imported mesh has been preferred to reduce the computational effort.

EM analyses results
During a MD event, eddy currents are induced in the Diagnostic EPP. The current loops lie mainly on horizontal and vertical planes because the major changes involve the poloidal field component. Fig. 6 shows the eddy currents density induced in the EPP at the end of the MD UP LIN 36 ms event. The currents are induced in loops to compensate the vertical and horizontal components variation of poloidal field.
It is worthy to note that the position of the TF coil in the model of the present analysis does not allow the reproduction of the ripple effect. Fig. 7(a) and (b) shows the results of the ZA implementation in terms of EM forces and moments acting on the Diagnostic EPP. The EM loads are evaluated for both the MD events referring to a local coordinate system, illustrated in Fig. 2, with origin placed in the center of the EPP structure flange. The radial force is the dominant component; its absolute value reaches the maximum at the thermal quench. In the MD UP EXP 16 ms the radial force becomes zero about 10 ms in advance compared to the linear disruption case. This is mainly due to the fact that the exponential law governing the current quench implicates a more rapid field variation at the beginning of the current decay, hence, a more important EM force in this phase.
The EM moments are reported in Fig. 7(b). The radial component is the biggest one and the maximum is reached in proximity of the end of the current quench. The maximum radial moment is bigger in case of exponential current decay again because the EM transient is faster in the first phase of current quench.

ZA and GA results comparison
In order to evaluate the computational efficiency of the ZA, the results of the previous EM analyses have been compared with a study carried out by PPPL [5].
In the PPPL analysis, the GA has been implemented in ANSYS Maxwell modelling a 20°sector of the ITER vacuum vessel with the SE array placed around the plasma region. The TF coils are properly modelled. Concerning the mesh, the number of elementary volumes in the GA model is much higher than in the ZA: the difference is of 75%.   The results comparison has been made for the MD UP LIN 36 ms case. The EM moments acting on the Diagnostic EPP are compared in Fig. 8. The radial moment has comparable values in both cases while the poloidal and toroidal moments have the same trends but lower values in the ZA. In correspondence of the first moment peak, the different values are justified by a different choice of the time steps at the thermal quench; in the ZA they have been tightly defined with the aim to capture the first fast field variation. The difference in the poloidal moment values is instead related to the absence of the ripple effect in the zooming model. In fact the presence of a TF coil in the middle of the sector creates a field axisymmetry responsible of the generation of forces with null poloidal moment if calculated respect to the Tokamak axis.

ZA and ferromagnetism
ANSYS Maxwell allows the accurate reproduction of a ferromagnetic material behavior simply importing the B-H curve material or defining the main curve parameters.
In order to evaluate the software capabilities in case of ferromagnetic effects, the previous model has been implemented in ANSYS Maxwell 3D. In this case, the DFWs of the EPP are assumed to be composed of EUROFER [6], a ferritic − martensitic steel which will be used for structures in EU test blanket modules. The ferromagnetic properties have been defined in ANSYS Maxwell 3D introducing the EUROFER B-H curve. Fig. 9 illustrates the EM forces acting on the ferromagnetic DFWs in case of MD UP LIN 36 ms. Due to the high background field, the ferromagnetic material reaches the saturation magnetization. The radial force component pushes both the DFWs towards the plasma because, according to the Maxwell's force [7], the radial component is directed as the toroidal field gradient. The Maxwell's force depends only by the magnetization and the external magnetic field, hence it is possible to note that high forces appear already before the transient starts.
The poloidal and toroidal components in upper and lower DFWs are instead opposite in sign due to the mutual repulsion. A ferromagnetic material in fact aligns its magnetic domains to the external magnetic field producing in this way a bipolar field inside the material. In the present case, the two DFWs are polarized in the same manner such that they are subjected to a mutual repulsive force.
The field amplification in a ferromagnetic body is also responsible of EM loads that are much bigger compared to those calculated on the same components in absence of ferromagnetic materials. Fig. 10 shows the difference in terms of forces values, which can be of some magnitude orders. It is also possible to note that the ferromagnetic material magnetization produces important forces already before the start of the EM transient.

Conclusions
The EM analyses performed in this work demonstrate that many advantages in EM loads evaluation arise from the implementation of the ZA in ANSYS Maxwell 3D. The Zooming represents a promising numerical procedure for the definition of EM loads acting on all the invessel components that is fundamental to guarantee the structural resistance during the tokamak operation.
The comparison with the GA implemented in the PPPL's study confirms that the zooming technique allows important time saving both in FE model preparation and in calculation time. However the absence of the ripple effect is responsible of some differences in the EM loads. For this reason, further developments of the ZA model are already planned, aiming at investigating the possible solutions for the correct reproduction of the ripple effect.
Concerning ANSYS Maxwell, it is demonstrated to be a very userfriendly software whose major advantage is the simplicity of the mesh   generation procedure, which allows a great reduction of the computational time.
The physical consistency of results proves that ANSYS Maxwell could be a very powerful tool also for the evaluation of the EM loads on the ITER Tokomak ferromagnetic components and the investigation of the ferromagnetic shielding efficiency.

Disclaimer
The views and opinions expressed herein do not necessarily reflect those of the ITER Organization.