Quasi-3D magnetic ﬁeld simulation of superconducting devices with translational symmetry

This work presents a quasi-three-dimensional (Q3D) approach for the magnetic ﬁeld simulation in superconducting devices. First-order two-dimensional ﬁnite-element edge functions in the model’s cross-section are combined with one-dimensional orthogonal polynomials along the longitudinal direction. The interﬁlament coupling currents arising in superconducting multi-ﬁlament materials are modelled by taking the associated magnetization into account. For this formulation, the Q3D ansatz is elaborated, veriﬁcated and applied to a superconducting cable model. In the end, the approach is compared to a conventional three-dimensional ﬁnite-element method against which the proposed Q3D method demonstrates a superior computational efﬁciency.


INTRODUCTION
Accelerator technology is driven by the desire for higher and higher particle energies. Managing such high-energy particle beams requires either larger accelerator circumferences to adapt to the larger particle gyroradius, or higher magnetic fields to be able to deflect the beam more strongly [1]. While conventional normal-conducting magnets based on an iron yoke produce magnetic fields limited to 1.5T due to iron saturation, superconducting magnets can generate much higher magnetic fields, for example, up to 9T are achieved at the Large Hadron Collider (LHC) at CERN using superconducting magnets based on NbTi and Nb 3 Sn [2]. Superconductors are materials which lose their electric resistivity upon reaching the superconducting state, that is, below a critical parameter set of temperature, magnetic field and current density. Then, they are able to carry high current densities leading to high magnetic fields without producing Joule losses [3]. A main issue of superconductors is the so-called quench phenomenon [4]: When the critical temperature is surpassed in a region of the superconducting magnet, this region shifts to normal-conducting state meaning that the electric resistivity jumps to a value above zero. In normal operation, quench phenomena arise locally in the magnet and vanish due to cooling. In the worst case, however, a quench can lead to a thermal runaway resulting into severe damage or even the destruction of the magnet, for example, reported in [5]. This makes a reliable quench protection system mandatory, which is triggered upon quench detection, that is, when predefined voltage or resistance thresholds are exceeded [6]. Today, these thresholds are equipped with a large safety margin leading to frequent and unnecessary shutdowns, which considerably reduce the availability of the accelerator machine [7]. To reduce the amount of false triggers and therefore increase the availability of the accelerator, quench phenomena need to be understood better and quench detection has to get preciser. For this, numerical simulations play a crucial role, but they also prove to be challenging for computational engineers due to the multiphysical magneto-thermal and multi-scale nature of superconducting accelerator magnets. A promising approach to tackle this deeply intertwined problem is STEAM, a hierarchical co-simulation framework for superconducting accelerator magnets [8]. The idea is to decompose the whole problem into manageable subproblems which are solved independently from each other and exchange information. In this context, this work deals with the subproblem of an efficient spatial discretization. Especially the multi-scale nature is problematic and concerns both geometry and physics: While the magnet's cross-section contains many geometrical details in the FIGURE 1 (a) The cross-section of a superconducting dipole at the LHC is depicted, where the iron yoke is colored in gray and the coils either in blue or red depending on the current directions. (b) The superconducting coil configuration is illustrated showing the rectangular superconducting Rutherford-type cables and the trapezoidal copper keystones. (c) Three cables are shown in detail, which are wrapped by a glass fibre insulation (black) and contain multiple superconducting strands (circles). Figure adapted from [19] mm-scale (c.f. Figure 1), and also quench phenomena take place in a mm-scale [9], the magnet is several meters long [10]. Consequently, the use of conventional three-dimensional (3D) finiteelement (FE) methods would demand for a fine meshing in the transversal plane to cover all geometrical details, but also in longitudinal direction to resolve the transient quench phenomena. Inevitably, this leads to unacceptably time-consuming simulations [11]. Two-dimensional (2D) [12] or even one-dimensional (1D) [13] approximations are able to catch some important features, but cannot accurately reflect the 3D quench propagation behaviour.
While the superconducting magnet exhibits a 3D physical behaviour, the geometry shows a translational invariance in longitudinal direction. This enables a quasi-3D (Q3D) ansatz, in which the transversal cross-section is discretized with a 2D FE method using a triangular mesh and the longitudinal direction is resolved with a 1D spectral-element (SE) method using orthogonal polynomials [14]. While the FE ansatz functions are of lowest order, the order of the orthogonal polynomials is arbitrary, making this approach similar to an anisotropic hpmethod, for example, employed in [15] and [16]. The choice of the orthogonal polynomials is crucial as it determines the sparsity of the resulting matrices and therefore the computational costs, for example, an approach with Legendre polynomials produces a dense stiffness matrix and therefore proves to be suboptimal [17]. In contrast, the beneficial properties of the modified Lobatto polynomials [18] lead to sparse matrices and an overall easy handling on the implementation part. Such a Q3D method using modified Lobatto polynomials has been employed successfully for the thermal subproblem in [19] to compute the heat conduction in superconducting models and to show the superior efficiency compared to a conventional 3D FE method.
Now, this work focuses on the magnetic subproblem of quench simulation. While the Q3D thermal subproblem of quench propagation is discretized with hybrid nodal shape functions [19], the Q3D magnetic subproblem is more involved, as the magnetoquasistatic (MQS) formulation in terms of the magnetic vector potential (MVP) necessitates the use of edge func- On the left, the cross-section of a superconducting cable and one of its multi-filamentary strands are shown [23]. On the right, the effects of an external time-varying magnetic field on the superconducting filaments are illustrated in case of a non-conducting matrix (a) and a conducting matrix (b), respectively tions [20] and gauging [21]. Furthermore, so-called interfilamentary coupling currents (IFCCs) arise in a superconducting cable which cannot be expressed by the standard MQS formulation but need special modelling [22].
This work is structured as follows: First, the IFCC phenomenon is examined and a modified MQS formulation is derived. Next, a summary of the 2D FE edge functions is given and the modified Lobatto polynomials are introduced for the SE part. Then, the FE and SE parts are combined into a Q3D method and the discrete Q3D magnetic formulation is presented. The incorporation of boundary and gauging conditions is discussed. Next, the method is verificated against a test problem and a convergence study is carried out. The method is then applied on a superconducting cable model. Lastly, a conclusion and an outlook to further work is given. Figure 2 (left) shows the structure of a Rutherford-type cable, which is a superconducting cable with multi-filamentary strands. This cable is used in accelerator magnets at the LHC and consists of several strands composed of many superconducting filaments which are embedded in a copper matrix [23]. When an external time-varying magnetic field is applied to a superconductor, screening currents arise to generate a magnetization counteracting the external field such that the magnetic field within the superconductor is expelled [3]. If the superconducting filaments were embedded in a non-conducting matrix, each filament would carry its own screening currents as illustrated in Figure 2a. However, the conductive copper matrix causes the screening currents to couple the neighboring filaments via IFCCs as depicted in Figure 2b, resulting into a completely different magnetic field distribution compared to the non-conductive case [4].

