Superconductor Science and Technology

PAPER

User defined elements in ANSYS for 2D multiphysics modeling of

superconducting magnets
To cite this article: Lucas Brouwer et al 2019 Supercond. Sci. Technol. 32 095011

View the article online for updates and enhancements.

This content was downloaded from IP address 70.115.114.155 on 28/10/2020 at 18:09

 Superconductor Science and Technology

Supercond. Sci. Technol. 32 (2019) 095011 (12pp) https://doi.org/10.1088/1361-6668/ab2e63

User defined elements in ANSYS for 2D

multiphysics modeling of superconducting
magnets
Lucas Brouwer1 , Diego Arbelaez1, Bernhard Auchmann2,3,
Lorenzo Bortot3 and Edvard Stubberud3
1
Lawrence Berkeley National Laboratory, Berkeley, CA 94720 United States of America
2
Paul Scherrer Institute, CH-5232 Villigen, Switzerland
3
CERN, CH-1211 Geneva 23, Switzerland

E-mail: [email protected]

Received 11 March 2019, revised 20 May 2019

Accepted for publication 1 July 2019
Published 2 August 2019

Abstract
Dynamic simulation of superconducting magnets is critical for the design of quench protection
systems to prevent potentially damaging temperatures and high voltage from developing after
magnet quench. Modeling these scenarios is challenging due to the many multiscale phenomena
which impact magnet behavior. These range from conductor scale effects of quench and
interfilament coupling currents up to the behavior of the magnet in its powering and protection
circuit. In addition, a strong coupling between electromagnetic and thermal domains is required
to capture temperature and field dependent material properties and quench behavior. We present
a finite element approach which integrates the various effects into the commercial software
ANSYS by means of programming new element types. This is shown capable of simulating the
strongly coupled transient electromagnetic, thermal, and circuit behavior of superconducting
magnets required for quench protection studies. A benchmarking study is presented which shows
close agreement between the new ANSYS elements and a COMSOL Multiphysics
implementation developed at CERN for dump resistor and coupling loss induced quench based
magnet protection of a Nb3Sn block dipole. Following this, the ANSYS implementation is
shown reproducing strongly coupled quench back behavior observed during the test of a Nb3Sn
superconducting undulator prototype at Lawrence Berkeley National Laboratory.

Keywords: superconducting magnets, multiphysics modeling, finite element, quench protection,

superconducting undulators
(Some figures may appear in colour only in the online journal)

1. Introduction material, current sharing and quench propagation within the

conductor, and inductive and resistive coupling of the various
Transient behavior of superconducting magnets is frequently effects to the magnet’s QPS circuit.
determined by multiscale and multiphysics phenomena. A Previous work has focused on simulating these challen-
common example of this occurs when a quench protection ging multiphysics problems using laboratory developed finite
system (QPS) is activated in an attempt to safely bring down element codes [1–4], lumped circuit element models [5],
the magnet current. To design a system which prevents customization of the commercial software COMSOL Multi-
damage to the magnet, it is critical to be able to accurately physics [6], and the coupling of multiple codes or softwares
simulate the temperature and voltage rise during the current using co-simulation [7, 8]. We present a new approach using
decay. This requires modeling phenomena such as quench only the commercial finite element software ANSYS [9]. It
back due to eddy currents in the conductor or structural has previously been demonstrated that ANSYS has the

0953-2048/19/095011+12$33.00 1 © 2019 IOP Publishing Ltd Printed in the UK

Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Figure 1. An overview of coupled electromagnetic, circuit, and thermal simulation in ANSYS with user defined elements is shown. Such an
approach allows for simulating the impact of interfilament coupling loss, quench, and structural eddy currents on magnet behavior while
including temperature and field dependent material properties. The independently meshed electromagnetic and thermal domains are coupled
using the Multi-field Solver as described in section 3.

capability to simulate some aspects of the complete problem, data taken during the test of a prototype Nb3Sn undulator
such as quench propagation [10–12] and the effect of eddy magnet at Lawrence Berkeley National Laboratory.
currents in mechanical structures on the current decay of a
circuit coupled magnet [13]. We show user defined elements
replicating these features and adding the additional missing 2. The finite element model
capabilities of: (1) modeling magnetization of the conductor
due to coupling currents and (2) combining all the effects into An electromagnetic and thermal model was developed for use
a single, coupled simulation with field and temperature in conductor regions where superconducting effects are
dependent material properties. desired (as illustrated in figure 1). These models were
The ability for a user to define their own element type is a implemented in ANSYS by the creation of two user elements.
documented feature of ANSYS, for which the authors are This approach integrates the desired effects at the point of
aware of two previous examples relevant to electromagnetic element matrix generation, no longer requiring manual
applications [14, 15]. The creation of a user element is updating of superconducting properties between a stop and
accomplished by writing the Fortran code which defines the restart of the solver as implemented in previous work [11, 12].
element’s properties and builds the finite element matrices. A The thermal model follows directly from ANSYS (see doc-
custom ANSYS executable is then compiled which allows for umentation in [17]), with extended capabilities of program-
the use of this element as if it were included in the standard mable material property fits and automatic quench checking
distribution (making it compatible with all geometry genera- and heat generation. The material properties are homogenized
tion, meshing, solving, and post-processing features). User during element matrix generation based on specified fractions
programmed generation of the element matrices gives full of conductor, superconductor, and insulation using a method
control over the choice of element shape functions, integra- similar to what is described in [6]. The electromagnetic model
tion points, material properties, and FEM formulation. is based on the vector potential approach used in ANSYS,
Two elements are used to implement a custom FEM with modifications to the formulation made to include the
approach, with the first being an electromagnetic element with effects of quench and IFCL. The following subsections
optional coupling to an external circuit. An equivalent mag- describe this model.
netization term is included in the vector potential formulation
to model interfilament coupling loss (IFCL) within the con-
2.1. Vector potential with equivalent magnetization for IFCL
ductor. A second, thermal element is used to model the
temperature rise due to quench induced ohmic heating and The default approach to modeling eddy currents in ANSYS
IFCL. These two elements are coupled using the Multi-field uses the A, V- A formulation with both vector potential A and
Solver in ANSYS such that magnetic field, temperature, and electric scalar potential V degrees of freedom in conducting
various loads are shared between them (see figure 1). We regions [18], requiring modeling and meshing of the con-
present an initial benchmarking study with results for a Nb3Sn ductive paths in which the induced currents flow. This is
dipole magnet compared between this approach and a similar impractical for the simulation of interfilament coupling cur-
implementation in COMSOL Multiphysics [16]. Following rents in a 2D magnet cross section due to filament sizes on the
this, behavior predicted by the user elements is compared to order of 5–100 μm and the 3D nature of the induced current

