Optimized Field/Circuit Coupling for the Simulation of Quenches in Superconducting Magnets

TL;DR

This paper introduces an optimized coupling method for simulating quench transients in superconducting magnets, enhancing convergence and efficiency in multi-physics, multi-scale simulations of magnetothermal phenomena.

Contribution

It develops a novel field/circuit coupling framework with an equivalent magnet model, improving simulation convergence for superconducting magnet transients.

Findings

Enhanced convergence of the iterative coupling scheme.

Effective simulation of magnetothermal transients during faults.

Validated approach on an accelerator dipole magnet example.

Abstract

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.

Equations

$$\nabla \times (\nu \nabla \times \vec{A})+\nabla \times \left({\nu }_{\mathrm{s}}{\tau }_{\mathrm{eq}}\nabla \times \frac{\partial \vec{A}}{\partial t}\right)={\vec{\chi }}_{\mathrm{s}}\mathbf{i}$$

$${\nu }_{\mathrm{s}}^{-1}{\vec{M}}_{\mathrm{s}}={\tau }_{\mathrm{eq}}\frac{\partial \vec{B}}{\partial t},$$

$$\mathbf{v}=\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{\Psi },\text{ with }\mathbf{\Psi }={\int }_{\Omega }{\vec{\chi }}_{\mathrm{s}}^{\top }\vec{A}\mathrm{d}\Omega .$$

$$\rho C_{\mathrm{p}}\frac{\partial T}{\partial t}-\nabla \cdot (k\nabla T)=P_{\mathrm{s}}+P_{\mathrm{Joule}},$$

$$P_{\mathrm{Joule}}=q_{\mathrm{flag}}{\sigma }^{-1}\Vert {\vec{J}}_{\mathrm{s}}{\Vert }^2$$

$$P_{\mathrm{s}}={\vec{M}}_{\mathrm{s}}\cdot \frac{\partial \vec{B}}{\partial t}.$$

$${\mathbf{v}}_{\text{t}}=\mathbf{Ri}\text{ with }\mathbf{R}=q_{\mathrm{flag}}{\int }_{{\Omega }_{\text{s}}}{\vec{\chi }}_{\mathrm{s}}^{\top }{\sigma }^{-1}{\vec{\chi }}_{\mathrm{s}}\mathrm{d}\Omega .$$

$$\displaystyle\mathbf{K}_{\mathrm{s}}\bigl{(}\mathbf{a}\bigr{)}\frac{\mathrm{d}\mathbf{a}}{\mathrm{d}t}+\mathbf{K}_{\nu}\bigl{(}\mathbf{a}\bigr{)}\mathbf{a}-\mathbf{X}\mathbf{i}$$

$$\mathbf{v}-{\mathbf{X}}^{\top }\frac{\mathrm{d}\mathbf{a}}{\mathrm{d}t}$$

$${\mathbf{M}}_{\rho }\frac{\mathrm{d}\mathbf{t}}{\mathrm{d}t}+{\mathbf{K}}_{\mathrm{k}}\mathbf{t}-\mathbf{q}(\mathbf{a})$$

$${\mathbf{v}}_{\text{t}}-\mathbf{R}(\mathbf{t},\mathbf{a},\mathbf{i})\mathbf{i}$$

$$\displaystyle\mathbf{A}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}+\mathbf{B}\bigl{(}\mathbf{x}\bigr{)}\mathbf{x}+\mathbf{P}\mathbf{i}$$

$${\mathbf{P}}^{\top }\mathbf{x}-\mathbf{v}$$

$${\mathbf{Q}}^{\top }\mathbf{x}-{\mathbf{v}}_{\text{t}}$$

$${\mathbf{A}}_{\mathrm{c}}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{x}}^{(k+1)}+{\mathbf{B}}_{\mathrm{c}}{\mathbf{x}}^{(k+1)}+\mathbf{P}{\mathbf{i}}_{\mathrm{c}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{c}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{t},\mathrm{c}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{c}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{t},\mathrm{c}}^{(k+1)}$$

