Optimized Field/Circuit Coupling for the Simulation of Quenches in
Superconducting Magnets
I. Cortes Garcia1 , S. Schöps1 , M. Maciejewski2,3 , L. Bortot2 , M. Prioli2 , B. Auchmann2,4 , and A.P. Verweij2
1 Technische

arXiv:1702.00958v2 [physics.comp-ph] 6 Jul 2017

Universität Darmstadt, Darmstadt, Germany
2 CERN, Geneva, Switzerland
3 Łódź University of Technology, Łódź, Poland
4 Paul Scherrer Institut, Villigen, Switzerland

In this paper, we propose an optimized field/circuit coupling approach for the simulation of magnetothermal transients in
superconducting magnets. The approach improves the convergence of the iterative coupling scheme between a magnetothermal
partial differential model and an electrical lumped-element circuit. Such a multi-physics, multi-rate and multi-scale problem requires
a consistent formulation and a dedicated framework to tackle the challenging transient effects occurring at both circuit and magnet
level during normal operation and in case of faults. We derive an equivalent magnet model at the circuit side for the linear and the
non-linear settings and discuss the convergence of the overall scheme in the framework of optimized Schwarz methods. The efficiency
of the developed approach is illustrated by a numerical example of an accelerator dipole magnet with accompanying protection
system.
Index Terms—Convergence of numerical methods, coupling circuits, eddy currents, iterative methods.

I. I NTRODUCTION
UPERCONDUCTING magnets produce high magnetic
fields used in high-energy particle accelerators for bending
particle beams. In order to reach the superconducting state, the
magnets are operated at very low temperatures (1.9 K). Since
the heat capacity is low at cryogenic temperatures, the magnets
are prone to quench due to a local energy deposition (coupling
losses in the superconducting cable, beam losses, cryogenic
malfunction, mechanical movement, etc.). A quench is a transition from the superconducting to the normal conducting state.
As a consequence, the release of the magnetic energy as Ohmic
losses may result in a catastrophic damage in the magnet and
circuit.
The simulation of quench initiation, propagation, and subsequent protective measures represents a challenge in terms of
the number of coupled physical domains, their highly nonlinear
behavior, their geometric scales, and their vastly different
time constants. Quench protection systems such as the quench
heaters [1] and the coupling-loss induced quench system [2]
are affecting both magnet and circuit. Therefore, their mutual
influence has to be studied carefully. To this end, field-circuit
coupling is inevitable. It brings together a magneto-thermal
problem based on partial differential equations and a circuit
model of the magnet’s powering circuit. In the paper, we are
presenting an equivalent superconducting magnet model on
the circuit side capturing the magnetic and thermal behavior
needed to represent a quench in the superconducting magnet.
The coupling of electric fields and circuits is a well understood technique. Differences among various approaches are
how the coupled problem is tackled. It is either solved ‘monolithically’ as one system, such as e.g. proposed in [3] or it is
coupled ‘weakly’ using additional iterations or exploiting temporary linearizations [4], [5], [6]. In most cases, the coupling is
established via current sources, voltage sources or inductances.

Iron yoke

S

B

F

v

Fig. 1. Coil surrounded by an iron yoke. The current flowing through the
coil into the plane is depicted with red and its return path is marked with
blue. Positively charged particle traveling in the middle of the coil’s aperture
subject to vertical component of the magnetic field will be deflected due to
the Lorentz force (green dashed line).

In [7], [8] the various approaches have been interpreted in
the context of waveform relaxation methods [9], [10] such
that convergence could be proven. This paper interprets the
coupling conditions in terms of optimized Schwarz methods,
where they are often called ‘transmission conditions’, e.g. [11],
[12]. In contrast to prior works, o numerical optimization is
carried out, since it is cumbersome in the case of multi-port
devices and it is limited to linear cases. Instead, it is shown
that inductances are a first order approximation to the optimal
condition.
The paper is structured as follows: Section 2 deals with
the modeling of superconducting magnets and their spatial
discretization, Section 3 introduces the idea of waveform
relaxation and discusses the optimized contraction factor. Section 4 underlines the findings by a real world simulation of
an accelerator magnet coupled to its surrounding protection
circuitry. The paper closes with conclusions and an outlook.

c 2017 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media,
including reprinting/republishing this material for advertising or promotional purposes, creating new collective works, for resale or redistribution to
servers or lists, or reuse of any copyrighted component of this work in other works.