2
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

path. An alternative approach, which does not require a mesh the resulting element matrices and load vector are
density on the same order of the induced currents, assumes a 1
predetermined current path producing an equivalent magne- [K AA] =
m0 ò [( ´ [NA]T )T ( ´ [NA]T )] dVelem (8 )
tization, with time constant τ, of
2t ¶B 2t
Me = -
m 0 ¶t
. (1 ) [C AA] =
m0 ò [( ´ [NA]T )T ( ´ [NA]T )] dVelem (9 )

This approach is found applied to the case of a twisted fila-

mentary composite in a uniform, changing transverse field
{J S} = ò {JS}[NA]T dVelem. (10)

(relevant for 2D simulations of multifilamentary super- Here the equations are given in general 3D form in terms of
conducting strands) in [19, 20] and has been implemented in the nodal potential (assuming {A} = {Ax , Ay , Az }). For the
several magnet modeling codes [1, 5, 6]. With assumptions 2D user element these simply reduce to a single component
about the filament layout within the strand and the resulting Az, and are evaluated using Gaussian quadrature with similar
current loops, an IFCL time constant of shape functions and integration points as used by ANSYS.
m 0 ⎛ L ⎞2
t= ⎜ ⎟ , (2 ) 2.2. Coupling to an external circuit as a stranded conductor
2ret ⎝ 2p ⎠
Regions modeled as coils are coupled to an external circuit
is written in terms of an effective transverse resistivity of the
following the method developed by ANSYS [22], with sev-
strand matrix ρet and filament twist pitch L along the length of
eral modifications to account for IFCL, quench effects, and
the strand. Limitations of this approximation and a more
separate effective lengths for coil resistance and inductance.
detailed approach can be found summarized in [21]. The
A stranded formulation is used which adds a current i and
induced currents deposit energy as heat within the strand
voltage e DOF to the vector potential given in equation (5).
matrix with a power per unit volume of
Both these new DOF are constrained to be single values for a
¶B modeled coil region, with i being the current per stranded coil
Pe = Me · , (3 )
¶t turn, and e being the voltage drop across the coil. The voltage
drop is made up of both a resistive eR and inductive eL
which in many cases leads to IFCL being an effective quench
contribution such that
back mechanism.
The equivalent magnetization approach includes the e = eR + eL . (11)
effects of eddy currents without the need for an additional
degree of freedom (DOF), and the finite element formulation For a general stranded coil of fixed cross section and resis-
in 3D is derived from the vector potential only tivity ρ, the resistive voltage is given by

 ´ n  ´ A -  ´ Me = Js , (4 ) ⎛ N ⎞2
eR = i ⎜ c ⎟
⎝ Sc ⎠ ò rdV , (12)
with Js as a source current density and n = m-1. Considering
the form of Me, the differential equation to be solved using
where Nc is the number of coil turns and Sc is the modeled
the FEM is
coil area. As will be described in section 2.3, quench and
¶A current sharing effects are accounted for using a single
 ´ n  ´ A +  ´ 2tn  ´ = Js. (5 )
¶t parameter Ifcu representing the fraction of current assumed to
be flowing in the stabilizing material compared with the
Here it is seen the addition of the magnetization term intro-
superconductor. With this assumption, the contribution to the
duces a damping matrix of similar form (curl–curl) as the
resistive voltage is determined by the fractional area of sta-
stiffness matrix used for a typical vector potential element. To
bilizer and the amount of current assumed to be flowing in
implement the FEM, the weak integral of equation (5) is taken
this region. Adjusting equation (12) to account for these
with test functions chosen to be the same as the shape func-
assumptions, the resistive voltage drop is given by
tions carrying the DOF within the element [18]. This leads to
¶ ⎛ N ⎞2 I fcu
[K AA]{A} + [C AA]
¶t
{A} = {J s}, (6 ) eR = iL c ⎜ c ⎟
⎝ Sc ⎠ fcond (1 - fsc ) ò rst dS, (13)

from which the element stiffness matrix [K AA], damping

