3D-FOX—A 3D Transient Electromagnetic Code for Eddy Currents Computation in Superconducting Magnet Structures: DTT TF Fast Current Discharge Analysis

Bulky metallic structures are needed in the toroidal field (TF) superconducting magnets for fusion applications to withstand the large Lorentz forces acting on the winding. The pulsed coil operation during a plasma scenario and the fast current discharge, or a plasma disruption event in off-normal operating conditions, cause transient magnetic fields, inducing eddy currents in the TF structures. The eddy currents generate heat in the structures heating in turn the winding pack, eroding the temperature margin: such power deposition is a key input for thermal-hydraulic (TH) analyses. However, the computation of eddy currents in fusion magnets is a challenging topic since a transient, fully 3D electromagnetic (EM) model is required. The EM problem is solved here by means of the finite element (FE) open source code FreeFEM++. First, the correct implementation of the EM problem is verified by means of suitable benchmarks against both simple analytical cases and the results obtained with state-of-the-art FE commercial codes on the DTT TF coil, used as a reference geometry. Then the EM code is applied to the evaluation of the magnetic fields and eddy currents induced in the same reference coil during the normal (static) and off-normal (transient) operation; the output of the EM analysis is used as input to the TH analysis carried out with the 4C code, aimed at computing the temperature margin evolution during the transient.


I. INTRODUCTION
Bulky stainless steel structures are adopted in superconducting magnets for the plasma confinement in nuclear fusion reactors [1], to withstand the large Lorentz forces resulting from the interaction between high currents (several tens of kA) and high fields (1-10+ T). Due to the transient magnetic fields, produced either by the normal operation of the pulsed coils or by off-normal transients (e.g. fast discharges The associate editor coordinating the review of this manuscript and approving it for publication was Diana Leitao . or plasma disruptions), eddy currents are generated within the structures. These currents generate heat by Joule effect heating in turn the winding pack (WP), contributing to the erosion of the temperature margin.
The 3D-FOX (3 Dimensional Finite element Open source Eddy Current Simulator) is a 3D transient electromagnetic (EM) finite element (FE) code developed using the open source software FreeFEm++ [2] aimed at the computation of the power deposited by the eddy currents. Thanks to the 3D nature of the model, the spatial distribution of deposited power can be evaluated with a high level of detail, differently from what was proposed in previous analyses where it was approximated with simplified equations providing only the evolution of the integral value of the power, with no information on its distribution within the coil structures [3]. Similar computations to those proposed in this work are presented in [4], in which eddy currents in the EAST TF casing have been computed, but adopting commercial tools (ANSYS) and without focusing on the TH aspect presented in this work.
In the first part of the paper the development of the code is described, followed by the presentation of three benchmarks that have been adopted for its verification.
A first application of the tool is also presented here, simulating a fast current discharge in a toroidal field (TF) coil of DTT (Divertor Tokamak Test) facility [5], aiming to reproduce the situation that will be tested in the DTT cold test facility in the next few years. All the Nb 3 Sn DTT coils will indeed be cold-tested in the dedicated Frascati Coil Cold Test Facility, being prepared at ENEA Frascati premises, to qualify the coil performance before their assembly in the tokamak pit.
The output of the EM model is used as input for the thermal-hydraulic (TH) simulations which are carried out with the state-of-the-art 4C code [6]. TH analysis of the fast current discharge in the DTT TF coil is presented in this paper too.

II. EM MODEL DESCRIPTION
The model has been developed to simulate both steady states and transients. The magnetostatic model is useful to obtain the magnetic field map during nominal operation of DC coils as the TF ones, while the dynamic one is needed for the eddy current evaluation.