MQS FORMULATION FOR SUPERCONDUCTING CABLES
It is not manageable to resolve the individual filaments in an overall cable model. Instead, the IFCC magnetization is modelled as where is the reluctivity in Am/Vs, is the IFCC time constant [9] and ⃗ B is the magnetic flux density in T [22], [24]. Using Faraday's law and by applying the rotation operator, the IFCC magnetization current density in A/m 2 is obtained, where ⃗ E denotes the electric field strength in V/m. The standard MQS formulation reads Hereby, the MQS formulation (3) becomes In the magnetostatic (MS) case with an impressed current density ⃗ J s , (3) reduces to To obtain the weak formulation, the formulation (6) and (7), respectively, are multiplied with test functions ⃗ A ′ and integrated over the domain V . Exemplarily, the weak formulation of the MS case (7) ultimately reads: where U is appropriately chosen as H (curl; V ) to correctly reflect the tangential continuity property of the MVP [25].

2D FE EDGE FUNCTIONS
For an appropriate discretization of (8), ansatz functions belonging to H (curl; V ) are needed. For the FE method, edge functions fulfill this requirement. In the 2D FE method on a Cartesian xy-plane, two types of edge functions are dis- The two types of edge functions in a 2D Cartesian mesh: Transversal edges (purple) are associated with the mesh edges e, while longitudinal edges (blue) perpendicularly emerge from the mesh nodes j tinguished: A transversal edge function is associated with an edge e = (i, j ) from a node i to a node j of the 2D mesh, while a longitudinal edge function is associated with a node j in an edge emerging perpendicularly from the xy-plane [26] as illustrated in Figure 3. The first-order edge functions can be expressed in terms of FE nodal functions N j (x, y) [27] as where the superscripts t and stand for the transversal and longitudinal type, respectively, z denotes the length of the geometry, and ⃗ e z is the unit vector in longitudinal direction.
Using the transversal edge functions (9) and adapting the Ritz-Galerkin procedure, that is, choosing the spaces of the test and ansatz functions identical, the transversal FE stiffness matrix and the transversal FE mass matrix are obtained, where Ω denotes the 2D cross-section of the computational domain. The transversal FE matrices have the dimensions N e × N e , where N e is the number of mesh edges. Analogously, the Ritz-Galerkin procedure applied on the longitudinal edge functions (10) yields for the longitudinal FE stiffness matrix which has the dimensions N n × N n with N n as number of mesh nodes. These matrices are computed on the reference trianglê T = {(0, 0), (1, 0), (0, 1)} and then mapped to the image triangles T . Here, it is crucial to use a H (curl; Ω)-conform mapping which is provided by a covariant Piola transformation [28]. Otherwise, the tangential continuity property of the MVP can not be maintained. For the discretization of (6), the MVP is decomposed into a transversal and a longitudinal component, where ⃗ A (t ) possesses only components parallel to the xy-plane and ⃗ A ( ) has only a z-component. Then, the MVP is discretized with the edge functions (9) and (10) as y). (15) The coefficients ⌢ a t e and ⌢ a j represent the line-integrated MVP in Wb on the transversal and longitudinal edges, respectively. This yields for (6) the semi-discrete system of equations where the vectors ⌢ a t and ⌢ a collect the coefficients in (15). Note that the transversal and longitudinal parts are decoupled from each other, resulting into two independent problems for the transversal and longitudinal MVP component, respectively.