where Lc is introduced as an effective length chosen to best
matrix [C AA], and load vector {J s} are extracted. With the
match the resistance of the 3D coil, fcond defines the fraction
vector potential within the element written in terms of the
of the total coil area which is conductor, fsc defines the
shape function matrix [NA] and the DOF values at the element
fraction of superconductor within this conductor area, and ρst
nodes {Ae }
is the resistivity of the stabilizing material. In case of multi-
strand cables the factor Lc is also accounting for the trans-
{A} = [NA]T {A e }, (7 ) position pitch.

3
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Table 1. Quench states. The conductor is assumed to have a critical surface

parametrized by magnetic field, temperature, and transport
Regime Fraction of current in stabilizer
current density beneath which it is fully superconducting.
Fully superconducting Ifcu=0 Strain dependence of this critical surface is neglected at this
Current sharing 0 < Ifcu < 1 time. Above the critical surface, some (current sharing) or all
Fully quenched Ifcu=1 (fully quenched) of the current is assumed to move from the
superconductor into the stabilizing material, leading to resis-
tive loss. Accurate assumptions for the variation of Ifcu during
current sharing is dependent on the superconducting material
The inductive voltage is given by the time derivative of and the significance of this quench state to magnet behavior
the linked flux being simulated. An exhaustive review of this topic is found
¶F ¶ in chapter 18 of [23]. A common assumption for the current
eL =
¶t
= Nc
¶t ò (tˆ · A) dS, (14)
sharing regime of Nb3Sn and Nb–Ti conductor is a linear
variation of Ifcu with temperature T
which for a 2D model with current restricted to flow only in
the z direction is given by TcB - T
I fcu = 1 - (18)
N ¶Az TcB - Tcs
eL = L i c t
Sc ò ¶t
dS. (15)
from the point at which the critical current equals the trans-
port current Tcs up to the temperature TcB at which the
Here the variable t is positive or negative one based on the
transport current is equal to zero in the superconductor [24].
direction of the current, and Li is introduced as an effective
length scaling of the inductance chosen to match the 3D
2.4. Homogenized joule heating
magnet.
Circuit coupling is accomplished by the addition of The two sources of element heating are resistive loss due to
equation (11) to the original vector potential formulation in quench and IFCL. These losses are homogenized to account
equation (6) (with the source current density now derived for fill factors of conductor and superconductor within the
from the current per turn DOF using the proper winding modeled coil region. The parameter fcond is used to define the
function relations in the K Ai matrix). This leads to the coupled fraction of the total coil area which is conductor, and fsc to
equations define the fraction of superconductor within this conductor
¶ area. Quench induced loss is assumed to occur only within the
[C AA] {A} + [K AA]{A} + [K Ai ]{i} = {0}, (16) stabilizer of the conductor, with the magnitude dependent on
¶t
the quench state (see table 1) and the resistivity of the sta-
¶
[C eA] {A} + [K ee]{e} + [K ei ]{i} = {0}, (17) bilizing material. If Je is the element current density, the
¶t power per unit volume of modeled conductor is then given by
where the matrices shown are matched to the DOF depend-
(I fcu Je )2
ence of each term. With respect to the original ANSYS Pres = rst (19)
model, a new matrix C AA is included for IFCL, and the fcond (1 - fsc )
voltage balance matrices K ei and C eA are derived from
for quench based loss, and
equations (13) and (15) to include the changes mentioned.
2
The discrete form for the circuit coupled element (similar to fcond dB
Pe = 2t (20)
equations (8)–(10)), and the method by which the e and i m0 dt
DOF’s are coupled to other external circuit elements can be
found in the appendix. for losses resulting from interfilament coupling currents (see
equation (3)). For a given quench state, the temperature and
field dependence of these losses is driven by the variation of
2.3. Current sharing and quench
the stabilizer’s resistivity ρst. For Nb–Ti and Nb3Sn this sta-
The stranded formulation assumes a uniform current density bilizer is typically a high RRR copper, resulting in ρst (and
and therefore does not solve for the distribution of current then also τ for which ρet in equation (2) also depends on ρst)
within the superconductor or stabilizing regions which showing strong temperature dependence and magnetoresistive
changes based on the superconducting properties and quench effects.
behavior (this could be done, for example, with A–V or H
formulations). The choice of a stranded approach is motivated
by modeling of multifilamentary conductors typical of Nb–Ti 3. Coupled simulation using the Multi-field Solver
and Nb3Sn magnets. In order to capture the impact of the
effect of quench on ohmic loss and resistance growth, a single The ANSYS Multi-field Solver allows for solving of
scaling parameter Ifcu is used to represent the fraction of sequentially coupled problems with independent meshes. A
current assumed to be in the stabilizing material. With this unique, meshed region is generated for each physics field
definition, three distinct quench states can be defined as and load coupling interfaces for which loads will be
shown in table 1. passed between them are specified. Each region is solved

4
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Figure 2. A stagger loop within the Multi-field Solver is shown for coupled electromagnetic and thermal fields (see figure 1 for an example of
how such a simulation is set up with the user elements). In this example, the loads transferred between fields are heat generation f HGEN and
temperature f TEMP. This approach loops over the stagger time step (from ts to tf) with a relaxation factor α applied to the load transfer until
convergence of the loads is achieved. Separation of the problem into sequentially defined stagger steps is used to simulate over the entire time
domain.