A. MAGNETOSTATIC MODEL
A current flowing within a conductor generates a magnetic field as described by Ampère's law where B is the magnetic induction, µ the magnetic permeability of the considered materials, J is the current density that flows in the conductor, generating the magnetic field, ε is the electric permeability of the considered materials and E is the electric field. In steady state applications, the transient term is set to zero and the equation to be solved is simplified as only the current density remains at the right hand side (rhs), representing the driver of the problem (Eq. 2).
For the solution the A-formulation is adopted, expressing the magnetic induction as a function of the magnetic vector potential A.
Substituting Eq. 3 in Eq. 2 the final equation to be solved by the model is For the FE solution of (Eq. 4) the Galerkin weak form must be used. Taking v as test function and simplifying the expression making use of suitable integration techniques the equation to be solved in FreeFem++ is the following: In the magnetostatic case nodal elements have been adopted for the solution decreasing the size of the algebraic problem with respect to the traditional use of edge elements. In order to assure the uniqueness of the solution the Coulomb gauge must be introduced in Eq. 5 [2], obtaining the final equation to be solved Using the FE each variable is approximated as where φ k are the basis functions, depending on the type of FE selected, w k are real numbers and M is the number of degrees of freedom (DoF) of the finite element space, which depends on the type of FE adopted and on the mesh size. Substituting this expression within the Galerkin equation and using proper numerical integration formulae, already implemented in FreeFem, the differential equation is finally translated into an algebraic system of equation whose size depends on the number of DoF (one equation must be solved for each DoF).
To withstand the big size of the algebraic system, it has been solved with the conjugate gradient (CG) method already implemented in FreeFem.

B. ELECTRO-DYNAMIC MODEL
The electro-dynamic model aims at solving the same problem, but with a time dependent driver. In principle, this would require the solution of the full set of Maxwell equations, but, as reported by Albanese and Rubinacci [7] and by Biddlecombe et al. [8], the displacement current (∂(εE)/∂t in the Ampère law) can be neglected whenever the time scale of the transient is longer than the relaxation time of the electric charge τ r = ε/σ , which is true for the typical cases in which the model will be adopted, where the electrical conductivity (σ ) is of the order of the 10 6 [S/m]. Therefore, these kind of problems are normally solved simply coupling the steady state Ampère's law (Eq. 2) and the Faraday's law (Eq. 8).
In this case, both the magnetic induction, the current density and the electric field are time dependent, and the current VOLUME 10, 2022 density includes both the prescribed current (the driver of the transient) and the generated eddy currents. Recalling Ohm's law, the overall current density can be expressed as Applying to Faraday's law (Eq. 8) the magnetic vector potential definition (Eq. 3) it is possible to express the electric field as function of A and use it in the current density definition expressed in Eq. 9, obtaining Thus the eddy current equation can be rewritten with only A as unknown as Using the same strategies adopted for the magnetostatic case, Eq. 11 can be translated into the corresponding Galerkin formulation to be solved by the FE software.
For the solution of the eddy current problem, the use of edge element is suggested [9]; moreover this choice automatically satisfies the gauge condition assuring, by definition, the continuity of the tangential component of the field.
In this case also the time domain must be discretized: the backward Euler method was adopted for simplicity given that no stability problem were expected depending on the time discretization scheme.
The resulting algebraic system is solved, in the eddy current case, by MUMPS (MUltifrontal Massively Parallel Solver) [10]. This is a direct solver which performs a Gaussian factorization. This solver has been chosen in FreeFem given that it was easily adaptable for parallel computation, which has been preferred to sequential one for the large size of the problem.
Once the value of the magnetic vector potential has been computed at each time step, the magnetic field is automatically obtained simply from the potential definition (Eq. 3), and also the eddy currents are obtained from Once the eddy current distribution is known, applying Ohm's law, the local power density can be computed. Finally, by integration over a desired volume also the power value can be obtained.

III. MAGNETOSTATIC MODEL BENCHMARK
The 3D-FOX magnetostatic model has been adopted for the simulation of the steady state operation of the DTT TF coil both in stand alone and in tokamak configuration (the former configuration being relevant for the coil qualification tests in the cold test facility) in order to benchmark the results with those obtained by means of the state-of-the-art software COMSOL [11] in the same situation.
The tokamak configuration considers here only the contribution given by the entire set of 18 TF coils and neglects the contributions given by the Central Solenoid, the Poloidal Field coils and the plasma.
For the magnetostatic simulation, the structures can be neglected as they are not expected to play any role, since no current will flow in the casing contributing to the generation of the magnetic field. The WP has been modeled on the base of the geometrical information reported in [12]. The WP is considered as a single rectangular conductor with an external surface corresponding to the envelope of all the turns on which a uniform current density, whose integral is equal to the total operative current, is imposed. This solution simplifies the geometry, with respect to the separate model of each single turn, and still does not influence the evaluation of the magnetic field outside the winding pack.
The resulting geometry of the WP and its cross section are shown in figure 1.
The magnetic permeability of the superconducting cables is prescribed to be constant, uniform and equal to the vacuum one, namely µ = 4π · 10 −7 H/m. The air is modeled in different ways according to the physical situation to be simulated. In this work the stand alone configuration is simulated using a full sphere as external world whereas in the tokamak configuration 1/18th of the sphere is used, to exploit the symmetry of the torus by means of homogeneous Neumann boundary conditions (BCs) on the cutting surfaces. In figure 2 the two configurations described are shown.
The driver is the current flowing in the coil which is defined as components of the current density vector supposed to be uniformly distributed on the conductor section with value imposed such that the integral current is preserved, knowing that each turn carries I = 42.5 [kA].