II. M ODELING
Starting from Maxwell’s equations and assuming a magnetoquasistatic (MQS) setting, i.e., neglecting displacement
currents, the following partial differential equation
!
~
∂
A
~ + ∇ × νs τeq ∇ ×
∇ × (ν∇ × A)
=χ
~ si
(1)
∂t
~ 0 at time
is obtained on a domain Ω ∈ R with an initial value A
t0 and with Dirichlet and Neumann boundary conditions on
Γdir ∪Γneu = Γ := ∂Ω, respectively. The eddy currents are taken
into account in the coil domain Ωs by homogenization [2],
~ is the magnetic vector potential, ν and νs
[13], see Fig. 2. A
are inverse magnetic permeabilities that may depend on B =
~ with B
~ = ∇×A
~ the magnetic flux density and i the
kBk
lumped currents through each coil. In Ωs , the time constant
~ s satisfies
τeq is defined such that the cable magnetization M
the relation
~
~ s = τeq ∂ B ,
νs−1 M
∂t

Γneu
Γs
Ω

3

where τeq also depends on B. Finally, χ
~ s : Ω 7→ R3×ns , with ns
the number of coils, is the stranded-conductor winding function
(see [14]), which homogeneously distributes the current in the
domain Ωs such that the current density J~s = χ
~ s i. The field
equation (1) is coupled to a circuit equation model via the flux
linkage
Z
d
~
v = Ψ, with Ψ =
χ
~>
(2)
s A dΩ.
dt
Ω
The temperature in the domain Ωs can be obtained from the
heat balance equation
ρCp

∂T
− ∇ · (k∇T ) = Ps + PJoule ,
∂t

(3)

with suitable boundary and initial conditions, where ρ is the
mass density, Cp the heat capacity, k the thermal conductivity
and T the temperature. The Joule losses PJoule are defined as
PJoule = qflag σ −1 kJ~s k2
with σ being the homogenized non-linear conductivity of
the cables in the case of a quench. It is activated by the
sigmoid-type activation function qflag that depends on time, the
~ and current density J~s . Ps is the power
magnetic flux density B
density in the superconductor and relates the heat equation with
the magnetoquasistatic field solution through
~
~ s · ∂B .
Ps = M
∂t
In case of quenching, i.e., qflag > 0, the superconducting coils
have an Ohmic resistance with the voltage drop
Z
vt = Ri with R = qflag
χ
~ s> σ −1 χ
~ s dΩ.
(4)
Ωs

The non-linear resistance R inherits the dependencies on t,
~ =∇×A
~ and J~s = χ
B
~ s i.

Ωs
Γdir
Fig. 2. Sketch of the computational domain and subdomains with boundaries

After spatial discretization of (1-4) as e.g. obtained by the
Finite Element (FE) Method [15], the semi-discrete system

 da
+ Kν a a − Xi = 0
(5a)
Ks a
dt
da
v − X>
=0
(5b)
dt
dt
Mρ
+ Kk t − q(a) = 0
(5c)
dt
vt − R(t, a, i)i = 0
(5d)

is obtained, where a denotes the degrees of freedom for
the magnetic vector potentials, t the discrete temperatures,
q(a) the discretized Joule losses, M? and K? the discretized
material and differential operator matrices and X is the discretization of χ
~ s . Initial conditions v0 , vt,0 , a0 and t0 at time
t0 are set. The system is differential algebraic, as Ks is a
singular matrix due to the fact that τeq is only non-zero in the
coil region (and additionally in the 3D case as a consequence
of the null space of the curl-curl operator).
The equations describing the overall behavior of the circuit
can be written in an abstract form as the following system

dx
A
+ B x x + Pi = f (t)
(6a)
dt
P> x − v = 0
(6b)
Q> x − vt = 0.

(6c)

with initial conditions x(t0 ) = x0 and i(t0 ) = i0 . In the case
of modified nodal analysis (see [16]), which is used in most
modern circuit simulators, the degrees of freedom x contain
the electric potentials at the circuit nodes and currents through
the circuit branches associated with voltage sources, inductors
and magnets (coupled via the voltages v and vt ).
III. WAVEFORM R ELAXATION
Following [7], waveform relaxation is used to simulate the
coupled problem (5-6). In this approach, the magnetothermal
problem (5) is simulated separately from the circuit (6). Information between the two systems is exchanged only at certain
synchronization points t̄0 , . . . , t̄n . Among other advantages,
this exploits the different time scales at which significant