independently with its own time stepping and solution accelerator magnets. This model is shown in figure 3, and a
options. The solver transfers the loads across the defined list of high level parameters are given in table 2.
interfaces (even with dissimilar meshes), and iterates between
each physics field in sequence until the transfer of loads
converges for a user defined ‘stagger’ time step as shown in 4.2. Protection with a dump resistor
figure 2.
This solver has been successfully used for fully coupled A first comparison between ANSYS and COMSOL was
simulations including the user elements (for example see the performed for a dump resistor extraction exhibiting strongly
verification study in section 4). To do this, two physics fields coupled electromagnetic and thermal behavior. A simple
are created which are shown labeled as ‘electromagnetic’ and circuit was built using CIRCU124 elements consisting of a
‘thermal’ in figure 1. A load transfer interface is specified resistor, voltage source, and a stranded coil element coupled
between meshed coil regions and any structural regions with in e and i to the coil region in the meshed model (see
eddy currents. This allows for passing Joule heat loads from figure 3). A dump resistor value of 30 mΩ was used for all
the electromagnetic region to the thermal region, and passing tests. The simulation begins with the magnet operating in a
temperature back. Both temperature and Joule heating are static condition at 4.5 K and 13.8 kA. This operating point is
standard loads which may be transferred with the Multi-field slightly less than 90% of the magnet’s short-sample limit. At
Solver. To allow for thermal material properties to also vary 5.0 ms the voltage source is ramped down to zero over 0.1 ms
with magnetic field and quench state, a workaround using a to effectively put the magnet in series with the dump resistor
shared module was implemented to pass these two variables only. The full details of this simulation including the material
in a non-standard use of the solver. property fits are found in [26].
The behavior of the magnet was studied with increasing
levels of detail as outlined in table 3. A first simulation was
performed with no IFCL or quench losses, making the current
4. Benchmarking with COMSOL decay dependent only on the magnet’s inductance. A second
simulation added IFCL which influences the current decay by
A verification study was completed comparing results from changing the differential inductance of the magnet. For this
the ANSYS user elements to a similar 2D FEM imple- case, coupling to a thermal model was also included to cap-
mentation in COMSOL developed at CERN [6] within the ture effects on τ due to changes in material properties from
STEAM project [25]. The results from the full study are heating of the conductor. For the final simulation, current
found in [26]. Effects such as quench resistance, yoke sharing and quench were added to the previous case, allowing
saturation, IFCL, and structural eddy currents were compared for IFCL induced quench resistance growth.
across several models with good agreement found. We pre- Agreement between ANSYS and COMSOL results for
sent the results of one such study which focuses on a sim- the magnet current decay is shown in figure 4 for the final
plified Nb3Sn dipole magnet protected with a dump resistor simulation case outlined in table 3. This magnet exhibits
and coupling loss induced quench (CLIQ) [27]. strong quench back behavior with IFCL heating the coil to
quench, after which the coil resistance growth rapidly
increases rate of magnet current decay. This is further illu-
4.1. The Nb3Sn dipole model
strated by figure 5 where the rapid growth of hotspot temp-
A dipole model was designed to allow for comparison of erature and coil resistance in the final, fully coupled case
results in a regime representative of realistic Nb3Sn (IFCL w/Quench) is compared between codes. The energy

5
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Figure 3. The two conductor layers and iron yoke of the Nb3Sn dipole ANSYS model are shown with the dump resistor circuit (a) and the
CLIQ circuit (b). The DOF of the meshed conductor regions of these models are coupled to their respective circuits, allowing for consistent
simulation of both electromagnetic and circuit effects. A similar model consisting of only the conductor region is used for the thermal domain
during the Multi-field solution.

loss for each mechanism is compared to the total change in

Table 2. Overview of the Nb3Sn dipole model. energy of the system in table 4.
Parameter Value Unit
Inner layer turns (per quadrant) 14 4.3. Protection with CLIQ
Outer layer turns (per quadrant) 18
Turn width 1.5 mm A second study compares results for a CLIQ discharge [27] in
Turn height 15 mm a layer-layer configuration as shown in figure 3. CLIQ is a
Strands (per turn) 40 protection scheme which seeks to quench large portions of the
Strand diameter 0.75 mm magnet coil in a distributed fashion in order to rapidly bring
Filament twist pitch 14 mm down the magnet current using coil quench resistance. At the
feff with ρet=feff ρst 1.0 detection of a quench, a capacitor bank is discharged across
Cu RRR 200 one or more sections of the magnet coil which generates an
Non-Cu fraction 0.4 oscillation of current in the coil sections about the nominal
Nb3Sn Jc (4.5 K, 12 T) 2040 A mm−2 magnet transport current (see figure 6). The resulting field
Short-sample current (4.5 K) 15.37 kA
oscillations induce interfilament coupling currents which heat
Short-sample cond. field (4.5 K) 11.7 T
Lc: effective res. coil length 10.11 m
portions of the coil to quench. This approach is typically
Li: effective ind. coil length 9.2 m considered for large inductance accelerator magnets where
dump resistor based protection is no longer possible due to
voltage and hotspot temperature constraints. CLIQ protection
scales well to large inductance magnets because the current
Table 3. Dump resistor verification tests. oscillations are induced between sections of the magnet, and
Yoke Thermal the relevant inductance for inducing these currents is much
saturation IFCL coupling Quench lower than total magnet inductance.
In this comparison, the dump resistor in the main circuit
No Loss x
is set to zero, and a CLIQ resistor of 15 mΩ and capacitor of
IFCL Only x x x
IFCL w/Quench x x x x
35 mF are added. For the results shown, the capacitor was
charged to 350 V. The CLIQ system is activated between

6
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Figure 6. The induced current oscillation in the two coil sections