A. RESULTS
The 3D nature of the tool allows to obtain detailed local information concerning the evaluated magnetic field. In figure 3 the magnetic field map on the vertical cross section of the TF coil is shown for both configurations at hand.
The stand alone configuration behaves like a single conductor spire with current, with uniform magnetic field in the coil center. On the contrary, in tokamak configuration the effect of the other coils in the TF magnetic system is evident in the non symmetric behavior of the field, showing a peak at the inboard leg where the effect of the other coils is stronger; in that region, indeed, all the inboard straight legs of the TF coils are close to each other.

B. BENCHMARK
For benchmarking purposes these results have been compared with those computed with COMSOL [11] at the center of each conductor on the equatorial plane. Both the radial and the toroidal distribution of the field amplitude have been compared for both analyzed configurations. The toroidal distribution has been compared along the first, the fifth and the ninth turn, while the radial one along the sixth, eighth and tenth pancake, both at inboard and outboard legs on the equatorial plane. As an example, the comparison of the toroidal distribution at the inboard leg is reported here (figure 4). The majority (96%) of the points shows a relative error between the 3D-FOX and COMSOL smaller than 0.75% and only two isolated points at very low field show a relative error of almost 1%. At high fields, the agreement between the two tools increases, with errors within ∼ 0.1%. Even though the agreement could be improved by further refining the mesh (and increasing the computational cost), this accuracy is considered sufficient for the application at hand.

IV. ELECTRO-DYNAMIC MODEL BENCHMARK
As for the magnetostatic model, also the electro-dynamic one has been benchmarked against other codes for verification purposes. In this case there are no available results on the same geometry adopted in the static case (the DTT TF coils), so two simpler test cases have been adopted. The first benchmark, taken from the literature, is the Felix brick problem (TEAM 4) [13]; the second one is a simplified coil model for which comparative results are obtained with the state-ofthe-art FE software COMSOL [11].

A. BENCHMARK 1: THE FELIX BRICK PROBLEM
The Felix brick problem is a classical literature benchmark in which a conductive aluminum brick with a central hole is subjected to a time varying magnetic field that is imposed as BC. For geometry, material properties and driver evolution see [13].
The same problem has been simulated with the 3D-FOX, imposing the time varying magnetic field as a correspondent time varying magnetic potential at the boundaries of the environment modeled as 1-m-side cube.
The integral value of the power deposited by eddy currents in the entire brick and the axial profile of the magnetic induction are the selected outputs of the calculation. The results obtained with the 3D-FOX are compared with those collected in [13], showing excellent agreement between different codes outcome as reported in figure 5. Indeed, the relative difference between the reference results and those obtained with the 3D-FOX for what concerns the peak power is of ∼ 3.6%(∼ 4.6 W) which translates into a relative difference of ∼ 1.4%(∼ 25 mJ) on the energy deposition during the transient. Moreover, the relative difference between the reference and the computed magnetic field is always well below 1%.

B. BENCHMARK 2: CIRCULAR CONDUCTOR WITH JACKET
The 3D-FOX is expected to work with a driver that is a time varying current that generates a time varying magnetic field. For this reason, another benchmark problem is proposed: a circular conductor surrounded by a stainless steel jacket in which eddy currents are generated due to the time-varying nature of the current J(t) defined in the conductor itself (Eq. 15).
where J 0 is the initial value of the current and t is the time.
The relevant input data of the simulation are presented in table 1. The integral value of the deposited power due to eddy currents in the jacket evaluated with the 3D-FOX is compared with that computed with COMSOL, showing a very good agreement, see figure 6. Indeed, a relative difference of ∼ 4.4% (∼ 0.67 mW) has been found on the peak power, which translates into a relative difference of ∼ 0.9% (∼ 0.11 mJ) on the energy deposited in the jacket during the transient.

