I Introduction
Superconducting accelerator magnets are used in the Large Hadron Collider (LHC) at CERN to achieve high magnetic fields while keeping the energy consumption low [1]. This comes at the cost of having to deal with quenches [2], which can lead to extensive damages of the magnet, e.g. reported in [3], making a reliable quench protection system mandatory. Quench protection systems are triggered upon quench detection, i.e. when predefined voltage or resistance thresholds are exceeded [4]. However, today, these thresholds are equipped with a large safety margin resulting into an overly sensitive quench detection system leading to frequent and unnecessary shutdowns, which considerably reduce the availability of the whole LHC machine due to false triggers [5]. This issue becomes even more critical for the newly developed magnets based on the NbSn technology [6].
Numerical simulations play a crucial role for understanding quench phenomena and therefore for improving quench detection. Yet, computational engineers struggle with the problem’s complexity, which is the reason for the recent development of a hierarchical cosimulation framework for superconducting accelerator magnets called STEAM [7]. A major numerical challenge arises due to the multiscale nature of both the magnet’s geometry and the transient effects of a quench: While a superconducting accelerator magnet is over m long, its diameter only measures a fraction of that [8], and the crosssection contains many geometrical details (see Fig. 1a). Additionally, quenches are very local effects leading to steep temperature fronts in longitudinal direction. Therefore, conventional threedimensional (3D) finiteelement (FE) simulations are computationally too expensive, as they would need very fine meshes to resolve the geometry and physical phenomena sufficiently, resulting into extremely huge systems of equations. This is further aggravated by the fact that it is a nonlinear magnetothermal coupled problem, altogether leading to timeconsuming simulations, e.g. demonstrated in [9] and recently in [10].
In the context of STEAM [7], this work presents an alternative approach by exploiting the geometrical translational invariance in longitudinal direction, which allows to use a quasi3D (Q3D) ansatz. Herein, the transversal crosssection of the magnet is discretized with a twodimensional (2D) FE method, while the longitudinal direction is resolved with a onedimensional (1D) spectralelement (SE) method using orthogonal polynomials. In this way, the system’s size can be significantly reduced. Similar approaches have been made e.g. in [11] for a geometry with translation symmetry using higherorder FE methods, and in [12][13] for cylindrical structures using a Fourier expansion in azimuthal direction.
Here, the idea is illustrated for a thermal benchmark problem with linear material characteristics. The benchmark model represents a stack of Rutherford cables used for the coil in the superconducting magnet [14]. The method is validated against a 3D simulation using the commercial software COMSOL [15] and their computational efficiency is compared.
This work is structured as follows. First, the benchmark model is introduced and the problem setup is explained in Sec. II. In Sec. III, the Q3D formulation for the heat conduction equation is elaborated. In Sec. IV, the simulation results are shown and discussed. Lastly, a conclusion of this work is given in Sec. V.
Ii Benchmark Model
Geometrical dimensions  Thermal conductivity (W/m/K)  