about the starting current of 13.8 kA from a CLIQ discharge is
shown for a no loss case. The oscillation period is similar to the
IFCL time constant in the high field region.

Figure 4. ANSYS and COMSOL results are compared for a dump

resistor extraction exhibiting strongly coupled electromagnetic and
thermal behavior. Effects are built up from no losses, to only IFCL,
and finally to a fully coupled case where quench and IFCL are
considered. A slight difference in the initial current decay between
the ‘No Loss’ and other two cases is seen up to about 20 ms. Here
the coil is not yet quenched, but the influence of coupling currents on
the magnet’s differential inductance causes the current to decay
faster. At 20 ms, quench resistance growth begins to play a dominant
role as ‘IFCL w/Quench’ quickly decays away from ‘IFCL Only’.
This is an example of quench back, where the rapid field change
initiated by the dump resistor induces IFCL which heats the coil to
quench. The resistance rise due to quench drives the current down
much faster than a case where the effects of IFCL and quench are not
considered.

Figure 7. The currents in the different coil sections following a CLIQ

discharge are compared between the ANSYS and COMSOL models.
The induced field oscillations generate IFCL heating which leads to
quenching of a large portion of the coil (also see figure 8).

Figure 5. The rise in coil resistance and peak coil temperature is

compared between ANSYS and COMSOL for the final, fully
coupled simulation (IFCL w/Quench).

Figure 8. The rise in coil resistance and peak coil temperature is

Table 4. Energy loss comparison for IFCL w/Quench case. compared between ANSYS and COMSOL following the CLIQ
discharge.
Location ANSYS COMSOL
E (kJ) % E (kJ) % 5.00 and 5.01 ms with the magnet previously set up in a static
Dump Res. 505.96 34.45 497.25 33.80 condition at 4.5 K and 13.8 kA (90% of short-sample). The
IFCL 14.32 0.98 14.24 0.97 CLIQ discharge induces current oscillation between the
Coil Res. 948.24 64.56 959.55 65.23 magnet layers on a similar scale of the IFCL time constant
Note. % energy deposition is based on a total energy (see the no loss case in figure 6). Figures 7 and 8 show how
change between 5 and 500 ms. the IFCL heating induced by this oscillation drives the coil to

7
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Table 5. Overview of the Nb3Sn SCU short model. effects from 2D to 3D, were chosen to match the physical
length of the coil and the linked flux of a 3D ANSYS model.
Parameter Value Unit
The use of a large scaling factor feff between the matrix
Period length 19 mm resistivity and an effective resistivity for all coupling losses is
Pocket width 6.32 mm based on measurements made on a similar Nb3Sn strand
Pocket height 4.67 mm [29, 30]. For each case, a single static solution was first solved
Powered turns per pocket (left to right 7-35-56-56 with the voltage source set to produce the initial current
in sym.)
matching the test. This source was then ramped to zero over
Strand architecture 132/169
RRP®
0.1 ms to effectively place the magnet in series with a fixed
Strand diameter 0.6 mm resistance dump resistor for current extraction. The solution
Filament twist pitch 12 mm proceeds with transient effects using the Multi-field Solver for
Cu RRR 250 coupling of independently meshed electromagnetic and ther-
Non-Cu fraction 0.45 mal domains as described in section 3.
feff with ρet=feff ρst 4.0 To match data taken during the test, current extractions
Effective resistive coil length (Lc) 97.1 mm with no initial quench were simulated from 400 to 800 A
Effective inductive coil length (Li) 90 mm using a dump resistor of 48.1 mΩ. Figure 10 shows the cur-
Nb3Sn Jc (4.5 K, 10 T) 2880 A mm−2 rent decay curves from these simulations including cases with
Short-sample current (4.5 K) 965 A and without IFCL and quench effects. When normalized to
Short-sample cond. field (4.5 K) 5.2 T
peak current, the ‘No Loss’ cases show only a small variation
with initial current which is the result of nonlinear magneti-
zation of the iron core with field level. At the 400 A level, the
quench, with the resulting coil resistance growth bringing ‘IFCL w/Quench’ case remains superconducting due to IFCL
down the magnet current. Comparison between the ANSYS heating not being able to overcome the large margin to
and COMSOL results in these figures show excellent agree- quench. This results in a current decay which shows little
ment for this fully coupled electromagnetic, thermal, and deviation from the no loss case. As the initial current
circuit simulation. increases up to 800 A, the margin is reduced and IFCL
heating grows. At a certain level this begins to induce quench,
adding resistance to the circuit and driving the current down
5. Comparison to Nb3Sn undulator test data faster (quench back). As expected, the degree of quench back
is seen increasing with initial current (this is particularly clear
A series of Nb3Sn superconducting undulators (SCUs) were when comparing the normalized current decay).
built and tested at Lawrence Berkeley National Laboratory as A sensitivity study of the ANSYS results to input para-
part of a R&D program for LCLS-II and future free electron meters and material property fits was performed to allow for a
lasers [28]. These magnets consist of two ferromagnetic iron more detailed comparison to test data. Table 6 summarizes
cores with machined pockets along the length. The pockets this study and the results. The deviation from the nominal
are wound with superconducting strand such that the current case was evaluated on a parameter by parameter basis by
polarity changes from one pocket to the next along the length. comparing both the quench integral and peak coil hotspot
The two cores are assembled together with a gap left between temperature for each initial current level. The quench integral
them for an electron beam. The interaction of the beam with is the time integral of the square of magnet current from the
the alternating fields generated by the SCU serves as the start to end of the decay. This is a material property inde-
radiation source in free electron lasers. pendent measure of the total energy deposited during the
These magnets are well suited for benchmarking the user decay which would generate joule heating at a quench loca-
elements due to: (1) the use of single strands eliminating cable tion. The quench integral is typically used for an estimation of
based coupling currents found in many other Nb3Sn magnets quench location hotspot temperature in the adiabatic limit
and (2) exhibiting strong quench back behavior changing as a with the material properties considered [19], and in this case
function of field level. The test of a short prototype SCU was serves as a metric for the degree of quench back when
selected for a first comparison due to existing data for a series compared to a no loss case. In the future, a more detailed
of pre-quench, dump resistor extractions at increasing levels metric could be used along with an exploration of combined
of initial magnet current. This allows for benchmarking of the effects due to multi-parameter deviation from the nom-
user elements with test data over a wide range quench back inal case.
behavior. An overview of the properties of this magnet is The results show the most sensitive fits and parameters
found in table 5 and a picture of the cross section can be seen are those associated with the resistivity of the copper matrix
in figure 9. This prototype corresponds to a short, 80 cm, ( feff, Cu RRR, and Cu resistivity fit). This is not unexpected,
length single core, whereas the final SCU magnet includes as these impact the IFCL time constant and coil resistance. A
two cores of 1.4 m length assembled together. range of behavior about the nominal case due to changes in
The ANSYS model was matched to the test configuration Cu resistivity was created to visualize this sensitivity when
of a single, short prototype magnet as seen in figure 9. The comparing to test data. This range is bounded by two curves
effective lengths Lc and Li, which scale resistive and inductive which represent the extremes of the material property fits and