V. ELECTRO-DYNAMIC MODEL APPLICATION -DTT TF COIL FAST DISCHARGE
After the verification, the model has been applied to the simulation of the fast current discharge in a single DTT TF  coil. This transient is relevant since it will be tested in the DTT cold test facility.

A. GEOMETRY AND MATERIALS
The considered geometry includes the TF WP, in which the driver (the time-varying current) is defined, and the coil structures, in which eddy currents are generated. A simplified version of the coil structures is proposed, see figure 7, neglecting all the tiny details, e.g., the bolts and the inter-coil structures.
The structure material is stainless steel and the magnetic permeability has been fixed for the entire domain to be equal to the vacuum value (µ = 4π ·10 −7 H/m), while the electrical conductivity has been assumed to be constant within the structures and equal to the value evaluated at 4.5 K, i.e., the initial temperature. The value is obtained from [14] and it is equal to σ = 2.0 · 10 6 S/m.

B. SIMULATION SETUP
The fast current discharge is a fast decrease of the coil current from its nominal value to zero. For the DTT TF, different evolutions are proposed. In this case, the exponential dump (with time constant τ = 5 s) is analyzed. The nominal value of the current density is obtained defining it uniformly The external air is modeled as a sphere with a 10 m radius on whose external surface ''far boundaries'' BCs have been imposed, namely prescribing that A ×n = 0. Making use of the edge finite elements basis function as weighting, this condition is satisfied by definition [9].

C. MESH GENERATION
The flexibility of 3D-FOX is enhanced by the mesh generation strategy. The meshes are generated with the open-source mesher GMSH [15], that takes as input any kind of geometry coming from any CAD software. The resulting mesh is then imported in FreeFem++, where the problem is solved.
The mesh has been generated with different refinements in the different regions, carefully discretizing the structures region in which eddy currents are generated. The different adopted refinement for structures, WP and external world are shown in figure 8.

D. RESULTS
A non-uniformity of the power deposition between inboard and outboard can be observed in figure 9, with a peak at the inboard, as well as the non-uniform power deposition within the same section, being the peak located at the back side (with respect to the plasma side) of the casing, close to the WP.
The non-uniformity of the power deposition was not easy to be predicted a priori due to the several different aspects that should be considered. Indeed both local electrical resistance and the magnitude of the induced eddy currents, which depends on the time derivative of the magnetic field, play a role in the power deposition computation.
In this case a higher power density is detected at the inboard leg section, where the casing cross section is smaller, resulting in a larger resistance for the same induced eddy current value. Indeed, the casing cross section at the inboard leg is ∼ 0.063 m 2 , while at the outboard leg is ∼ 0.096 m 2 . The simulation of a current discharge within a simplified D-shape coil with uniform casing section along the whole generatrix has also been performed (not shown here), confirming the dependence on the structure cross section and finding an almost identical power deposition at both equatorial inboard and outboard cross sections.
It may often be necessary to have the integral value or an average information on the power deposition (e.g. to define an adequate input for the 4C code). In figure 10 the evolution of the overall integral power deposited in the entire casing is reported. That evolution is in agreement with expectations, since it is directly correlated with the driver evolution. Indeed, combining Ampère law, Faraday law and Ohm law it is possible to state that P ∝ ( ∂J ∂t ) 2 ∝ e − 2t τ . So, a decreasing exponential behaviour with half the time constant of the driver is expected. Indeed the integral power evolution compares well with an exponential curve with the expected time constant, as shown in figure 10.
As described better in section VI-A, for the TH simulation the averaging of the 3D power deposition obtained with the FE tool in the structure segments must be provided as input. The tool itself is able to perform this averaging on arbitrary subvolumes to be defined by the user. As an example, the evolution of the average power density deposited in segments 3  (inboard leg, equatorial) and 11 (outboard leg, equatorial) in the single DTT TF coil is reported in figure 11.
The average power density in the two segments recalls what has been already observed from the power deposition maps of figure 9, with a lower power deposited in the outboard leg.