effects occur in each subproblem, as the time stepper of each
problem can choose its step sizes independently. A GaussSeidel-type scheme on a time window Ij = [t̄j , t̄j+1 ] using
current and voltages as ‘transmission conditions’ is given by
System 1 (circuit, cf. (6))
Ac

d (k+1)
x
+ Bc x(k+1) + Pi(k+1)
= f(t)
c
dt
(k+1)
vc
= P> x(k+1)
(k+1)
vt,c

= Q> x(k+1)

(7a)
(7b)
(7c)

A. Optimized waveform relaxation
The following analysis is carried under several simplifying
assumptions: only the coupling between electromagnetic field
(8a-8b) and circuit (7a-7b) is considered, the thermal case can
be treated analogous. Moreover, the reluctivities ν and νs , as
well as the time constant τeq are kept constant. It is shown
later that these steps do not limit the applicability in practice.
Let us rewrite the system in frequency domain using the
angular frequency ω = 2πf . Then
Ac jωx(k+1) + Bc x(k+1) + Pi(k+1)
= g(ω)
c

with the transmission conditions
(k)
vc(k+1) = vm
(k+1)

vt,c

(7d)

(k)

(7e)

= vt,m

d (k+1)
a
+ Kν a(k+1)
dt
d
X> a(k+1)
dt
d (k+1)
Mρ t
+ Kk t(k+1)
dt
(k+1)
vt,m

= Xi(k+1)
m

(8a)

(k+1)
vm



(8b)

= q a(k+1)



= R(k+1) i(k+1)
m

(8c)
(8d)

with the transmission condition
i(k+1)
= i(k+1)
.
m
c

(8e)

The algorithm consists of the following steps
(0)

(0)

0) set vm (t) = v0 and vt,m (t) = vt,0 on t ∈ I0
(k)

=P x
=

(k)
vm

(k)

1) solve (7) for given inputs vm (t) and vt,m (t) on Ij to
(k+1)
(t)
obtain the output ic
(k+1)

(t) on Ij to obtain the
2) solve (8) for given input ic
(k+1)
(k+1)
outputs vm (t) and vt,m (t)
3) if not converged, set k = k+1 and go to step 1), otherwise
(k+1)
(0)
(k+1)
(0)
proceed with vm (t) = vm (t̄j ), vt,m (t) = vt,m (t̄j ),
for t ∈ Ij+1 , j = j + 1 and go to step 1).
Here, initialization of the voltages at each waveform relaxation
window is carried out via constant extrapolation of its value
at the last time step of the previous window. However, this
initialization can also be implemented differently, e.g. by using
higher order extrapolation.
It shall be noted that the excitation of the field problem (8) by
a known voltage source instead of a current source increases
the numerical stability, as the differential algebraic index of
the system is minimized [17]. However, for this paper this
aspect is disregarded and the field problem is considered in the
circuit as a current-driven voltage source. This corresponds,
from a Schwarz-type methods point of view, to using ‘zeroorder’ Dirichlet transmission conditions. As in [11] and [12],
the waveform relaxation scheme can be studied from the
point of view of an optimized Schwarz method, such that
the transmission conditions are optimized to speed up the
convergence in terms of k.

(9a)
(9b)
(9c)

and
>

=

> (k+1)

Ks jωa(k+1) + Kν a(k+1) = Xi(k+1)
m

and System 2 (magnetothermal, cf. (5))
Ks

vc(k+1)
vc(k+1)

(k+1)

=

i(k+1)
m

=

X jωa

(k+1)
vm
i(k+1)
.
c

(10a)
(10b)
(10c)

For the following optimization of convergence, the first transmission condition (9c) is generalized to the linear combination
(k)
vc(k+1) = αi(k+1)
− αi(k)
c
m + vm ,

(11)

while the second transmission condition (10c) is kept in its
original form since its optimization is cumbersome for arbitrary
circuitry.
The following section discusses how the weighting factor α
influences the convergence of the waveform relaxation scheme.
B. Optimization of contraction factor in frequency domain
Waveform relaxation is a fixed point iteration in time domain
and therefore its convergence can only be guaranteed if there
is a contraction of the error, see e.g. [10]. This is expressed by
a contraction factor ρ(α) < 1, such that
||vc(k+1) − vc(k) || = |ρ(α)| ||vc(k) − vc(k−1) ||.