8
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Figure 9. (a) The electromagnetic ANSYS model is shown in a dump resistor circuit. The meshed regions consist of ferromagnetic low
carbon steel (light blue), un-powered corrector turns (magenta), powered turns (teal), glass filled epoxy (orange), and air which extends
beyond the top and left boarders for far field boundary conditions (dark blue). This is coupled to a similar thermal region (without air) using
the Multi-field Solver to simulate IFCL induced quench back in an external dump resistor protection circuit. (b) A cross-section of the short
prototype Nb3Sn undulator with the symmetric region marked.

Figure 10. Magnet current decay curves from ANSYS are compared with and without the effects of IFCL and quench for absolute (a) and
normalized current (b). The impact of quench back is seen increasing with initial current level. The large thermal margin and smaller rate of
field change (driving IFCL) at the 400 A current level results in the magnet staying superconducting and showing minimal deviation from the
no loss case. As the initial current increases up to 800 A, IFCL induced quench plays a larger role, with coil resistance growth driving the
current down more quickly than the no loss case.

Table 6. Sensitivity study of ANSYS results for the SCU.

Variation of resultsa: ΔQint (%), ΔThot (K)

Parameter/Fitb Nominal Changed 400 A 500 A 600 A 700 A 800 A
feff 4.0 2.0 −4.6, 3.0 −6.4, 2.1 −4.7, 0.7 −3.0, 0.1 −1.1, 0.0
6.0 2.3, −1.3 4.5, −5.9 5.6, −1.3 3.8, −0.5 2.4, −0.2
Cu RRR 250 150 1.9, −1.0 3.2, −0.8 0.4, 1.3 −2.3, 1.7 −4.2, 1.8
350 −1.1, 0.4 −1.2, 0.0 0.0, −0.6 1.1, −0.8 2.1, −0.9
Cu resistivity NIST CUDI 0.9, −0.5 1.9, −0.8 0.4, 0.7 −1.6, 1.4 −3.7, 2.2
MATPRO −0.1, 0.0 0.0, −0.2 0.5, −0.6 1.4, −1.0 2.4, −1.4
Cu heat capacity NIST CUDI 0.0, 0.0 0.0, 0.0 0.0, −0.1 0.0, −0.1 0.1, −0.2
Cu thermal cond. NIST CUDI 0.0, 0.0 0.0, 0.0 0.0, 0.0 0.0, 0.0 0.0, 0.0
Nb3Sn heat capacity NIST CUDI 0.0, −0.4 1.3, −1.0 1.9, 0.0 1.8, 0.5 1.7, 1.4
G10 heat capacityc NIST Fermilab 0.0, 0.6 −2.6, 1.5 −3.4, 0.8 −3.7, 0.4 −3.8, 0.0
G10⊥ therm. cond. CryoComp NIST 0.0, −0.4 0.9, −1.1 0.9, −0.9 0.9, −0.7 1.0, −0.7
a
Based on deviation of the results from the nominal case for the variation of the single listed parameter. ΔQint is the change in
quench integral and ΔThot is the change in peak hotspot temperature of the coil;
b
All material property fits can be found in [31];
c
G10 properties (normal) are used for the epoxy impregnated glass fiber region between strands.

9
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

Figure 11. Comparison of ANSYS results with test data shows both exhibiting similar quench back behavior over a range of initial currents.
Deviation from the no loss cases is minimal at low current and increases with current level as quench back becomes more prominent. To
visualize the sensitivity of the ANSYS model to the parameters and material fits affecting Cu resistivity, a range is plotted around the nominal
case corresponding to the limiting behavior found in the sensitivity study (as described in table 6 and section 5).