VI. FULL EM + TH APPLICATION: EVALUATION OF THE TEMPERATURE MARGIN EVOLUTION FOLLOWING A DTT TF COIL FAST DISCHARGE
The power deposition obtained as output of the electrodynamic model has been adopted as input to the 4C code to evaluate the temperature margin evolution during the same transient.

A. TH -EM COMMUNICATION
The outcome of the electro-dynamic model is used in this section as input for the thermal-hydraulic (TH) simulation, performed with the 4C code, aiming to evaluate the evolution of the temperature margin considering also the heat load contribution given by eddy currents.
The 3D results of the electro-dynamic tool must be adapted to be used in the 2D discretization of the coil structures adopted by the 4C code. Indeed to simulate the entire coil the 4C code subdivides the structures in different segments [16] (an example taken from the DTT TF is shown in figure 12a), each of them represented by its central 2D (radial-toroidal) cross section, with heat loads given by the average value on the segment itself. For this reason the average value of power density due to eddy current evolution is extracted, by means of Boolean controls, in each of the above-mentioned segments and used as input for the TH simulation. With respect to that reported in figure 12a, adopted in [17] for the nuclear heat load calculation with Monte Carlo tools, a different (e.g. more detailed) segmentation of the coil structures is possible, e.g. to separate sub-portions of the casing (plasma facing side vs. back side); however, this is beyond the scope of the present work. A schematic representation of the communication logic is shown in figure 13.

B. SIMULATION SETUP
Concerning the TH simulation setup, the strategy proposed in [18] has been adopted.
Here, the TH simulation has been performed with two alternative evolutions of the current, namely the exponential dump (discussed above) as well as a linear ramp discharge, see figure 14. The two different fast discharge solutions have been proposed since they are assumed to deposit, during the transient, the same amount of energy in the resistors. Indeed the integral of RI 2 (t) is equal in the two cases with a ramp length of 7.5 s, i.e. 3/2 of the exponential time constant (5 s).
The two different current evolutions are shown in figure 14.