MODIFIED LOBATTO POLYNOMIALS
Considering a general 1D differential equation a discrete system can be obtained by the SE method [14], where denote the spectral stiffness matrix, the spectral damping matrix, the spectral mass matrix and the spectral load vector, respectively, in terms of polynomials (k) q of q-th order on the k-th SE on an interval I . For the load vector, the right-hand side function f (z ) is forward transformed in terms of the orthogonal polynomials, which is expressed by the interpolant  k N for the k-th SE and polynomial order N . The solution is expanded as where N is the maximal polynomial order. In contrast to the FE method, the spectral coefficientsũ k q do not generally represent a physical value, but the weight of the q-th polynomial mode in element k. The polynomial's Vandermonde matrix [29] is used to backward transform a spectral solution into a physical one, or to forward transform a physical solution into a spectral one.
There are several possibilities for choosing the polynomials q (z ) which lead to different properties of the spectral matrices [14]. In this work, the modified Lobatto polynomials [18] are employed which are defined on the orthogonality intervalÎ = [−1, 1] aŝq Here,ẑ(z ) ∈Î is a mapped parameter and LO q (ẑ ) the Lobatto polynomial of q-th order [29]. One can map between the image interval I and reference intervalÎ with the usual transformation rules [14]. The modified Lobatto polynomials offer many advantages: The interior modes q = 2, … , N are modal, that is, the inner stiffness matrix is diagonal. Furthermore, the inner mass and damping matrix happen to be tridiagonal and two-banded, respectively. Thus, the spectral matrices in (18) are all sparse. Meanwhile, the boundary modes q = 1 and q = N + 1 are nodal, that is, the corresponding coefficientsũ 1 andũ N +1 directly represent the physical values at the SE boundariesẑ = −1 andẑ = 1, respectively. This makes the incorporation of boundary conditions straightforward and the continuity over SE interfaces easy to impose. Furthermore, one can find closed forms for the computation of (19)-(21) [18] which eliminates the need for numerical integration and speeds up the assembly.