RRR. One curve corresponds to a RRR of 150 and the fit meshing, solving, and post-processing capabilities of the
from CUDI, and the other a RRR of 350 and the fit from standard distribution. These new elements were shown
MATPRO. Figure 11 compares this range and the nominal benchmarked against existing codes for a Nb3Sn dipole and
case to test data. Current decay curves are shown from 400 to compared to test data for a Nb3Sn prototype undulator. In
800 A, along with the quench integral of these decays from 5 both cases the user elements were shown predicting IFCL
to 50 ms. The ANSYS simulations are seen reproducing the induced quench back, demonstrating for the first time that
trend seen in the measured data of larger deviation from the ANSYS can be used to simulate this strongly coupled beha-
no loss case at higher current. The source of the remaining vior required for accurate modeling of many superconducting
difference between the ANSYS predictions and test data magnets.
could originate from the 2D elements not including long- This work is part of a larger effort within the US Magnet
itudinal quench propagation or 3D effects on peak field, Development Program and the Berkeley Center for Magnet
inductance, and coupling loss. In addition, further study of the Technology to advance analysis and modeling capabilities for
accuracy of equivalent magnetization models for coupling superconducting magnets [32]. It is our goal that this work
currents in this application may be revealing. Despite room becomes a tool usable by the magnet design community. The
for further study, the level of agreement between ANSYS and effort to make these elements available is underway, and
test data is a promising sign the user elements can be used to interested parties are encouraged to contact the authors for
understand and predict quench back effects in Nb3Sn SCU’s. more information. Future work is focused on the extension of
this approach to 3D and towards simulation of HTS coated
conductors and bulk superconducting devices by imple-
menting the E–J power law model within the A–V
6. Conclusion formulation.

User defined elements in ANSYS allow for programming the

generation of the finite element matrices, providing the
opportunity to implement new formulations and material Acknowledgments
property fits. Two user defined elements were created as a
new method for simulating coupled, multiphysics behavior of Lucas Brouwer would like to thank Emmanuele Ravaioli of
superconducting magnets. These elements extend the cap- CERN (previously LBNL) for many interesting and helpful
abilities of ANSYS to now include superconducting specific discussions about the simulation of dynamic effects in
phenomenon of quench and IFCL, while maintaining the superconducting magnets and CLIQ based magnet protection.

10
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

This work was supported by the Director, Office of Science, The user element is coupled to an external circuit as
High Energy Physics, and US Department of Energy under voltage source using the standard distribution, circuit element
contract No. DE-AC02-05CH11231 in the context of broader CIRCU124 [22]. With key option one set to select a stranded
collaboration with the US Magnet Development Program. coil, this element consists of three nodes labeled i, j, and k.
The first two nodes are connected to adjacent circuit elements
and each carry a single voltage DOF. The third node carries
Appendix both a current and voltage drop DOF, and is chosen as one of
the nodes in the meshed coil region to make it part of the
The matrices used for the 2D electromagnetic element with coupled set. The stiffness matrix for the stranded coil
circuit coupling found in equations (16) and (17) are the result CIRCU124 element is given by
the following steps. The K AA and C AA matrices in 2D are
found by simply reducing the general form already given in ⎡0 0 1 0 ⎤ ⎡ Vi ⎤ ⎡ 0 ⎤
⎢0 ⎢ ⎥
equations (8) and (9). For the stranded conductor, the source 0 -1 0 ⎥ ⎢V j ⎥ = ⎢ 0 ⎥ ,
⎢ ⎥ ⎢ ⎥ (A.8)
term is now supplied by the circuit and determines the form of ⎢- 1 1 0 s ⎥ ⎢ ik ⎥ ⎢ 0 ⎥
K Ai. Considering the weak integral of the source term in 2D ⎣0 0 0 0 ⎦ ⎢⎣ ek ⎥⎦ ⎣ 0 ⎦
( J is now a scalar Jz), and that the current density of the
stranded conductor is derived from the current per turn i using where s is a factor to account for modeling of a symmetric
N
Jz = Sc ti leads to region. This couples the stranded coil into an external circuit
c which may be made up of additional coil regions or generic
Nc circuit elements selected using other key options for
ò {N} Jz dS = ò {N} Sc t {N} T {i} dS, (A.1)
CIRCU124.

where {N} is a vector of shape functions carrying the DOF