Assuming that (Ks jω + Kν ) is invertible, which is necessary
for the solvability of the subsystem, we obtain
−1

ρ(α) = I + αx−1
α − Z(ω) x−1
P (ω)
P (ω),

where we assume that the transfer impedance xP (ω) =
P> (Ac jω+Bc )−1 P exists and is invertible and the impedance
of the magnetoquasistatic system
Z(ω) = jωX> (Ks jω + Kν )−1 X.

(12)

Optimal convergence is attained for ρ(α) = 0, that is
α = Z(ω).
Therefore, the optimized transmission condition for the circuit
is given by
(k)
vc(k+1) = Z(ω)i(k+1)
− Z(ω)i(k)
c
m + vm .

In frequency domain, the impedance Z is computable and
no further (numerical) optimization is needed to obtain an
optimal transmission condition. However, when proceeding to
non-linear systems in time domain, Z contains time derivatives
which need to be approximated.

C. Approximation of the optimal condition by inductances
Let us assume that Kν is invertible. Then, we can rewrite
the impedance as


−1
−1
− 1 −1
−1
Kν 2 X,
Z(ω) = jωX> Kν 2 I + jωKν 2 Ks Kν 2
which can be recast in a Neumann series
− 12

Z(ω) = jωX> Kν

∞ 

X
−1
−1 l
−1
−jωKν 2 Ks Kν 2 Kν 2 X,
l=0

if the operating frequency ω is low, as it must hold
−1

−1

κ := || − jωKν 2 Ks Kν 2 || < 1.
To analyze the smallness of κ, we make some further assumptions: let the magnetoquasistatic system be discretized,
such that the curl-curl operators can be factorized as K? =
C> M? C with diagonal material matrices M? . Then, it can be
shown that
1
fmax <
2πτmax
holds, since
−1

−1

κ ≤ ωmax ||Kν 2 C> Mτmax Mν CKν 2 ||,
with Mτmax = max(Mτ )ii I = Mτ + M̄τ , and M̄τ a

D. Improved algorithm in time domain
When treating non-linear field/circuit coupled problems with
waveform relaxation, the inductance in the transmission condition (14) is time dependent due to the magnetic saturation in
(5), i.e.


a(k) (t) X.
L(k) (t) := X> K−1
ν
Following this approach rigorously requires to continuously
update the inductance, e.g. by repeating the inductance extraction described in the previous section after each successful
time step. In contrast, we propose a simplified procedure: the
differential inductance
Lm = X> K−1
ν,d (a? )X

is extracted at a working
 point a? from the differential curld
curl matrix Kν,d = da Kν a a and it is kept constant for
the subsequent waveform relaxation, i.e., the condition (9c) is
replaced by
vc(k+1) (t) = Lm

Finally, the optimal coefficient (12) can be approximated by a
finite number of terms of the Neumann series, e.g. for l = 0,
it follows that
Z(ω) ≈ jωL := jωX> K−1
ν X.

(13)

Using the coefficient α = jωL and dividing the transmission
condition (11) by jω yields
(k)
Ψ(k+1)
= Li(k+1)
− Li(k)
c
c
m + Ψm ,

(14)

1
v? . The inductance
with the magnetic flux linkage Ψ? = jω
can be obtained by solving the magnetostatic problem, i.e.
(1), with τ = 0 excited by 1 A in each coil. This extraction
comes with moderate additional numerical costs, since only
one (linear) system of equations has to be solved with several
right hand sides. The transmission condition (14) corresponds
to considering the field model in the circuit as an inductance
with a correction term as proposed in [18]. Considering further
terms of the Neumann series is possible but has not been investigated. This would require higher regularity of the solution
since higher order time derivatives are needed.
Although the analysis above has been carried out for a linearization of the original system and under restrictive assumptions on the discretization, in practice, it has been observed
that the convergence of the waveform relaxation scheme with
the transmission condition (14) is substantially improved. The
following section discusses the non-linear case.

d (k+1)
(k)
i
(t) + ∆vm
(t),
dt c

(16)

