Proceedings of the 7th GACM Colloquium on Computational Mechanics
for Young Scientists from Academia and Industry
October 11-13, 2017 in Stuttgart, Germany

Reduced Order Modelling for the Simulation of
Quenches in Superconducting Magnets

arXiv:1710.04859v1 [] 13 Oct 2017

Sebastian Schöps1 *, Idoia Cortes Garcia1 , Michal Maciejewski2 and Bernhard Auchmann2,3
Micro Abstract
This contributions discusses the simulation of magnetothermal effects in superconducting magnets
as used in particle accelerators. An iterative coupling scheme using reduced order models between a
magnetothermal partial differential model and an electrical lumped-element circuit is demonstrated.
The multiphysics, multirate and multiscale problem requires a consistent formulation and framework
to tackle the challenging transient effects occurring at both system and device level.
1

Graduate School CE, Technical University of Darmstadt, Darmstadt, Germany

2

Machine Protection and Electrical Integrity Group, CERN, Geneva, Switzerland

3

Magnet Group, Paul Scherrer Institut, Villigen, Switzerland

* Corresponding author: schoeps@temf.tu-darmstadt.de

1 Introduction
The LHC (Large Hadron Collider) at CERN, Geneva uses superconducting magnets for focussing
and bending particle beams. To ensure their superconductivity, they are cooled to very low
temperatures (1.9 K). However, in the case of a malfunction the superconducting material
becomes resistive and quenches. The quenched area quickly overheats, the heat propagates
within the magnet and from one magnet to another causing more magnets to quench. Eventually,
the whole network of magnets can be damaged. This is avoided by quench detection and
protection systems, e.g. Quench heaters [3] or CLIQ (Coupling-Loss Induced Quench), [11].
Quench protection analysis can be further improved by coupling controller and circuit models [9].

2 Magnet Models
The electromagnetic field inside each magnet is described by Maxwell’s equations [8], i.e.,
∇×E=−

∂B
,
∂t

∇×H=

∂D
+J,
∂t

∇·D=ρ,

∇·B=0

(1)

on a bounded domain Ω ⊂ R3 and time interval I = (t0 , tend ] with specified boundary and
initial conditions on ∂Ω and t0 . They are interlinked by the material relations
D = εE ,

H = νB ,

J = σE + Js .

(2)

Here, E is the electric field, D the electric flux density, H the magnetic field, B the magnetic
flux density, J, Js are the total, source current densities and ρ is the electric charge density. The
positive (semi-)definite material tensors σ  0, ε ≻ 0 and ν ≻ 0 define the electric conductivity,
the electric permittivity and the magnetic permeability, respectively. For simplicity of notation,
we consider only the case of linear material laws, see e.g. [2] for a generalization.
∂
D = 0 and introducing the magnetic
Assuming a magnetoquasistatic approximation with ∂t
vector potential B = ∇ × a, the parabolic, semi-elliptic partial differential equation

σ

∂a
+ ∇ × (ν∇ × a) = χim
∂t

(3)

is obtained, where the current density Js (r, t) = χ(r)im (t) is expressed in terms of scalar
currents im . In practice, one often replaces the term σ ∂a/∂t by a homogenized model to
avoid the resolution of small geometric details in the discretization, e.g. single strands within
cables [1]. Following [2], equation (3) is coupled via the voltage vm to a circuit using the flux
linkage
Z
χ(r) · a(r, t) dΩ.

Φ(t) =

(4)

Ω

The temperature in the domain Ωs can be obtained from the heat balance equation
ρCp

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

(5)

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 kJs k2 and are activated in case of a quench by the sigmoid-type activation
function qflag that depends on time, the magnetic flux density B = ∇ × a and current density
Js . The power density Ps in Ωs relates the heat equation with the field solution [2]. In case of
quenching, i.e., qflag > 0, the superconducting coils have a resistance with the voltage drop
Z
d
χ ⊤ σ −1 χ dΩ.
(6)
vm = Φ + Rt im with Rt = qflag
dt
Ωs
The resistance Rt inherits the dependencies on t, B and Js . As temperature increases, thermal
stress is caused which leads to displacement. The problem of linear elasticity reads
∇ · (C : ǫ) = F(T, J, B)