across the element based on the values at the element nodes
ORCID iDs
(i={N}T{ie} for example). Here clearly the K Ai matrix is
given by
Lucas Brouwer https://orcid.org/0000-0003-2170-7278
tN
[K Ai ] =- c
Sc ò {N}{N} T dSelem. (A.2)
References
The 2D matrices in the voltage balance are derived from
equations (11), (13), and (15). The e term is re-written as an [1] ROXIE: a Computer Code for the Integrated Design of
area integral with the shape functions added, such that Accelerator Magnets URL https://cern.ch/roxie
[2] Guo X L, Wang L and Green M A 2012 Cryogenics 52 420
Nc ¶ 1 [3] Marinozzi V, Sorbi M, Manfreda G, Bellina F, Bajas H and
t
Sc
Li ò {N} T ¶t {Az} dS - Sc ò {N} T {e} dS Chlachidze G 2015 Phys. Rev. Spec. Top. Accel. Beams 18
032401
⎛ N ⎞2 I fcu [4] Breschi M, Cavallucci L, Ribani P L, Gavrilin A V and
+ Lc ⎜ c ⎟ r
⎝ Sc ⎠ fcond (1 - fsc ) st ò {N} T {i} dS = 0. (A.3)
Weijers H W 2017 IEEE Trans. Appl. Supercond. 27
4301013
ANSYS assembles the three DOF using the following form [5] Ravaioli E, Auchmann B, Maciejewski M, ten Kate H H J and
Verweij A P 2016 Cryogenics 80 346
⎡ ⎤ [6] Bortot L, Auchmann B, Cortez-Garcia I,
⎡C AA 0 0 ⎤ ⎢ ¶A ⎥ ⎡ K AA 0 K Ai ⎤ ⎡ A ⎤ ⎡ 0 ⎤ Fernandez-Navarro A M, Maciejewski M, Prioli M,
⎢ eA ⎥ ⎢ ¶t ⎥ + ⎢ ⎥ ⎢ ⎥
0 K ee K ei ⎥ ⎢ e ⎥ = ⎢ 0 ⎥. (A.4)
Schops S and Verweij A 2018 IEEE Trans. Appl.
⎢C 0 0⎥
⎢ 0 ⎥ ⎢ ⎢
⎣ ⎥
⎦ Supercond. 54 7000404
⎣ 0 0 0⎦ ⎣ 0 0 0 ⎦ i ⎣0⎦
⎣ 0 ⎦ [7] Schops S, Gersem H D and Bartel A 2010 IEEE Trans. Magn.
46 3233
To match the sizing of the submatrices here, the vectors in [8] Bortot L, Auchmann B, Cortez-Garcia I,
Fernandez-Navarro A M, Maciejewski M, Mentink M,
equation (A.3) are expanded to square matrices using the
Prioli M, Ravaioli E, Schops S and Verweij A 2018 IEEE
outer product with an identity vector {I}. This leads to 2D Trans. Appl. Supercond. 28 4900706
element matrices of [9] ANSYS https://ansys.com
[10] Caspi S, Chiesa L, Ferracin P, Gourlay S A, Hafalia R,
Nc
[C eA] = t
Sc
Li ò {I}{N} T dSelem, (A.5) Hinkins R, Lietzke A F and Prestemon S 2003 IEEE Trans.
Appl. Supercond. 13 1714
[11] Yamada R, Marscin E, Lee A and Wake M 2004 IEEE Trans.
1 Appl. Supercond. 14 291
[K ee] = -
Sc ò {I}{N} T dSelem, (A.6) [12] Ferradas-Troitino J et al 2019 IEEE Trans. Appl. Supercond.
29 4701306
[13] Brouwer L, Arbelaez D, Caspi S, Marchevsky M and
⎛ N ⎞2 I fcu Prestemon S 2018 IEEE Trans. Appl. Supercond. 28
[K ei ] = L c ⎜ c ⎟ r
⎝ Sc ⎠ fcond (1 - fsc ) st ò {I}{N} T dSelem. (A.7) 4001006
[14] Hauser A 1997 IEEE Trans. Magn. 33 1572

11
Supercond. Sci. Technol. 32 (2019) 095011 L Brouwer et al

[15] Testoni P 2003 Implemenation in the ANSYS finite element [26] Brouwer L, Auchmann B, Bortot L and Stubberud E 2019
code of the electric vector potential T-Ω, Ω formulation and Crosscheck of the ANSYS-COMSOL 2D FEM
its validation with the magnetic vector potential A-V, A Implementations for Superconducting Accelerator Magnets
formulation PhD Thesis University of Cagliari, Italy SU-1010-4841,R1 https://usmdp.lbl.gov/scpack-code/
[16] COMSOL Multiphysics https://comsol.com Lawrence Berkeley National Laboratory
[17] 2016 ANSYS Mechanical APDL Theory Reference [27] Ravaioli E 2015 CLIQ. A new quench protection technology
Release 17.1 for superconducting magnets PhD Thesis University of
[18] Biro O and Preis K 1989 IEEE Trans. Magn. 25 3145 Twente, Netherlands
[19] Wilson M 1983 Superconducting Magnets (Oxford: Oxford [28] Emma P et al 2014 A plan for the development of
University Press) superconducting undulator prototypes for LCLS-II and future
[20] Morgan G 1970 J. Appl. Phys. 41 3673 FELs FEL 2014 Conf. Proc. (Basel, Switzerland) p THA0
[21] Louzguiti L, Zani L, Ciazynski D, Turck B and Topin F 2016 [29] Zhou C, van Lanen E, Veldhuis D, ten Kate H, Dhalle M and
IEEE Trans. Appl. Supercond. 26 4700905 Nihuis A 2011 IEEE Trans. Appl. Supercond. 21 2501–4
[22] Wang J S 1996 IEEE Trans. Magn. 32 1071 [30] Zhou C, Miyoshi Y, van Lanen E, Dhalle M and Nihuis A
[23] Russenschuck S 2010 Field Computation for Accelerator 2012 Supercond. Sci. and Technol. 25 065018
Magnets: Analytical and Numerical Methods for [31] Manfreda G 2011 Review of ROXIE’s Material Properties
Electromagnetic Design and Optimization (Wiley-VCH) Database for Quench Simulation 1178007 CERN
[24] Stekly Z J J and Zar J L 1965 IEEE Trans. Nucl. Sci. 12 367 [32] Gourlay S, Prestemon S, Zlobin A, Cooley L and
[25] STEAM: Simulation of Transient Effects in Accelerator Larbalestier D 2016 The US Magnet Development Program
Magnets https://espace.cern.ch/steam/ Plan https://osti.gov/scitech/servlets/purl/1306334

12