which corrects the misfit due to disregarding the non-linearity
and the higher order terms of the Neumann series by a lumped
voltage source

i

positive semidefinite diagonal matrix. Mν = Mνs + M̄νs , with
M̄νs also positive semidefinite, as (Mνs )ii = (Mν )ii in the
coil domain and zero everywhere else. Therefore the series
converges if
ωmax τmax < 1.

(15)

(k)
∆vm
(t) :=

d (k)
d
Ψ (t) − Lm i(k)
(t).
dt m
dt m

(17)

Inspired by the inductance coupling in the field-circuit case,
resistances are now passed from the thermal equation to the
circuit instead of the voltages in equation 8d through
(k+1)

vt,c

= R(k) (t)ik+1
,
c

(18)

(k)

where R(k) (t) = R(t, a(k) , im ). The Gauss-Seidel-type
scheme using these new transmission conditions reads
(0)

0) set k = 0, ∆vm (t) = ∆vm,0 (t) and
R(0) (t) = R(t0 , a0 , i0 ) and extract Lm for a? = a0
using (15)
1) solve (7a-7c) with conditions (16) and (18) for given
(k)
(k+1)
inputs ∆vm (t) and R(k) (t) to obtain the output ic
(t)
(k+1)

2) solve (8) for given input ic
(t) to obtain the output
(k)
voltage ∆vm (t) according to (16) and R(k+1) (t) as
before
3) if not converged, set k = k + 1 and go to step 1),
(0)
(k+1)
otherwise proceed with ∆vm (t) = ∆vm (t̄j ),
(0)
(k+1)
R (t) = R
(t̄j ), j = j + 1 and go to step 1).
Obviously, the algorithm can be improved by appropriate
heuristics to update Lm if the number of iterations k grows.
IV. N UMERICAL E XAMPLES
Let us first discuss the convergence of the Neumann series in
practice and then give an example for the field/circuit quench
simulation of the D1 accelerator magnet [19].

kL
kL
kL
kL
kL

+
−

10

∆Um (t) Rm (t)

kL Lm

REE = 0.1 Ω
Fig. 3. Schematic of the protection circuit with a lumped model representing
the magnet (dashed box), where Rm (t) is defined in equation (21), ∆vm (t)
in (17) and Lm in (19). Finally, im (t) is the current through the magnet.

Inner iteration

im (t)
8

=1
= 0.9
= 0.8
= 0.7
= 0.6

6

4

2
1

0

2

3

∆vm [V]

Outer time window
Fig. 5. Convergence comparison for first 3 time windows of the magnetostatic
case

−100
kL
kL
kL
kL
kL

−200
0

0.2

0.4

=1
= 0.9
= 0.8
= 0.7
= 0.6

0.6

0.8

t [s]
Fig. 4. Voltage differences ∆vm correcting the mismatch between field and
circuit model in the magnetostatic case

A. Discussion of the time constant
The time constant τeq depends on the cable parameters and
is mathematically expressed by [13]
µ0
τeq (B) =
2



lf
2π

2

1
,
(c0 + c1 B)feff,x

with lf being the filament twist-pitch, µ0 the vacuum permeability, feff,x the fraction of superconductor in the matrix, c0
and c1 being characterized by the resistivity of the material
used in the matrix and 0 < B ≤ 10 being the magnitude of
the magnetic flux density.
Typical values for those parameters are given in [2]. For the
case of a Nb-Ti dipole accelerator magnet, one has
lf = 1.5 · 10−2 m,
−11

c1 = 4.2 · 10

ΩmT

−1

,

c0 = 1.7 · 10−10 Ωm

feff,x = 1

and thus it follows
τeq <

µ0
2



lf
2π

2

1
< 0.0211 s,
c0 feff,x

which yields convergence of the Neumann series for frequencies below fmax < 7.5 Hz.

B. Simulation of dipole magnet
The described method is illustrated for a numerical example
of the single aperture dipole magnet D1 [19]. The magnet
is designed to be used for beam separation for the High
Luminosity upgrade of the Large Hadron Collider. The dipole
magnet is represented by the strongly coupled magnetothermal
system (5) in a 2D setting. Only one quarter is modeled
with homogeneous Dirichlet boundary conditions at the outer
boundaries and Neumann conditions at the symmetry planes