(7)


where ǫ = 21 ∇u + (∇u)⊤ denotes the strain tensor and u the displacement vector. Boundary
conditions still to be defined. The right hand side F depends on the thermal stress and Lorentz
forces. In turn, the computed displacement deforms the computational domain Ω.

3 Networks of Magnets
Let us assume that we have Nm magnets connected within a network. Then each magnet
n = 1, ..., Nm is given by the same multiphysical model which is described by the coupled
system of partial differential equations (3)-(7) defined above. We address the solution of each
model abstractly by the operator Ψn , such that vm,n = Ψn (im,n ) (n = 1, ..., Nm ), where the
internal variables an , Tn , un are suppressed as the electrical coupling to neighboring magnets
is established only via currents im,n and voltages vm,n . Let us collect all solution operators as
vm = Ψ(im )

(8)

with currents i⊤m = [im,1 , im,2 , ..., im,Nm ] and voltages v⊤m = [vm,1 , vm,2 , ..., vm,Nm ]. Let an electrical network consist of simple elements and described by the modified nodal analysis [7], i.e.,
AC CA⊤C

d
ϕ + AR GA⊤R ϕ + AL iL + AV iV +Am im = −AI is (t) ,
dt
d
L iL − A⊤L ϕ = 0 ,
dt
A⊤V ϕ = vs (t) ,
A⊤m ϕ = vm (t) .

(9)
(10)
(11)
(12)

The position of each device in the network is encoded by a sparse incidence matrix A⋆ which
consists only of entries from the set {−1, 0, 1}. The corresponding branch voltage is denoted by

v⋆ = A⊤⋆ ϕ and the current through the device is i⋆ . The simplest network elements are inductors,
capacitors and resistors; they are determined by lumped constitutive parameters, collected in
the matrices of capacitances C, inductances L and conductances G. Time dependent voltage
and current sources are given as vector-valued functions vs (t) and is (t). The unknowns of
the system are the scalar potentials ϕ and a few currents, i.e. iL and iv . More complicated
devices can be modeled by controlled sources, see e.g. [4]. Finally, the terms enclosed in the two
boxes indicate the additional (non-standard) contributions due to the field/circuit coupling via
equation (8).

4 Iterative Coupling Accelerated by Reduced Models
Instead of solving the large coupled system (8)-(12) monolithically, we propose an iterative
scheme that fits into the framework of waveform relaxation or dynamic iteration, [5, 13]. This
allows us to use a different discretization in time for each system and more importantly, even
different software packages can be used.
(k)

Let us assume that we have an initial guess of the time-transient behavior of the voltages vm (t)
for t ∈ I and k = 0. Then, in the simplest case, the following scheme can be performed
(k)

1. solve the electric network equations (9)–(12) by considering the voltages vm to be known
(k+1)
and extract the currents im
(k)

2. solve n = 1, ..., Nm electrothermal problems (3)-(6) in parallel, each on domain Ωn excited
(k+1)
(k+1)
(k+1)
.
by the current im,n . Extract voltages vm,n and temperatures Tn
(k)

3. solve n = 1, ..., Nm elasticity problems (7) in parallel, consider the temperature Tn to be
(k+1)
(k+1)
.
known, compute displacements um,n and deform the computational domains Ωn
4. if not converged: increase counter k = k + 1 and return to step 1.
Due to the smallness of the displacements it is common practice to exclude the computation
of the deformations from the iteration and solve the elasticity problem (7) only once in a postprocessing step, i.e., after the iteration of steps 1-2 and 4 has converged. Convergence for the
remaining electrothermal scheme is well understood and has been discussed based on Banach’s
fixed point theorem and an analysis of the contraction factor in [12], [2].
The iterative method simulates each magnet in parallel but the computational burden might
still be heavy. This can be mitigated by using reduced order models to decrease the number of
iterations and eventually to reduce simulation costs, [10]. Let us assume that there is a reduced
order model available that approximates the multiphysical problem (3)-(7). Then, the solution
of the n-th magnet can be rewritten as
vm,n = Ψn (im,n ) + ∆vm,n