Q3D DISCRETIZATION
In the Q3D setting, the transversal xy-cross-section Ω located at z = 0 is discretized with a 2D FE mesh with triangular elements  (9) and (10) with the polynomials (24), yielding with the transversal and longitudinal Q3D edge functions ⃗ v Then, the 3D MVP is discretized by Here, the use of a tilde instead of a bow indicates the spectral nature of the coefficients. To calculate the Q3D stiffness matrix for (6), the curl operator needs to be applied on (25), yielding Note the missing −1 z factor for the longitudinal edge functions, which is omitted due to the polynomial expansion in zdirection. Assuming that both material properties and can be expressed as a product of a transversal and longitudinal factor, that is, (x, y, z ) = (t ) (x, y) ( ) (z ), ∈ { , }, the Q3D stiffness matrix is obtained by with where ⊗ denotes the Kronecker tensor product [30] and is a matrix which couples transversal and longitudinal edges. In contrast to (16), here the transversal and longitudinal MVP components are coupled, such that a 3D behaviour is indeed described. Finally, the Q3D semi-discrete system for (6) reads where the vectorã collects the coefficients of (28). With the exception of (34), all occurring matrices in (30) are well known from the 2D FE and the 1D SE method. Furthermore, the Kronecker tensor product enables an efficient assembly.

BOUNDARY CONDITIONS
As combination of a 2D FE and a 1D SE method, the Q3D method distinguishes three boundary parts: For the 2D FE part, there is the boundary of the 2D mesh, which extrudes to a mantle surface Ω M in the Q3D setting, while for the 1D SE part, there are the end points of the 1D interval, which correspond to two end surfaces Ω F and Ω B at front and back of the geometry as illustrated in Figure 5. Boundary conditions (BCs) can only be assigned at MVP components tangential to the domain boundary. At Ω M , Dirichlet BCs can affect both the unknownsã t,k eq andã ,k jw where e is an edge and j is a node at Ω M . The mantle Dirichlet data ⃗ A(x, y, z )| Ω M = g M (x, y, z ) have to be forward transformed w.r.t. the modified Lobatto polynomials and need to be applied on all modes q, w = 1, … , N + 1 for each affected node j and edge e. In the special case that the data is constant along z, it follows that the inner modes q, w = 2, … , N are set to zero.
At Ω F and Ω B , only the tangential MVP component is tangential to those boundary surfaces and can therefore be affected by Dirichlet BCs, and only the boundary modes at the outermost SEs,ã t,1 e,1 andã t,N SE e,N +1 , need to be manipulated since they alone determine the solution at Ω F and Ω B . The Dirichlet data is incorporated through a system reduction [31], where the final system Ax = b is partitioned into known and unknown coefficients, where the subscripts u and k stand for "unknown" and "known", respectively. As the known coefficients x k do not need to be calculated, the corresponding system rows are discarded and the known term is moved to the right-hand side, yielding the constrained system As usual, homogeneous Neumann BCs being natural BCs do not need to be enforced explicitely.

GAUGING
The curl-curl operator in (3) leads to a gauge invariance of the MVP making it non-unique, as one can add an arbitrary gradient field in performing the gauge transformation ⃗ A * → ⃗ A + ∇Λ while the physically observable magnetic field does not change. This is reflected in the discrete system (35) by the equations being indefinite. To alleviate this problem, one can either add an explicit gauge condition to (35), or one can use an iterative solver whose properties implicitly gauge the system [32]. In terms of computational time and precision, a direct solver is to be preferred over an iterative solver. Therefore, an explicit Coulomb gauge following [21] is employed in this work, enforcing ∇ ⋅ ⃗ A = 0. This yields a saddle-point system where M FE,n is the nodal FE mass matrix The vector˜collects the Lagrange multiplier coefficients which are set to zero on the domain's boundary [21].

VERIFICATION AND CONVERGENCE STUDY
To verify the Q3D method, two studies are carried out: First, the consistency of the method is checked. Consider a box domain of pure vacuum, in which a constant magnetic field with an arbitrary orientation is imposed. The corresponding MVP is linear in space; therefore, the Q3D discretization should be able to represent the MVP exactly down to machine accuracy, independently of the number of elements or the maximal polynomial order and also independently of the field orientation. Indeed, when discretizing the analytical MVP solution and computing the numerical magnetic energy, the difference between numerical and analytical energy is in the order 10 −15 , so in range of machine accuracy. Thus, the Q3D method proves to be consistent. Next, a convergence study is carried out. In the Q3D setting, there are multiple ways to gain a precise discretization: Either one refines the transversal or longitudinal mesh, thus increasing the number of prism elements and edges, or one increases the maximal polynomial order N in longitudinal direction. These approaches are called hand p-refinement, respectively [33]. Here, the MS formulation (7) is solved in a cubic domain V = [−1, 1] 3 with appropriate boundary conditions and right-hand side for the manufactured solution for varying mesh discretizations and polynomial order sets. For each refinement, the numerical magnetic energy is compared to the analytical one. The resulting relative errors are plotted against the characteristic mesh length h (h-refinement) and the maximal polynomial order N (p-refinement) in Figures 6 and  7, respectively. As expected from theory, the energy error converges quadratically with decreasing mesh length h, that is, as the mesh becomes finer. Meanwhile, the energy error decreases rapidly with growing polynomial order N , indicating the exponential convergence known from spectral methods [14]. As a consequence, the p-refinement quickly reaches a point in which the error due to the mesh discretization gets dominant causing the total discretization error to stagnate with further increasing N (c.f. Figure 7). As a recommendation for the computational engineer, it follows that a small N is already sufficient and should be chosen Convergence study for the Q3D method in terms of a p-refinement. The relative energy error is plotted for different numbers of FEs against the maximal polynomial order N for the sake of computational efficiency, as the mesh discretization error outweighs the comparatively small improvements that would be generated by a p-refinement. Alternatively, one could use second-order FE ansatz functions to achieve a better convergence in the h-refinement.