$${\mathbf{K}}_{\mathrm{s}}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{a}}^{(k+1)}+{\mathbf{K}}_{\nu }{\mathbf{a}}^{(k+1)}$$

$${\mathbf{X}}^{\top }\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{a}}^{(k+1)}$$

$${\mathbf{M}}_{\rho }\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{t}}^{(k+1)}+{\mathbf{K}}_{\mathrm{k}}{\mathbf{t}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{t},\mathrm{m}}^{(k+1)}$$

$${\mathbf{i}}_{\mathrm{m}}^{(k+1)}$$

$${\mathbf{A}}_{\mathrm{c}}j\omega {\mathbf{x}}^{(k+1)}+{\mathbf{B}}_{\mathrm{c}}{\mathbf{x}}^{(k+1)}+\mathbf{P}{\mathbf{i}}_{\mathrm{c}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{c}}^{(k+1)}$$

$${\mathbf{K}}_{\mathrm{s}}j\omega {\mathbf{a}}^{(k+1)}+{\mathbf{K}}_{\nu }{\mathbf{a}}^{(k+1)}$$

$${\mathbf{X}}^{\top }j\omega {\mathbf{a}}^{(k+1)}$$

$${\mathbf{i}}_{\mathrm{m}}^{(k+1)}$$

$${\mathbf{v}}_{\mathrm{c}}^{(k+1)}=\alpha {\mathbf{i}}_{\mathrm{c}}^{(k+1)}-\alpha {\mathbf{i}}_{\mathrm{m}}^{(k)}+{\mathbf{v}}_{\mathrm{m}}^{(k)},$$

$$||{\mathbf{v}}_{\mathrm{c}}^{(k+1)}-{\mathbf{v}}_{\mathrm{c}}^{(k)}||=|\rho (\alpha )|||{\mathbf{v}}_{\mathrm{c}}^{(k)}-{\mathbf{v}}_{\mathrm{c}}^{(k-1)}||.$$

$$\rho(\alpha)=\left(\textbf{I}+\alpha\mathbf{x}_{\mathrm{P}}^{-1}(\omega)\right)^{-1}\big{(}\alpha-\mathbf{Z}(\omega)\big{)}\mathbf{x}_{\mathrm{P}}^{-1}(\omega),$$

$$\mathbf{Z}(\omega )=j\omega {\mathbf{X}}^{\top }({\mathbf{K}}_{\mathrm{s}}j\omega +{\mathbf{K}}_{\nu }{)}^{-1}\mathbf{X}.$$

$$\alpha =\mathbf{Z}(\omega ).$$

$${\mathbf{v}}_{\mathrm{c}}^{(k+1)}=\mathbf{Z}(\omega ){\mathbf{i}}_{\mathrm{c}}^{(k+1)}-\mathbf{Z}(\omega ){\mathbf{i}}_{\mathrm{m}}^{(k)}+{\mathbf{v}}_{\mathrm{m}}^{(k)}.$$

$$\mathbf{Z}(\omega )$$

$$\kappa :=||-j\omega {\mathbf{K}}_{\nu }^{-\frac{1}{2}}{\mathbf{K}}_{\mathrm{s}}{\mathbf{K}}_{\nu }^{-\frac{1}{2}}||<1.$$

Full text

I Introduction

Superconducting 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 [2] and the coupling-loss induced quench system [3] 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 [4] or it is coupled ‘weakly’ using additional iterations or exploiting temporary linearizations [5, 6, 7]. In most cases, the coupling is established via current sources, voltage sources or inductances. In [8, 9] the various approaches have been interpreted in the context of waveform relaxation methods [10, 11] 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. [12, 13]. 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.

II Modeling

Starting from Maxwell’s equations and assuming a magnetoquasistatic (MQS) setting, i.e., neglecting displacement currents, the following partial differential equation

[TABLE]