(13)

where ∆vm,n is the deviation between the reduced model and the full one. This additional
reduced model allows us to define a waveform relaxation in terms of voltage differences, i.e.,
AC CA⊤C

d (k+1)
(k+1)
(k+1)
+Am i(k+1)
+ AV iV
ϕ
+ AR GA⊤R ϕ(k+1) + AL iL
= −AI is (t) ,
m
dt
d (k+1)
− A⊤L ϕ(k+1) = 0 ,
L iL
dt
A⊤V ϕ(k+1) = vs (t) ,
(k)
) = ∆vm
.
A⊤m ϕ − Ψ(i(k+1)
m

(14)
(15)
(16)
(17)

This scheme still requires the simulation of the full model Ψ but may significantly reduced the
number of evaluations. On the other hand, sophisticated reduced models Ψ may also increase
the computational costs (particularly in the offline stage). However, in the simplest case the
reduced model is given by an improved transmission condition in terms of the optimized Schwarz
waveform relaxation [6] [2].

5 Outlook
The full contribution will show simulations of the mechanical response of a magnet during quench
protection and magnet discharge. For CLIQ-protected magnets, there are two phenomena
in the coil: an imbalanced current profile leading to inhomogeneous Lorentz forces and an
almost homogeneous temperature profile [11]. We quantify the impact of Lorentz forces and
temperature gradients during a discharge as large stresses may lead to cable degradation over
time.

Acknowledgements
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.

References
[1] L. Bortot, M. Maciejewski, M. Prioli, A. Fernandez Navarro, S. Schöps, I. Cortes Garcia,
B. Auchmann, and A. Verweij. Simulation of electro-thermal transients in superconducting
accelerator magnets with comsol multiphysics. In Proceedings of the European COMSOL
Conference 2016, Munich, 2016.
[2] I. Cortes Garcia, S. Schöps, L. Bortot, M. Maciejewski, M. Prioli, A. Fernandez Navarro,
B. Auchmann, and A. Verweij. Optimized field/circuit coupling for the simulation of
quenches in superconducting magnets. JMMCT, 2(1):97–104, 2017.
[3] 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, volume 5, pages 3200–3202, 1999.
[4] D. Estévez Schwarz and C. Tischendorf. Structural analysis of electric circuits and consequences for MNA. Int. J. Circ. Theor. Appl., 28(2):131–162, 2000.
[5] M. J. Gander. Waveform relaxation. In B. Engquist, editor, Encyclopedia of Applied and
Computational Mathematics. Springer, 2012.
[6] M. J. Gander and A. E. Ruehli. Optimized waveform relaxation methods for RC type
circuits. IEEE Trans. Circ. Syst., 51(4):755–768, 2004.
[7] C.-W. Ho, A. E. Ruehli, and P. A. Brennan. The modified nodal approach to network
analysis. IEEE Trans. Circ. Syst., 22(6):504–509, 1975.
[8] J. D. Jackson. Classical Electrodynamics. Wiley and Sons, New York, 3rd edition, 1998.
[9] M. Maciejewski, I. Cortes Garcia, S. Schöps, B. Auchmann, L. Bortot, M. Prioli, and A. P.
Verweij. Application of the waveform relaxation technique to the co-simulation of power
converter controller and electrical circuit models. In MMAR 2017. IEEE, 2017.
[10] M. Rathinam and L. R. Petzold. Dynamic iteration using reduced order models: A method
for simulation of large scale modular systems. SIAM J. Numer. Anal., 40:1446–1474, 2002.
[11] E. Ravaioli. CLIQ – A new quench protection technology for superconducting magnets. PhD
thesis, University of Twente, 2015.

[12] S. Schöps, H. De Gersem, and A. Bartel. A cosimulation framework for multirate timeintegration of field/circuit coupled problems. IEEE Trans. Magn., 46(8):3233–3236, 2010.
[13] J. K. White, F. Odeh, A. L. Sangiovanni-Vincentelli, and A. E. Ruehli. Waveform relaxation: Theory and practice. Trans. Soc. Comput. Sim., 2(1):95–133, 1985.