~ Γ × ~n = 0 Wb/m,
~
×~n = 0 A/m
A|
ν∇ × A
dir
Γneu
and k∇T |Γs · ~n = 0 W/m2 , where ~n is the outward normal vector. The partial differential equations (5a)-(5d) are
discretized and solved using COMSOL M ULTIPHYSICS R
[20]. The model consists of first order and second order
nodal elements for the thermal and the magnetic problem,
respectively, adding up to 9871 degrees of freedom in total.
The semi-discrete problem is time-discretized by a backward
differentiation formula (variable order, maximum time step size
δt = 1 ms). All parameters, e.g. non-linear materials laws for
the reluctivity ν and time constant τ , are defined as specified
in [19], [21]. We consider a magnet operating at 5 kA and
protected by a resistor REE = 0.1 Ω as depicted in Fig. 3. The
circuit is modeled and simulated in O R CAD PS PICE R using
the trapezoidal method (maximum time step δt = 10 µs) [22].
The waveform relaxation is carried out on I = (0 s, 0.75 s] with
~ 0 , T0 and x0 that correspond to a system ramped
initial values A
up to 5 kA. The ramped-up values are obtained by initializing
the circuit and magnetothermal system with zero conditions
and then solving both problems independently, with current
sources that are linearly increased from 0 to 5 kA. The cosimulation was established within CERN’s in-house coupling
tool STEAM (Simulation of Transient Effects in Accelerator
Magnets) [23].
The aim of this study is to analyze the influence of the
differential inductance estimation (15) on the waveform relaxation convergence. For that purpose, we introduce a scaling

coefficient kL to be multiplied by the (scalar) differential
inductance
Lm = X> K−1
(19)
ν,d (a0 )X,

Rm = qflag X> M−1
σ X

(21)

is the time-dependent coil resistance appearing after a quench
and ∆vm (t) is the (scalar) voltage source accounting for
differences between the induced voltage of a FE model and the
voltage across the magnet’s differential inductance according to
(16). In the case of the dipole, the voltage source is determined
by
d
d
ψ(t) = kL Lm im (t) + ∆vm (t).
dt
dt
This is represented by a time-dependent voltage source in
the circuit that linearly interpolates time-discrete data obtained
from the field model.
We will study two simulation scenarios for the magnetostatic
and the magnetoquasistatic case with varying
kL = [1, 0.9, 0.8, 0.7, 0.6, 0.5].
This reduces the inductance extracted for a? = a0 and can be
interpreted as an increase of magnetic saturation with respect
to the situation optimized for t0 .
Both scenarios account for the non-linear iron yoke characteristic, whereas only the latter incorporates the inter-filament
coupling currents, i.e., non-vanishing τeq . The coupling is
performed by means of the improved waveform relaxation
scheme described in Section III-D. The time interval is split
into 38 windows Ij = [t̄j , t̄j+1 ]. Each window has fixed length
Hj = t̄j+1 − t̄j = 20ms and communication between solver
occurs at time stamps t̄j . Iterations k of each window have
been carried out until two subsequent currents are below a
threshold
R t̄j+1 (k)
(k−1)
|im (t) − im
(t)|dt
t̄j
≤ 10−3
R t̄j+1 (k)
|i
(t)|dt
m
t̄j
C. Magnetostatic model
The first simulation scenario considers the magnetostatic
case, i.e., τeq ≡ 0 in the partial differential equation (1). This
simplified case disregards eddy currents and the iron yoke
is not saturated at the considered current level such that the
magnetic material behaves linearly. These simple models have
been traditionally considered as a first approximation in the
simulation of quench protection systems [24]. The simplification facilitates the extraction of an equivalent inductance Lm
according to (15), which perfectly describes the behavior of
the partial differential model in the circuit.
Fig. 4 summarizes five simulations with varying kL =
[1, 0.9, 0.8, 0.7, 0.6] and the given communication time window
size of 20 ms. The reference case kL = 1 is characterized by

∆vm [V]

where im (t) is the current flowing through the magnet,

0

−100

kL
kL
kL
kL
kL
kL

−200

0

0.2

0.4

=1
= 0.9
= 0.8
= 0.7
= 0.6
= 0.5

0.6

0.8

t [s]
Fig. 6. Voltage differences ∆vm correcting the mismatch between field and
circuit model in the magnetoquasistatic case

kL
kL
kL
kL
kL
kL

15