is obtained on a domain $\Omega\in\mathbb{R}^{3}$ with an initial value $\vec{A}_{0}$ at time $t_{0}$ and with Dirichlet and Neumann boundary conditions on $\Gamma_{\text{dir}}\cup\Gamma_{\text{neu}}=\Gamma:=\partial\Omega$, respectively. The eddy currents are taken into account in the coil domain $\Omega_{\text{s}}$ by homogenization [3, 14], see Fig. 2. $\vec{A}$ is the magnetic vector potential, $\nu$ and $\nu_{\text{s}}$ are inverse magnetic permeabilities that may depend on $\mathit{B}=\|\vec{B}\|$ with $\vec{B}=\nabla\times\vec{A}$ the magnetic flux density and $\mathbf{i}$ the lumped currents through each coil. In $\Omega_{\text{s}}$, the time constant $\tau_{\mathrm{eq}}$ is defined such that the cable magnetization $\vec{M}_{\mathrm{s}}$ satisfies the relation

[TABLE]

where $\tau_{\mathrm{eq}}$ also depends on $B$. Finally, $\vec{\chi}_{\mathrm{s}}:\Omega\mapsto\mathbb{R}^{3\times n_{\text{s}}}$, with $n_{\text{s}}$ the number of coils, is the stranded-conductor winding function (see [15]), which homogeneously distributes the current in the domain $\Omega_{\text{s}}$ such that the current density $\vec{J}_{\mathrm{s}}=\vec{\chi}_{\mathrm{s}}\mathbf{i}$. The field equation (1) is coupled to a circuit equation model via the flux linkage

[TABLE]

The temperature in the domain $\Omega_{\text{s}}$ can be obtained from the heat balance equation

[TABLE]

with suitable boundary and initial conditions, where $\rho$ is the mass density, $C_{\mathrm{p}}$ the heat capacity, $k$ the thermal conductivity and $T$ the temperature. The Joule losses $P_{\mathrm{Joule}}$ are defined as

[TABLE]

with $\sigma$ being the homogenized non-linear conductivity of the cables in the case of a quench. It is activated by the sigmoid-type activation function $q_{\mathrm{flag}}$ that depends on time, the magnetic flux density $\vec{B}$ and current density $\vec{J}_{\mathrm{s}}$. $P_{\mathrm{s}}$ is the power density in the superconductor and relates the heat equation with the magnetoquasistatic field solution through

[TABLE]

In case of quenching, i.e., $q_{\mathrm{flag}}>0$, the superconducting coils have an Ohmic resistance with the voltage drop

[TABLE]

The non-linear resistance $\mathbf{R}$ inherits the dependencies on $t$, $\vec{B}=\nabla\times\vec{A}$ and $\vec{J}_{\text{s}}=\vec{\chi}_{\mathrm{s}}\mathbf{i}$.

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

[TABLE]

is obtained, where $\mathbf{a}$ denotes the degrees of freedom for the magnetic vector potentials, $\mathbf{t}$ the discrete temperatures, $\mathbf{q}(\mathbf{a})$ the discretized Joule losses, $\mathbf{M}_{\star}$ and $\mathbf{K}_{\star}$ the discretized material and differential operator matrices and $\mathbf{X}$ is the discretization of $\vec{\chi}_{\text{s}}$. Initial conditions $\mathbf{v}_{0}$, $\mathbf{v}_{\mathrm{t,0}}$, $\mathbf{a}_{0}$ and $\mathbf{t}_{0}$ at time $t_{0}$ are set. The system is differential algebraic, as $\mathbf{K}_{\mathrm{s}}$ is a singular matrix due to the fact that $\tau_{\mathrm{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

[TABLE]

with initial conditions $\mathbf{x}(t_{0})=\mathbf{x}_{0}$ and $\mathbf{i}(t_{0})=\mathbf{i}_{0}$. In the case of modified nodal analysis (see [17]), which is used in most modern circuit simulators, the degrees of freedom $\mathbf{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 $\mathbf{v}$ and $\mathbf{v}_{\text{t}}$).

III Waveform Relaxation

Following [8], 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 $\bar{t}_{0},\ldots,\bar{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 Gauss-Seidel-type scheme on a time window $\mathcal{I}_{j}=[\bar{t}_{j},~{}\bar{t}_{j+1}]$ using current and voltages as ‘transmission conditions’ is given by System 1 (circuit, cf. (6))

[TABLE]

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

[TABLE]

The algorithm consists of the following steps

set $\mathbf{v}_{\mathrm{m}}^{(0)}(t)=\mathbf{v}_{0}$ and $\mathbf{v}_{\mathrm{t,m}}^{(0)}(t)=\mathbf{v}_{\mathrm{t},0}$ on $t\in\mathcal{I}_{0}$

• 1)

solve (7) for given inputs $\mathbf{v}_{\mathrm{m}}^{(k)}(t)$ and $\mathbf{v}_{\mathrm{t,m}}^{(k)}(t)$ on $\mathcal{I}_{j}$ to obtain the output $\mathbf{i}_{\mathrm{c}}^{(k+1)}(t)$

• 2)

solve (8) for given input $\mathbf{i}_{\mathrm{c}}^{(k+1)}(t)$ on $\mathcal{I}_{j}$ to obtain the outputs $\mathbf{v}_{\mathrm{m}}^{(k+1)}(t)$ and $\mathbf{v}_{\mathrm{t,m}}^{(k+1)}(t)$

• 3)