APPLICATION AND RESULTS
As application example, consider a small superconducting cable with four multi-filamentary strands as depicted in Figure 8. The time-varying magnetic field excitation is applied on the domain boundary and a time domain simulation is carried out. A transport current is not present. A transient simulation for one whole period T = f −1 is performed, dur- Geometry of the benchmark model: Small superconducting cable with four multi-filamentary strands (blue) in vacuum (orange). In the application scenario, the south-eastern strand is subject to a quench in the region Ω Q × K Q ing which the strand segment Ω Q × K Q is treated as quenched, which is realized by setting with SC ≠ 0. So, the time constant is zero in the normalconducting matrix Ω NC × (K SC ∪ K Q ) and especially in the quenched segment Ω Q × K Q . The system (35) is time-stepped using a backward implicit Euler scheme. In Figure 10, the magnetic field flux magnitude resulting from the Q3D simulation is shown at the time instance in which (42) reaches its maximum. As expected, the magnetic field is expelled from the superconducting strands while penetrating the vacuum and the quenched strand segment. The magnetic energy W mag of the system is calculated for each time step. The results of the Q3D method are compared with a highly resolved reference solution computed with a 3D FE method using the software GetDP [34]. There, the full 3D geometry is discretized with 480 164 tetrahedra resulting into 614 062 mesh edges coinciding with the degrees of freedom (DoFs). For the Q3D simulation, a transversal mesh with 2132 triangles is chosen leading to N n = 1125 mesh nodes and N e = 3256 mesh edges representing longitudinal and transversal edges, respectively. A longitudinal discretization with N SE = 3 SE and maximal polynomial order N = 6 is considered. Altogether, the total number of Q3D DoFs results into (N e + N n )(N SE N + 1) = 83 239, which is only a fraction of the number of 3D DoFs. The magnetic energy over time is plotted in Figure 9 at the top. The Q3D result aligns very well with the result of the GetDP simulation. At the bottom in Figure 9, the relative difference for each time step between the Q3D and 3D simulation is depicted. The relative difference is approximately in the order of 10 −2 with exception of the time range in which the energy approaches zero, for  which rounding errors are most probably to blame. Thus, the Q3D method achieves here an excellent approximation of the 3D behaviour with a much higher computational efficiency.

CONCLUSION AND OUTLOOK
In this work, a Q3D FE method for the simulation of magnetic problems-especially for the quench simulation in superconducting magnets-has been developed using standard linear 2D edge functions in the xy-plane and modified Lobatto polynomials in z-direction. The Q3D matrices have been derived and the insertion of boundary and gauging conditions has been explained. The method has been verificated and its convergence has been shown. Eventually, the method has been applied to simulate a quench scenario within a superconducting cable model. A comparison to a conventional full 3D FE method has demonstrated the superior computational efficiency of the Q3D method.
In future work, the Q3D method will be used to simulate the coupled magneto-thermal quench propagation problem, where the heat conduction is treated with a combination of nodal FE functions and modified Lobatto polynomials. For this task, several additional aspects have to be considered to obtain an accurate prediction of the quench propagation such as field-circuit coupling [35] and interstrand coupling currents [9], which cause a snapback effect during the injection phase [36].
Compared to a conventional 3D FE method, a superior computational efficiency is expected from the Q3D method while sufficiently being able to resolve the quench phenomena in the superconducting cable.