Inner iteration

calculated once for the initial value a0 = a(t0 ). The loop
voltage reads

im (t) REE + Rm (t) + ∆vm (t) + kL Lm ∂t im (t) = 0, (20)

100

10

=1
= 0.9
= 0.8
= 0.7
= 0.6
= 0.5

5

1

2

3

Outer time window
Fig. 7. Convergence comparison for first 3 time windows of the magnetoquasistatic case

a voltage source of ∆vm = 0 V. This verifies numerically that
α = (jωLm )−1 leads to optimal convergence for the magnetostatic model, i.e., it corresponds to an optimal contraction factor
of ρ(α) = 0. For the remaining cases, the decrease of kL is
compensated by an increased magnitude of ∆vm (t). Additional
effort needed to obtain convergence of the voltage generator
is reflected in the increase of number of iterations as shown
in Fig. 5. The case kL = 0.5 is not depicted, as the iterative
scheme diverged.
D. Magnetoquasistatic model
The second simulation scenario considers the thermal/magnetoquasistatic case using a non-linear time constant
τeq (B) as described in Section IV-A. The inter-filament coupling losses are effectively decreasing the flux linkage and as
a consequence the induced voltage in the field model.
Fig. 6 shows the simulated voltage differences. For the
nominal inductance, i.e. kL = 1, the voltage source ∆vm (t)

is initially positive and decreases to zero in the consecutive
time windows. Fig. 7 shows the number of iterations needed to
obtain a sufficiently accurate solution. It is worth noticing that
for kL = 0.9 the convergence is obtained in the least number
of iterations as compared to other cases. The reason for such
convergence lays in inter-filament coupling losses, which dissipate the magnetic energy in the coil. These processes results
can be attributed to the decrease of the effective differential
inductance and explains the need to inject missing energy by
means of ∆vm (t). On the other hand, the case of kL = 0.5
still converges and leads to correct results but requires up to
k = 16 iterations instead of the k = 2 ones for the optimal
case.
V. C ONCLUSION
This paper has discussed multi-physical field/circuit waveform relaxation for the specific eddy-current model used in
quench simulation. In contrast to previous works, no numerical
optimization of the transmission conditions was carried out to
obtain optimal convergence rates. Instead, it has been proven
for low frequencies that the contraction factor is significantly
reduced by considering an equivalent inductance or higherorder terms of a Neumann series. Numerical simulations of
an aperture dipole magnet underline the importance of an
optimization of the coupling conditions, as iterations could
be reduced from 16 to 2 per window by using an inductive
coupling.
Future research will investigate improved circuit models for
loss prediction.
ACKNOWLEDGMENT
This work has been supported by the Excellence Initiative of
the German Federal and State Governments and the Graduate
School of CE at TU Darmstadt.
R EFERENCES
[1] K. Dahlerup-Petersen, R. Denz, J. L. Gomez-Costa,
D. Hagedorn, P. Proudlock, F. Rodrinuez-Mateos,
R. Schmidt, and F. Sonnemann, “The protection system for
the superconducting elements of the large hadron collider
at cern,” in Proceedings of the 1999 Particle Accelerator
Conference (Cat. No.99CH36366), vol. 5, 1999, pp. 3200–
3202.
[2] E. Ravaioli, “Cliq – a new quench protection technology
for superconducting magnets,” Ph.D. dissertation, University of Twente, 2015.
[3] T. Dreher and G. Meunier, “3d line current model of coils
and external circuits,” IEEE Trans. Magn., vol. 31, no. 3,
pp. 1853–1856, May 1995.
[4] G. Bedrosian, “A new method for coupling finite element
field solutions with external circuits and kinematics,” IEEE
Trans. Magn., vol. 29, no. 2, pp. 1664–1668, 1993.
[5] E. Lange, F. Henrotte, and K. Hameyer, “An efficient fieldcircuit coupling based on a temporary linearization of FE
electrical machine models,” IEEE Trans. Magn., vol. 45,
no. 3, pp. 1258–1261, 2009.