if not converged, set $k=k+1$ and go to step 1), otherwise proceed with $\mathbf{v}_{\mathrm{m}}^{(0)}(t)=\mathbf{v}_{\mathrm{m}}^{(k+1)}(\bar{t}_{\mathrm{j}})$, $\mathbf{v}_{\mathrm{t,m}}^{(0)}(t)=\mathbf{v}_{\mathrm{t,m}}^{(k+1)}(\bar{t}_{\mathrm{j}})$, for $t\in\mathcal{I}_{\mathrm{j+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 [18]. 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 ‘zero-order’ Dirichlet transmission conditions. As in [12] and [13], 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$.

III-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 $\nu$ and $\nu_{\text{s}}$, as well as the time constant $\tau_{\mathrm{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 $\omega=2\pi f$. Then

[TABLE]

and

[TABLE]

For the following optimization of convergence, the first transmission condition (9c) is generalized to the linear combination

[TABLE]

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 $\alpha$ influences the convergence of the waveform relaxation scheme.

III-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. [11]. This is expressed by a contraction factor $\rho(\alpha)<1$, such that

[TABLE]

Assuming that $(\mathbf{K}_{\mathrm{s}}j\omega+\mathbf{K}_{\nu})$ is invertible, which is necessary for the solvability of the subsystem, we obtain

[TABLE]

where we assume that the transfer impedance $\mathbf{x}_{\mathrm{P}}(\omega)=\mathbf{P}^{\top}(\mathbf{A}_{\mathrm{c}}j\omega+\mathbf{B}_{\mathrm{c}})^{-1}\mathbf{P}$ exists and is invertible and the impedance of the magnetoquasistatic system

[TABLE]

Optimal convergence is attained for $\rho(\alpha)=0$, that is

[TABLE]

Therefore, the optimized transmission condition for the circuit is given by

[TABLE]

In frequency domain, the impedance $\mathbf{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, $\mathbf{Z}$ contains time derivatives which need to be approximated.

III-C Approximation of the optimal condition by inductances

Let us assume that $\mathbf{K}_{\nu}$ is invertible. Then, we can rewrite the impedance as

[TABLE]

if the operating frequency $\omega$ is low, as it must hold

[TABLE]

To analyze the smallness of $\kappa$, we make some further assumptions: let the magnetoquasistatic system be discretized, such that the curl-curl operators can be factorized as $\mathbf{K}_{\star}=\mathbf{C}^{\top}\mathbf{M}_{\star}\mathbf{C}$ with diagonal material matrices $\mathbf{M}_{\star}$. Then, it can be shown that

[TABLE]

holds, since

[TABLE]

with $\mathbf{M}_{\tau_{\mathrm{max}}}=\max\limits_{i}(\mathbf{M}_{\tau})_{ii}\mathbf{I}=\mathbf{M}_{\tau}+\bar{\mathbf{M}}_{\tau}$, and $\bar{\mathbf{M}}_{\tau}$ a positive semidefinite diagonal matrix. $\mathbf{M}_{\nu}=\mathbf{M}_{\nu_{\mathrm{s}}}+\bar{\mathbf{M}}_{\nu_{\mathrm{s}}}$, with $\bar{\mathbf{M}}_{\nu_{\mathrm{s}}}$ also positive semidefinite, as $(\mathbf{M}_{\nu_{\mathrm{s}}})_{ii}=(\mathbf{M}_{\nu})_{ii}$ in the coil domain and zero everywhere else. Therefore the series converges if

[TABLE]

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

[TABLE]

Using the coefficient $\alpha=j\omega\mathbf{L}$ and dividing the transmission condition (11) by $j\omega$ yields

[TABLE]

with the magnetic flux linkage $\boldsymbol{\Psi}_{\star}=\frac{1}{j\omega}\mathbf{v}_{\star}$. The inductance can be obtained by solving the magnetostatic problem, i.e. (1), with $\tau=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 [19]. 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.

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

[TABLE]

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

[TABLE]

is extracted at a working point $\mathbf{a}_{\star}$ from the differential curl-curl matrix \mathbf{K}_{\nu,\mathrm{d}}=\frac{\mathrm{d}}{\mathrm{d}\mathbf{a}}\Bigl{(}\mathbf{K}_{\nu}\bigl{(}\mathbf{a}\bigr{)}\mathbf{a}\Bigr{)} and it is kept constant for the subsequent waveform relaxation, i.e., the condition (9c) is replaced by

[TABLE]

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

[TABLE]

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

[TABLE]

where $\mathbf{R}^{(k)}(t)=\mathbf{R}(t,\mathbf{a}^{(k)},\mathbf{i}_{\mathrm{m}}^{(k)})$. The Gauss-Seidel-type scheme using these new transmission conditions reads

set $k=0$, $\Delta\mathbf{v}_{\mathrm{m}}^{(0)}(t)=\Delta\mathbf{v}_{\mathrm{m,0}}(t)$ and and R =R() $\mathbf{R}^{(0)}(t)=\mathbf{R}(t_{0},\mathbf{a}_{0},\mathbf{i}_{0})$ and extract $\mathbf{L}_{\text{m}}$ for $\mathbf{a}_{\star}=\mathbf{a}_{0}$ using (15)

• 1)

solve (7a-7c) with conditions (16) and (18) for given inputs $\Delta\mathbf{v}_{\text{m}}^{(k)}(t)$ and $\mathbf{R}^{(k)}(t)$ to obtain the output $\mathbf{i}_{\mathrm{c}}^{(k+1)}(t)$

• 2)

solve (8) for given input $\mathbf{i}_{\mathrm{c}}^{(k+1)}(t)$ to obtain the output voltage $\Delta\mathbf{v}_{\text{m}}^{(k)}(t)$ according to (16) and $\mathbf{R}^{(k+1)}(t)$ as before

• 3)

if not converged, set $k=k+1$ and go to step 1), otherwise proceed with $\Delta\mathbf{v}_{\mathrm{m}}^{(0)}(t)=\Delta\mathbf{v}_{\mathrm{m}}^{(k+1)}(\bar{t}_{\mathrm{j}})$, $\mathbf{R}^{(0)}(t)=\mathbf{R}^{(k+1)}(\bar{t}_{\mathrm{j}})$, $j=j+1$ and go to step 1).