Cable width:  mm  Homogenized cable:  235.6 
Cable height:  mm  Glass fibre:  0.01 
Insulation thickness:  m  
Total width:  mm  Volumetric heat capacity (J/m/K)  
Total height:  mm  Homogenized cable:  314.1 
Model length :  m  Glass fibre:  750 
As benchmark model, a stack of three Rutherford cables wrapped by glass fibre insulation is considered, which is a component of the coil in superconducting magnets, see Fig. 1. Each Rutherford cable contains a number of wires consisting of superconducting NbSn filaments embedded in a copper matrix. These wires are not modeled explicitly. Instead, the cables are considered to be solid and the material properties are homogenized. Table I summarizes the geometrical dimensions and material characteristics of the benchmark model. Using a Cartesian coordinate system , the direction is defined as the longitudinal direction such that the transversal crosssection lies in the plane.
In the investigated scenario, the left cable quenches at a location m. The resulting heat excitation is described by a timeconstant Gaussian function
(1) 
where W/m is the excitation’s amplitude,
m is the standard deviation and
is a characteristic function that is equal to one for coordinates
in the left cable and zero otherwise. Furthermore, isothermal boundary conditions (BCs) with a fixed temperature K are set at the back and front sides of the model, and adiabatic BCs are applied to the hull . In the initial state, all cables are cooled down to K. Then, the heat conduction in the benchmark model is described by the boundaryvalue problem(2)  
(3)  
(4)  
(5) 
Herein, is the thermal conductivity in W/m/K, is the temperature in K, is the volumetric heat capacity in J/m/K, is the spatial coordinate in m, and is the time in s. For the thermal steadystate solution, one has to solve (2)(4) with .
Iii Quasi3D Formulation
Iiia Spatial discretization and basis functions
In the Q3D setting, the transversal crosssection is discretized with triangular FEs while the longitudinal direction is resolved with a 1D SE mesh, i.e. with line elements . As a result, the model’s volume is assembled from triangular prism elements, see Fig. 2. The temperature is discretized by
(6) 
where is the number of FE nodes, is the number of SEs and is the maximal polynomial degree. The coefficients
collect the degrees of freedom (DoFs). In contrast to the FE method, these coefficients do not entirely live in the
physical space, i.e. they do not represent the temperature solution directly, but they form the expansion coefficients for the polynomials and thus represent something like the frequency of the solution. Therefore, it is said that they live in the spectral or frequency space [16]. The shape functions are expressed as a product(7) 
of being the standard FE linear 2D nodal shape functions [17], and are chosen as modified Lobatto polynomials, which are defined on the orthogonality interval globally as
(8) 
with a mapped parameter and being the Lobatto polynomial of th order [18]. The choice of this particular polynomials is motivated by the desire to achieve sparse matrices: The FE nodal shape functions are chosen to be nodal, i.e. each node only interacts with its neighboring nodes. Similarly, the SE polynomial shape functions are chosen to be modal, i.e. their interactions in the frequency space are reduced [19]. An exception are the boundary modes , which are nodal and therefore represent the physical solution. In this way, although the SE matrices are not purely diagonal, the nodal boundary modes will make it easy to impose interface continuities and BCs. By using the Vandermonde matrix of the modified Lobatto polynomials [18], it is possible to forward transform a physical solution into a spectral one, or to backward transform a spectral solution into a physical one.
IiiB Discrete system
Applying the RitzGalerkin method to (2)(4) and choosing the test functions identically to the shape functions, , the system of equations
(9) 
arises. The Q3D matrices and vectors are constructed by the 2D Cartesian FE and 1D SE matrices and vectors by
(10)  
(11)  
(12) 
Herein, denotes the Kronecker tensor product. Fig. 3 sketches this tensorproduct structure of the Q3D matrices. The FE quantities are found as
(13)  
(14)  
(15) 
with as the 2D crosssection and standing for and , respectively. The elementwise SE quantities read
(16)  
(17)  
(18) 
where
is the modified Lobatto interpolation operator. Furthermore, it is assumed that the excitation function can be separated according to
. The global SE matrices are built by ordering the elementwise SE matrices diagonally like illustrated in Fig. 3 on the left, where it is important to overlap the elementwise matrices’ corners to ensure continuity across the SE interfaces [19]. The semidiscrete system (9) is discretized in time with an implicit backward Euler method. Eventually, the physical solution is obtained by a backward transform of the spectral solution.IiiC Adaptive spectral meshing
By choosing the SE interfaces appropriately, only a few polynomial orders suffice to resolve the quench phenomenon in direction. However, as the quench fronts move over time it becomes necessary to adapt the SE interfaces to them, because otherwise more polynomials would be required to represent the temperature distribution well, see Fig. 4. Even worse, steep temperature fronts within a SE may produce the Gibbs phenomenon or spectral ringing [20]. Therefore, an adapted spectral meshing of the longitudinal direction is mandatory to achieve accurate results with affordable computational costs.
This is done by regularly checking the solution during the simulation and shifting the old SE interfaces (Fig. 4, dashed black lines) according to the quench fronts (red dashed lines) when necessary. The physical solution is then forward transformed onto the new spectral mesh, and the matrices are reassembled. In this process, the structure and size of the matrices are maintained, as only the mapping coefficients need to be recalculated.
Iv Simulation Results
The benchmark problem described in Sec. II is computed with the proposed Q3D method and validated against a conventional 3D simulation using COMSOL [15]. The temperature distribution after ms is visualized in Fig. 5, where the model has been sliced at enabling to see the transversal temperature distribution. As a consequence of the glass fibre insulation, the heat exchange is faster in longitudinal than in transversal direction. As a result, there are sharp temperature gradients from one cable to another, while the temperature varies smoothly in longitudinal direction.
Fig. 6 shows the resulting temperature curves for the hotspot point at in the quenched cable as well as in the neighboring cable. Here, the Q3D solution (red) aligns well with the reference solution (black) leading to relative differences below at the end of the simulation time. When adapting the spectral mesh after halftime of the simulation to the new quench front locations, the Q3D solution (blue) is slightly improved and delivers relative differences below .
To compare the spatial accuracy of the Q3D and 3D method, steadystate simulations are carried out to exclude time stepping errors. The computational efficiency is visualized in Fig. 7, where the dimension of the system matrix is plotted against the relative error of the hotspot temperature w.r.t. a fine COMSOL reference solution with 252 triangular elements in the crosssection. Clearly, the Q3D method yields more accurate results with less computational effort than the 3D simulation, e.g. for a crosssection resolution with 65 triangles, the 3D simulation requires a system matrix of size , while the Q3D method achieves the same transversal resolution and a similar error with a smaller system matrix of size . It must be noted that the benchmark model has a length of m, while the real geometry of a superconducting accelerator magnet is over m long. Consequently, it is expected that the efficiency gain obtained by the Q3D approach will be much more pronounced for the real magnet geometry.
V Conclusion
To deal with the numerical challenges arising in the simulation of quench propagation in superconducting magnets, a Q3D FESE method has been presented as an alternative to the conventional 3D FE method. As a proof of concept, the heat conduction in a Rutherford cable stack has been computed and successfully validated against a 3D FE simulation using COMSOL. In terms of computational efficiency, the Q3D approach has proved to be superior to the conventional 3D FE method. Therefore, the presented Q3D method is a promising and efficient numerical tool to accurately simulate quench propagation phenomena in superconducting accelerator magnets.
References
 [1] S. Russenschuck, Field Computation for Accelerator Magnets. Weinheim, Germany: WileyVCH, 2010.
 [2] M. N. Wilson, Superconducting Magnets. New York, USA: Oxford University Press, 1983.
 [3] M. Bajko et al., Report of the Task Force on the Incident of 19 September 2009 at the LHC. Geneva, Switzerland: CERN, 2009.
 [4] J. Steckert et al., "Application of the New Generic Quench Detection System for LHC’s 11 T Dipole Magnet," IEEE Trans. Appl. Supercond., vol. 29, no. 5, August 2019.
 [5] A. Apollonio, "Machine Protection: Availability for Particle Accelerators", Ph.D. dissertation, Technische Universität Wien, Wien, Austria, March 2015.
 [6] G. Ambrosio, "Nb3Sn High Field Magnets for the High Luminosity LHC Upgrade Project," IEEE Trans. Appl. Supercond., vol. 25, no. 3, June 2015.
 [7] L. Bortot et al., "STEAM: A Hierarchical CoSimulation Framework for Superconducting Accelerator Magnet Circuits," IEEE Trans. Appl. Supercond., vol. 28, no. 3, 2017.
 [8] O. S. Bruning et al., LHC Design Report Vol.1: The LHC Main Ring. Geneva, Switzerland: CERN, 2004.
 [9] S. Caspi et al., "Calculating Quench Propagation With ANSYS®," IEEE Trans. Appl. Supercond., vol. 13, no. 2, June 2003.
 [10] J. F. Troitino et al., "3D ThermalElectric Finite Element Model of a NbSn Coil During a Quench," IEEE Trans. Appl. Supercond., vol. 29, no. 5, August 2019.
 [11] H. V. Jorks, E. Gjonaj, T. Weiland, and O. Magdun, "Threedimensional simulations of an induction motor including eddy current effects in core lamination," IET Sci. Meas. Technol., vol. 6, no. 5, 2012.
 [12] L. A. M. D’Angelo, and H. De Gersem, "Quasi3D FiniteElement Method for Simulating Cylindrical InductionHeating Devices," IEEE J. Multiscale and Multiphys. Comput. Techn., vol. 2, August 2017.
 [13] L. A. M. D’Angelo, Y. SpäckLeigsnering, and H. De Gersem, "Electroquasistatic quasi3D finite element simulation of a graded surge arrester," Int. J. Numer. Model., 2019.
 [14] A. P. Verweij, "Electrodynamics of Superconducting Cables in Accelerator Magnets," Ph.D. dissertation, Universitaeit Twente, the Netherlands, 1995.
 [15] COMSOL Multiphysics® v. 5.3. Stockholm, Sweden: COMSOL AB.
 [16] J. Shen, T. Tang, and L. Wang, Spectral Methods – Algorithms, Analysis and Applications. BerlinHeidelberg, Germany: Springer, 2011.
 [17] S. Brenner, and R. Scott, The Mathematical Theory of Finite Element Methods. New York, USA: Springer, 2008.
 [18] C. Pozridikis, Introduction to Finite and Spectral Element Methods Using MATLAB. CRC Press, 2014.

[19]
F. FakharIzadi, and M. Dehghan, "A spectral element method using the modal basis and its application in solving secondorder nonlinear partial differential equations,"
Math. Methods Appl. Sci., vol. 38, Jan. 2015.  [20] J. P. Boyd, Chebyshev and Fourier Spectral Methods. USA: Dover Publications, 2001.
Comments
There are no comments yet.