[6] P. Zhou, D. Lin, W. N. Fu, B. Ionescu, and Z. J. Cendes,
“A General Co-Simulation Approach for Coupled FieldCircuit Problems,” IEEE Trans. Magn., vol. 42, no. 4, pp.
1051–1054, Apr. 2006.
[7] S. Schöps, H. De Gersem, and A. Bartel, “A cosimulation
framework for multirate time-integration of field/circuit
coupled problems,” IEEE Trans. Magn., vol. 46, no. 8,
pp. 3233–3236, Jul. 2010.
[8] A. Bartel, M. Brunk, M. Günther, and S. Schöps, “Dynamic iteration for coupled problems of electric circuits
and distributed devices,” SIAM J. Sci. Comput., vol. 35,
no. 2, pp. B315––B335, Mar. 2013.
[9] E. Lelarasmee, A. E. Ruehli, and A. L. SangiovanniVincentelli, “The waveform relaxation method for timedomain analysis of large scale integrated circuits,” IEEE
Trans. Comput. Aided. Des. Integrated Circ. Syst., vol. 1,
no. 3, pp. 131–145, 1982.
[10] K. Burrage, Parallel and sequential methods for ordinary
differential equations. Oxford: Oxford University Press,
1995.
[11] M. Al-Khaleel, M. J. Gander, and A. E. Ruehli, “Optimization of transmission conditions in waveform relaxation techniques for RC circuits,” SIAM J. Numer. Anal.,
vol. 52, no. 2, pp. 1076–1101, 2014.
[12] J. d. D. Nshimiyimana, F. Plumier, P. Dular, and
C. Geuzaine, “Co-simulation of of finite element and
circuit solvers using optimized waveform relaxation,” in
IEEE International Energy Conference (ENERGYCON)
2016, 2016, pp. 1–6.
[13] A. P. Verweij, “Electrodynamics of superconducting cables in accelerator magnets,” Ph.D. dissertation, Universiteit Twente, Twente, The Netherlands, 1995.
[14] S. Schöps, H. De Gersem, and T. Weiland, “Winding
functions in transient magnetoquasistatic field-circuit coupled simulations,” COMPEL, vol. 32, no. 6, pp. 2063–
2083, Sep. 2013.
[15] P. Monk, Finite Element Methods for Maxwell’s Equations. Oxford: Oxford University Press, 2003.
[16] C.-W. Ho, A. E. Ruehli, and P. A. Brennan, “The modified
nodal approach to network analysis,” IEEE Trans. Circ.
Syst., vol. 22, no. 6, pp. 504–509, Jun. 1975.
[17] A. Bartel, S. Baumanns, and S. Schöps, “Structural
analysis of electrical circuits including magnetoquasistatic
devices,” APNUM, vol. 61, pp. 1257–1270, Sep. 2011.
[18] S. Schöps, “Multiscale modeling and multirate timeintegration of field/circuit coupled problems,” Dissertation,
Düsseldorf, May 2011, VDI Verlag. Fortschritt-Berichte
VDI, Reihe 21. [Online]. Available: http://elpub.bib.
uni-wuppertal.de/servlets/DocumentServlet?id=2132
[19] T. Nakamoto, M. Sugano, Q. Xu, H. Kawamata,
S. Enomoto, N. Higashi, A. Idesaki, M. Iio, Y. Ikemoto,
R. Iwasaki, N. Kimura, T. Ogitsu, N. Okada, K. i. Sasaki,
M. Yoshida, and E. Todesco, “Model magnet development
of d1 beam separation dipole for the hl-lhc upgrade,” IEEE
Trans. Appl. Super., vol. 25, no. 3, pp. 1–5, Jun. 2015.
[20] COMSOL, COMSOL Multiphysics Reference Manual,
5th ed., 2016. [Online]. Available: http://www.comsol.com

[21] “Lhc design report,” CERN, Report CERN-2004-003-V1, 2004.
[22] OrCAD, PSpice User’s Guide, 17th ed., 2016.
[Online]. Available: http://www.orcad.com/products/
orcad-pspice-designer/overview
[23] L. Bortot, M. Maciejewski, A. M. Prioli, Marco Fernandez Navarro, J. B. Ghini, B. Auchmann, and A. P. Verweij,
“A consistent simulation of electro-thermal transients in
accelerator circuits,” IEEE Trans. Appl. Super., 2016.
[24] L. Rossi and M. Sorbi, “Qlasa: A computer code for
quench simulation in adiabatic multicoil superconducting
windings,” Istituto Nazionale di Fisica Nucleare, Milano,
Italy, Report INFN/TC-04/13, Jul. 2004.