Obviously, the algorithm can be improved by appropriate heuristics to update $\mathbf{L}_{\text{m}}$ if the number of iterations $k$ grows.

IV Numerical Examples

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 [20].

IV-A Discussion of the time constant

The time constant $\tau_{\mathrm{eq}}$ depends on the cable parameters and is mathematically expressed by [14]

[TABLE]

with $l_{\mathrm{f}}$ being the filament twist-pitch, $\mu_{0}$ the vacuum permeability, $f_{\mathrm{eff,x}}$ the fraction of superconductor in the matrix, $c_{0}$ and $c_{1}$ being characterized by the resistivity of the material used in the matrix and $0<B\leq 10$ being the magnitude of the magnetic flux density.

Typical values for those parameters are given in [3]. For the case of a Nb-Ti dipole accelerator magnet, one has

[TABLE]

and thus it follows

[TABLE]

which yields convergence of the Neumann series for frequencies below $f_{\mathrm{max}}<7.5$ Hz.

IV-B Simulation of dipole magnet

The described method is illustrated for a numerical example of the single aperture dipole magnet D1 [20]. 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

[TABLE]

and $k\nabla T|_{\Gamma_{\text{s}}}\cdot\vec{n}=0$ W/m2, where $\vec{n}$ is the outward normal vector. The partial differential equations (5a)-(5d) are discretized and solved using COMSOL Multiphysics® [21]. 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 $\delta t=1$ ms). All parameters, e.g. non-linear materials laws for the reluctivity $\nu$ and time constant $\tau$, are defined as specified in [20, 22]. We consider a magnet operating at $5$ kA and protected by a resistor $R_{\text{EE}}=0.1\,\Omega$ as depicted in Fig. 3. The circuit is modeled and simulated in OrCAD PSpice® using the trapezoidal method (maximum time step $\delta t=10\,\mu$s) [23]. The waveform relaxation is carried out on $\mathcal{I}=(0\,\mathrm{s},0.75\,\mathrm{s}]$ with initial values $\vec{A}_{0}$, $T_{0}$ and $\mathbf{x}_{0}$ that correspond to a system ramped 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 [math] to $5$ kA. The co-simulation was established within CERN’s in-house coupling tool STEAM (Simulation of Transient Effects in Accelerator Magnets) [24].

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 $k_{L}$ to be multiplied by the (scalar) differential inductance

