Quasi-3D Thermal Simulation of Quench Propagation in Superconducting Magnets

09/18/2019 ∙ by Laura A. M. D'Angelo, et al. ∙ Technische Universität Darmstadt 0

To deal with the multi-scale nature of the quench propagation problem in superconducting magnets, this work presents a quasi-three-dimensional (Q3D) approach combining a two-dimensional finite-element method (FEM) in the transversal cross-section of the magnet for resolving the geometrical details, with a one-dimensional spectral-element method based on orthogonal polynomials in longitudinal direction for accurately and efficiently representing the quench phenomena. The Q3D formulation is elaborated and the idea is illustrated on a thermal benchmark problem. Finally, the method is validated against a conventional 3D simulation carried out by a commercial software. In terms of computational efficiency, it is shown that the proposed Q3D approach is superior to the conventional 3D FEM.



There are no comments yet.


page 1

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 co-simulation framework for superconducting accelerator magnets called STEAM [7]. A major numerical challenge arises due to the multi-scale 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 cross-section contains many geometrical details (see Fig. 1a). Additionally, quenches are very local effects leading to steep temperature fronts in longitudinal direction. Therefore, conventional three-dimensional (3D) finite-element (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 magneto-thermal coupled problem, altogether leading to time-consuming simulations, e.g. demonstrated in [9] and recently in [10].

Figure 1: In (a), the cross-section of the main dipole magnet in the LHC is shown. The coils are depicted in blue and red depending on the current direction, the iron yoke is colored in gray. In (b), the shell-type configuration of the coils is illustrated in detail showing the rectangular Rutherford cables and the trapezoidal copper keystones. The benchmark model (c) represents a stack of three Rutherford cables (white) wrapped by glass fibre insulation (black).

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 quasi-3D (Q3D) ansatz. Herein, the transversal cross-section of the magnet is discretized with a two-dimensional (2D) FE method, while the longitudinal direction is resolved with a one-dimensional (1D) spectral-element (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 higher-order 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
Table I: Benchmark model characteristics.

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 cross-section lies in the -plane.

In the investigated scenario, the left cable quenches at a location m. The resulting heat excitation is described by a time-constant Gaussian function


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 boundary-value problem


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 steady-state solution, one has to solve (2)-(4) with .

Iii Quasi-3D Formulation

Iii-a Spatial discretization and basis functions

Figure 2: In (a), the triangular FE mesh in the transversal cross-section (green) as well as the 1D SE mesh in longitudinal direction (red) are shown. This discretization results into triangular prism elements (b).

In the Q3D setting, the transversal cross-section 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


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


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


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.

Iii-B Discrete system

Figure 3: Tensor-product structure of the Q3D matrices: On the top level, there are block matrices for each SE which consist of a mode-wise block matrix pattern. These lowest level block matrices are the nodal 2D FE matrices.

Applying the Ritz-Galerkin method to (2)-(4) and choosing the test functions identically to the shape functions, , the system of equations


arises. The Q3D matrices and vectors are constructed by the 2D Cartesian FE and 1D SE matrices and vectors by


Herein, denotes the Kronecker tensor product. Fig. 3 sketches this tensor-product structure of the Q3D matrices. The FE quantities are found as


with as the 2D cross-section and standing for and , respectively. The element-wise SE quantities read



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 element-wise SE matrices diagonally like illustrated in Fig. 3 on the left, where it is important to overlap the element-wise matrices’ corners to ensure -continuity across the SE interfaces [19]. The semi-discrete 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.

Iii-C Adaptive spectral meshing

Figure 4: Temperature distribution in the quenched cable in -direction for two different points in time. The dashed lines indicate the SE interfaces and their adaptation to the moving quench fronts.

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

Figure 5: Temperature distribution in the benchmark model after ms. The geometry has been sliced at the hot-spot level .

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.

Figure 6: Temperature curve over time at for the quenched cable and its neighboring cable. The reference solution by COMSOL is depicted as black dashed line. The red curve shows the Q3D solution without making a spectral mesh adaption, while the blue curve depicts the Q3D solution after adapting the spectral mesh after half-time of the simulation (indicated with a blue dot).

Fig. 6 shows the resulting temperature curves for the hot-spot 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 half-time of the simulation to the new quench front locations, the Q3D solution (blue) is slightly improved and delivers relative differences below .

Figure 7: The dimension of the system matrix in the Q3D (red) and 3D COMSOL (black) case are plotted against the relative error w.r.t. a fine reference solution computed by COMSOL. The numbers below the marks denote the number of triangular elements resolving the transversal cross-section. The reference solution itself is resolved with 252 triangles.

To compare the spatial accuracy of the Q3D and 3D method, steady-state 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 hot-spot temperature w.r.t. a fine COMSOL reference solution with 252 triangular elements in the cross-section. Clearly, the Q3D method yields more accurate results with less computational effort than the 3D simulation, e.g. for a cross-section 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 FE-SE 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.


  • [1] S. Russenschuck, Field Computation for Accelerator Magnets. Weinheim, Germany: Wiley-VCH, 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 Co-Simulation 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., "3-D Thermal-Electric 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, "Three-dimensional 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, "Quasi-3D Finite-Element Method for Simulating Cylindrical Induction-Heating Devices," IEEE J. Multiscale and Multiphys. Comput. Techn., vol. 2, August 2017.
  • [13] L. A. M. D’Angelo, Y. Späck-Leigsnering, and H. De Gersem, "Electroquasistatic quasi-3D 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. Berlin-Heidelberg, 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. Fakhar-Izadi, and M. Dehghan, "A spectral element method using the modal basis and its application in solving second-order 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.