C. RESULTS
The evolution of the power deposition by eddy currents depends on the kind of discharge. The integral power deposited in the entire coil casing is shown in figure 15.
The exponential dump shows a higher peak with respect to the linear ramp case, due to the higher time derivative of the coil transport current in the first instants of the transient; on the other hand, the linear ramp shows a constant power deposition given by the constant time derivative of the current in the WP.
Even though the two cases are assumed to deposit the same amount of energy in the resistors, this does not imply that the same energy must be deposited in the coil structures too; indeed, the energy deposited in the casing during the transient is higher in the linear ramp case of almost 35%, indeed ∼ 625 kJ are deposited during the linear discharge, versus ∼ 465 kJ deposited in the structures during the exponential discharge.
It is possible to explain this behavior considering the WP and the casing as two coupled RL electrical circuits, whose mutual interaction is described by the following system of   As it is possible to see from the second equation, the current induced in the casing (I c ) is a function of the time derivative of the current flowing in the WP (I WP ). As a consequence, with two different kind of discharge, the induced eddy currents in the casing (and thus also the power they deposit) are different. In particular, the time integral of the derivative of the linear ramp (with a duration 3/2 of the exponential time constant, as in the case at hand) turns out to be 3/2 of the same integral in the case of the exponential dump; this is not far from the comparison between the different energy deposition mentioned above. The power deposition causes a temperature increase in the structures. Assuming the structures as 0D and constant properties, the structure temperature increase in the linear ramp case should be ∼35% higher with respect to that resulting from the exponential dump. Looking at the 4C code results, this expectation is not respected, as it can be seen from figure 16. Indeed, the casing maximum temperature in case of linear discharge is higher, but with an increase which is smaller with respect to that estimated by lumped parameter assumptions. This depends on the fact that the lumped parameter model adopted is adiabatic and thus neglects the heat exchange with the He. The power deposited in the structures is partly removed by the casing cooling channels, and the remaining is transferred to the WP contributing, together with the AC losses, to erode the temperature margin.
The superconductor temperature increase is monitored by means of the hot-spot temperature, whose evolution is reported in figure 17, in both the analyzed cases, for pancake 1 and 6.
The substantial increase of the temperature is detected in the first few seconds of the transient, when the current is decreasing fast, thus depositing the largest eddy current losses in the structures and AC losses in the superconductor. The evolution of the hot spot temperature in pancake 1 and pancake 6 is different and this is due to the different way in which the two pancakes are thermally coupled to the casing. Indeed, pancake 1 is in contact with the casing for all its length, while pancake 6, which is in the middle of the WP, is thermally coupled to the casing only with its first and last turn. As a consequence, the power which is deposited in the first pancake due to the casing heating is higher with respect to that deposited in pancake 6 and this is translated into an higher temperature reached in the conductor. Moreover, in the central pancakes (P6 in the plot) the evolution shows two separate peaks, which are not detected in the side pancakes (P1 in the plot). Again, this is a consequence of the thermal coupling to the casing. Indeed, the first pancake receives from the casing a load which is almost uniform along its length and so the hot spot temperature evolution shows a single peak which is due to the heating and than re-cooling of the magnet. On the contrary the central pancakes experiment localized heating from the casing, on the first and last turns, and this causes the generation of two hot spots which are responsible of the above mentioned two separate temperature peaks: the first one is generated by the power which is deposited in the last turn and is therefore immediately appreciable, while the second one is the peak generated by the power deposited in the first turn and that is then transported along the entire pancake by the helium flow before exiting the coil. This thesis is furthermore confirmed by the delay between the two peaks, which is comparable with the helium transit time.
Finally it is possible to observe also that, coherently with the higher energy deposition, in the linear ramp case the hotspot temperature is higher with respect to the exponential dump one.
Comparing the strand temperature with the local current sharing temperature, the temperature margin is evaluated and the evolution of its minimum is plotted in figure 18. A zoom on the first 20 s is proposed to show that the strand temperature increase is slower than the increase of the current sharing temperature due to current (and magnetic field) reduction, and so a margin increase is actually detected during the discharge. Note that the temperature margin shown here is much larger than that expected in a TF coil operated in a tokamak (5 K vs. 1-2 K). This is due to the fact that in the CTF the overall magnetic field is smaller (only one coil is tested, and charged, at a time) and no nuclear heat load is present.
The minimum margin is never reaching, neither in the linear nor in the exponential dump case, the minimum requirement for DTT TF coils, fixed at 1.4 K [16]. A fortiori, this means that no quenches are foreseen for this transient, expected to be tested in the DTT Cold Test Facility. The evolution shows, as expected, the minimum at the beginning of the transient, when the magnetic field and the current density are still high and consequently the current sharing temperature is low. Then the margin increases reaching an asymptotic value corresponding to the maximum current sharing temperature (obtained for zero magnetic field and current) and the new steady state strand temperature reached after the dump.

VII. CONCLUSION AND PERSPECTIVE
In this work, the development of a first-of-a-kind 3D transient electro-magnetic tool for eddy current computation in the structures of bulky SC magnets for fusion reactors (the 3D-FOX) has been described, including its utilization as input for the TH simulations performed with the 4C code. The 3D-FOX is open source and extremely flexible, potentially applicable to any kind of electro-magnetic transients in any type of coil geometry.
The tool has been first successfully verified through suitable benchmarks against analytical and literature test cases, as well as against other reference tools (COMSOL). Then the qualified model has been adopted for simulation of the DTT TF coil fast current discharge. Two different options being considered for the current evolution during the dump have been simulated, both from the EM and TH point of view. The 4C code results, adopting as input the results computed with the new tool, show that no quenches are expected during the current discharge in the cold test facility, even though a different hot spot temperature is computed to be reached in the two different options.
The tool will be applied to several different transients, ranging from the eddy current generated in the coil casing by a plasma disruption or by the pulsed coil operation during the tokamak pulsed operation to those induced by similar transients in the conductor jacket. Other magnets will also be analyzed, such as the TF coils of the EU DEMO. Moreover, improvement will be implemented to optimize the EM+TH coupling as well as to include the impact of the temperature increase on the electric resistivity of the structures and, in turn, on the eddy current generation.

ACKNOWLEDGMENT
Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.