[TABLE]

calculated once for the initial value $\mathbf{a}_{0}=\mathbf{a}(t_{0})$. The loop voltage reads

[TABLE]

where $i_{\text{m}}(t)$ is the current flowing through the magnet,

[TABLE]

is the time-dependent coil resistance appearing after a quench and $\Delta v_{\text{m}}(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

[TABLE]

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

[TABLE]

This reduces the inductance extracted for $\mathbf{a}_{\star}=\mathbf{a}_{0}$ and can be interpreted as an increase of magnetic saturation with respect to the situation optimized for $t_{0}$.

Both scenarios account for the non-linear iron yoke characteristic, whereas only the latter incorporates the inter-filament coupling currents, i.e., non-vanishing $\tau_{\mathrm{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 $\mathcal{I}_{j}=[\bar{t}_{j},\bar{t}_{j+1}]$. Each window has fixed length $H_{j}=\bar{t}_{j+1}-\bar{t}_{j}=20\mathrm{ms}$ and communication between solver occurs at time stamps $\bar{t}_{j}$. Iterations $k$ of each window have been carried out until two subsequent currents are below a threshold

[TABLE]

IV-C Magnetostatic model

The first simulation scenario considers the magnetostatic case, i.e., $\tau_{\mathrm{eq}}\equiv 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 [25]. The simplification facilitates the extraction of an equivalent inductance $L_{\text{m}}$ according to (15), which perfectly describes the behavior of the partial differential model in the circuit.

Fig. 4 summarizes five simulations with varying $k_{\text{L}}=[1,0.9,0.8,0.7,0.6]$ and the given communication time window size of $20$ ms. The reference case $k_{\text{L}}=1$ is characterized by a voltage source of $\Delta v_{\text{m}}=0$ V. This verifies numerically that $\alpha=(j\omega L_{\text{m}})^{-1}$ leads to optimal convergence for the magnetostatic model, i.e., it corresponds to an optimal contraction factor of $\rho(\alpha)=0$. For the remaining cases, the decrease of $k_{\text{L}}$ is compensated by an increased magnitude of $\Delta v_{\text{m}}(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 $k_{\text{L}}=0.5$ is not depicted, as the iterative scheme diverged.

IV-D Magnetoquasistatic model

The second simulation scenario considers the thermal/magnetoquasistatic case using a non-linear time constant $\tau_{\mathrm{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. $k_{\text{L}}=1$, the voltage source $\Delta v_{\text{m}}(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 $k_{\text{L}}=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 $\Delta v_{\text{m}}(t)$. On the other hand, the case of $k_{\text{L}}=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 Conclusion

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 higher-order 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.