Quasi-3-D Spectral Wavelet Method for a Thermal Quench Simulation

The finite element method is widely used in simulations of various fields. However, when considering domains whose extent differs strongly in different spatial directions a finite element simulation becomes computationally very expensive due to the large number of degrees of freedom. An example of such a domain are the cables inside of the magnets of particle accelerators. For translationally invariant domains, this work proposes a quasi-3-D method. Thereby, a 2-D finite element method with a nodal basis in the cross-section is combined with a spectral method with a wavelet basis in the longitudinal direction. Furthermore, a spectral method with a wavelet basis and an adaptive and time-dependent resolution is presented. All methods are verified. As an example the hot-spot propagation due to a quench in Rutherford cables is simulated successfully.


Introduction
A standard finite element method (FE method) allows to discretize and solve differential equations on arbitrary geometries [1]. Attaining a sufficient resolution with a standard three-dimensional (3-D) FE approach may be intractable within a limited computation time. For certain geometries, however, an adapted discretization and simulation method may reduce the computational expense. This work focuses on multi-scale problems featuring domains with a high ratio between the extend in one spatial direction, referred to as longitudinal direction, and the extend in the other spatial directions, referred to as cross-section, which are, moreover, the translationally invariant with respect to the longitudinal direction. If the cross-section requires a high resolution, because it is e.g. very detailed, performing a standard FE method would be computationally very expensive due to a large number of degrees of freedom. Usually, one would remedy this issue by reducing the model to a two-dimensional (2-D) [2] or even one-dimensional (1-D) [3] model by exploiting the symmetry planes. Unfortunately, this simplification is no longer possible as soon as the physical behavior breaks the symmetry, e.g. for a longitudinally asymmetrical heat source in a thermal heat conduction problem.
The described configuration occurs in superconducting magnets in particle accelerators. These magnets struggle with the quench phenomenon, which is a sudden and local transition from the superconducting to the normal state [4]. Quenches emerge randomly during continuous operation. The affected regions normally cool down again and return to a superconducting state in doing so. Sometimes, however,  a quench leads to a thermal runaway and the magnet can take irreversible damage if protection measures are not taken in time [4]. Sensitive safety margins lead to a high fault time of the accelerator [5]. The further improvement of the protection systems necessitates the study of the transient effects inside the magnet by numerical simulations [6].
An alternative to the FE approach is based on a quasi-3-D (Q3D) method, that becomes applicable by exploiting the geometrical symmetry of the model [7,8]. Here, a 2-D FE method in the cross-section is combined with a 1-D spectral method [9] in the longitudinal direction. While this approach utilizes FE ansatz functions of lowest order, the order of the spectral element ansatz functions is chosen arbitrarily or appropriately for the problem at hand. Hence, the Q3D method is similar to an anisotropic hp-method [10].
In earlier works, harmonic functions [7] and orthogonal polynomials [8] have been used as basis functions for the spectral method. In this work, however, wavelets and the corresponding scaling functions are used as an alternative to these. A comparison of orthogonal polynomials and wavelets is made. Wavelets are functions that are defined by properties within the multiresolution analysis (MRA) [11]. In this work, the Daubechies wavelets [12] are utilized. This paper is structured as follows. First, an introduction to wavelets is given and a spectral method using wavelets is derived. Then, the FE method and the spectral wavelet method are combined into a hybrid Q3D method. Next, an adaptive resolution strategy is discussed for the spectral method. Lastly, the methods are compared and employed to compute the heat conduction in a benchmark model representing a superconducting Rutherford cable as used in superconducting accelerator magnets [13,14], see Fig. 1. Here, the asymmetrical heat excitation generated by a quench is simulated by an artificial heat source. To obtain a realistic quench simulation, a magneto-thermal coupling is required, which would necessitate the use of vectorial edge shape functions [15]. These are more involved than nodal ones and beyond this work's scope [16].

General statements
The MRA [11] consists of closed spaces V j that are nested subspaces of each other, i.e. it holds that V j ⊂ V j−1 . Furthermore, the spaces W j := V ⊥ j ∩ V j−1 are defined as the orthogonal complement of V j in V j−1 . Thus, it holds that and By applying (1) to itself, one gets where ⊕ denotes the direct sum. There are scaling functions φ j,n and wavelets ψ j,n .
They are equipped with a scale j and a translation n and are defined as Within these definitions, the mother scaling function φ and the mother wavelet ψ are used. These coincide with the scaling function and wavelet, respectively, for j = 0 and n = 0. The scaling functions and wavelets at fixed scale j constitute an orthonormal basis to the space V j and W j , respectively. To be more precise, it holds that V j := span {φ j,n : n ∈ Z} , W j := span {ψ j,n : n ∈ Z} .

Daubechies-N kind wavelets
There is a great variety of different kinds of wavelets [12,17]. However, within the spectral methods of this work, scaling functions and wavelets of the DB-N kind are used, because they are compactly supported and constructed in a way to have a maximum number of vanishing moments [12]. To be more precise, the wavelets have a number of N vanishing moments, i.e. it holds x k ψ(x) dx = 0 for k = 0, . . . , N − 1. Moreover, they are constructed such that they have a minimal support width for the given number of vanishing moments [18]. The scaling functions used in this work are continuous for N > 1. In [19] one can find estimates for the Hölder regularity for the scaling functions with N = 2, . . . , 10. For N = 2, 3, 4, these are also determined in [12]. It follows that the scaling functions are continuously differentiable for N = 3, 4, 5 and two times continuously differentiable for N = 6, 7, 8. The wavelets have the same regularity properties. For N = 6, scaling function and wavelet are depicted in Fig. 2.
The differential equations to be solved in the subsequent sections are formulated on intervals. Consequently, for N ≥ 2 there are always functions that are partially inside and partially outside the interval. But a basis that is defined on the interval is required. The construction of such a basis as proposed in [20] is used in this work. Let I = [a, b] be an interval. For each boundary N boundary scaling functions and boundary wavelets are constructed, such that the number of vanishing moments is preserved. Inside the interval, the ordinary functions are used.
Because the boundary scaling functions and boundary wavelets can be written as a linear combination of scaling functions, they feature the same regularity. For a sufficient small scale j, the scaling functions on the interval I are defined bŷ The space of wavelets at scale j on the interval I, indicated by W I j , is defined analogously.

Spectral methods with wavelets-bases
As stated above, the Q3D methods is a combination of a 2-D FE method and a 1-D spectral method. The latter one, that uses a wavelet-basis, is examined in this section. It is exemplary developed for the 1-D heat equation. The partial differential equation (PDE) together with the boundary conditions and the initial condition reads as with (x, t) ∈ Ω t , Here, x is the spatial coordinate in m, t the time in s, θ(x, t) the unknown temperature in K, λ(x) the thermal conductivity in W/(m K), c V (x) the volumetric heat capacity in J/(m 3 K) and q(x, t) the volumetric heat flux density in W/m 3 . Furthermore, the problem is defined on the domain Ω t := Ω × (0, ∞), where Ω is the spatial domain. The boundary ∂Ω is divided in Γ dir and Γ neum , with ∂Ω = Γ dir ∪ Γ neum and Γ dir ∩ Γ neum = ∅, on which Dirichlet θ dir (x, t) in K and Neumann boundary conditions g neum (x, t) in W/m 2 are defined, respectively. The initial condition uses the initial value θ 0 (x) in K.
Let for this section be Ω = (0, L). The temperature θ(x, t) is then approximated byθ(x, t) in V Ω j for each time instance. The ansatz of method of lines leads to The Galerkin approach is applied to (2). Therefore, (2) is tested with the basis for m = 0, . . . , D − 1. After inserting (3) into (4), one receives the semi-discrete equation The used components are defined as Integration by parts is used in the calculation of the entries of A sf λ . For homogeneous material properties, M sf c V is a diagonal matrix due to the orthonormality of the scaling functions. Consequently, M sf c V is a regular matrix and (5) is a system of ordinary differential equations.
For the boundary conditions to be satisfied, one can find a matrix B sf ∈ R 2×D and a vector b sf ∈ R 2×1 such that the equation has to hold for every time instance. For discretizing (5) in time, an arbitrary time integration method can be used.

Definition of the domain and spaces
The Q3D method shall be derived by means of the 3-D heat equation where (x, t) ∈ Ω t . With a 2-D and open domain Ω S ⊂ R 2 and I = (0, L), the spatial domain defined as extends in z-direction from 0 to L. Moreover, it it assumed that the geometry is translationally invariant with respect to the z-direction. The boundary ∂Ω is divided in three different parts Next the related spaces need to be defined. Let be the space of polynomial finite elements of order m, where Π m is the set of all polynomials up to order m and T is a two-dimensional triangulation (see [1]). With this, Combining the bases of both spaces, S 1 and V I j , a basis for K j is given by with ν m ∈ S 1 ,φ j,n−1 ∈ V I j and I being the number of nodes in the triangulation T . To make the following derivation clearer, the notation is used to reduce the number of indices from 3 to 2. Thus, the basis for K j can be written as

Derivation of the Q3D method
The Q3D method is derived for the 3-D heat equation (6). The temperature is approximated in K j . Using the ansatz of method of lines leads to The heat equation is tested with κ j,m for m = 1, . . . , DI. Together with (9), one receives the semi-discrete equation similar to (5). It is assumed that the material parameters can be written as . Then, the matrices can be written as where ⊗ denotes the Kronecker product. The remaining matrices, i.e. the stiffness matrix of the 2-D FE method and mass matrices of the spectral and FE method, are defined as For homogeneous material properties, M fe c V is symmetric and positive-definite because of is non-singular due to the properties of the Kronecker product and thus, (10) is a system of ordinary differential equations.
Under the assumption that the volumetric heat flux density can be written as q(x, t) = q xy (x, y, t)q z (z, t), the right hand side of (10) reads as For the realization of the boundary conditions, the three different parts of the boundary, defined in (8), are considered individually. At the front and back side, S front and S back , this is done with the spectral method, and at the hull S side with the FE method.

Adaptive resolution
In the spectral method in Sec. 3, scaling functions at the same scale are used for the whole interval as basis functions. Thus, the resolution is constant for the whole interval as well. When considering a very long domain in z-direction, a high resolution becomes computationally very expensive due to the large number of degrees of freedom. With the adaptive resolution method (ARM) derived in this section, a resolution depending on the location is chosen utilizing the wavelet transform, such that the overall accuracy is sufficiently high. This can strongly decrease the required degrees of freedom.
The boundary basis functions (Sec. 2.2) are chosen such that the main properties of wavelets still hold. In particular, it holds that This allows to decompose the spaces V I j . Fig. 3 shows an example of such a decomposition. The space V I 0 is decomposed in V I 1 and W I 1 , V I 1 in V I 2 and W I 2 , and so on. Approximating a function in V I 0 is thus equivalent to approximating the same function in the boxed spaces and adding the solutions.
Let j max and j min denote the maximum and minimum scale used in the approximation, respectively. The approximation is then equivalent to one with scaling functions at scale j = j min − 1. The basis for the adaptive resolution is It contains all basis functions of V I jmax (first set) to ensure that every point in the interval I is covered by at least one function. Moreover, it contains subsets of the bases of the spaces W I j for j = j min , . . . , j max (second set). The information which of the latter are used is stored in the η j -tuples K j = K j,1 , . . . , K j,ηj . For these, K j,i < K j,i for i < i and 1 ≤ K j,i ≤ D j for all i. How to choose the tuples is discussed later on page 10.
The tuples are also used to define the truncation matrices where e Kj,n ∈ R Dj is the K j,n -th unit column vector. The basis Ξ jmax jmin is used within a spectral method for the 1-D heat equation from (2). This yields the semi-discrete Herein, the used vectors and matrices as well as B d , the matrix used for the boundary conditions, are determined with the help of A s , M s , B s , u s (t) and q s (t). The subscript s stands for a static resolution, whereas d stands for a dynamic resolution. One would receive the static matrices and vectors if all wavelet basis functions in Ξ jmax jmin were used, i.e. with K j = (1, . . . , D j ) for j = j min , . . . , j max . Because of (11), these matrices and vectors can be derived from those from Sec. 3 at scale j min − 1. With the transformation matrix where I Dj max ∈ R Dj max ×Dj max is the identity matrix and η := D jmax + jmax j=jmin η j , one obtains When changing the resolution, the matrices A d and M d have to be recalculated. Because A s and M s are known beforehand and T jmax jmin is sparse, the computational costs for recalculation are low.
It remains to determine the tuples K j , where · , · := |(· , ·)| is used. Note that · , · is not an inner product but an abbreviation. The initial condition θ 0 (x) is used at the beginning of the simulation:  The basis functionψ j,k is used if the corresponding coefficient is greater than a tolerance τ > 0, i.e. if the basis functionψ j,k has a significant contribution to the solution.
During the simulation it may occur that the resolution needs to be adapted. One has to distinguish between decreasing and increasing the resolution (see Fig. 4). For the former case, after every time step those basis functions whose coefficients are smaller than τ are removed from Ξ jmax jmin . For the latter case, two different causes are considered. An increase of the resolution may be necessary due to the source term q(x, t) or heat flow. Because the source term is known beforehand, the required resolution can be determined as for the initial condition. The functionψ j,k is used if q(x, t) ,ψ j,k > τ . Because the solution of the PDE is approximated with the basis Ξ jmax jmin , it is challenging to detect where a higher resolution caused by heat flow is needed. One possible approach is presented in the following.
Whenever a coefficient in u d (t) changes from one time step to another by a certain factor α, the resolution in the interval of the support of the corresponding basis functionψ j,k is increased. Therefore, the wavelets at scale j − 1, whose support is a subset of the support ofψ j,k , are added to Ξ jmax jmin , if j − 1 does not exceed the minimum scale j min , i.e. if j > j min . Assuming that j > j min , the scaling functions at scale j − 1 with the indices Λ := {−N + 1 + 2k, . . . , N + 2k} ∩ 0, 1, . . . , L2 −j+1 − 1 are added to Ξ jmax jmin , i.e.k + 1 ∈ K j−1 for allk ∈ Λ. The corresponding coefficients to these wavelets are initialized with 0. Hence, the added wavelets must not be removed from Ξ jmax jmin for a time interval of length β because the coefficients may need a certain time to exceed the tolerance τ . This approach has been implemented and is used for the results in Sec. 7.3.

Adaptive Quasi-3-D method
The ARM from the previous section is combined with a 2-D FE method into an adaptive Q3D (AQ3D) method. The adaptivity only takes place in the longitudinal direction, i.e. the 2-D mesh in the cross-section does not change. The method is derived for the 3-D heat equation (6) and on the domain Ω defined in (7). The space for basis and test functions is In contrast to the Q3D method, here the function g is chosen from the space Ξ jmax jmin . With {χ n : n = 1, . . . , η} being a basis of Ξ jmax jmin , a basis for K jmax jmin is given by with ν m ∈ S 1 and χ n ∈ Ξ jmax jmin . Similar to (10), the semi-discrete equation is Herein, the matrices can be written as With (12) and (13)  where I I ∈ R I×I is the identity matrix, the matrices can be written as The middle term is independent of the resolution and has to be calculated only once at the beginning of the simulation. Similar to the matrices, the right-hand side vector is The boundary conditions are realized at S front and S back with the ARM and at S side with the FE method. The process of changing the resolution is very similar to the ARM. Here, a separate resolution is calculated for every edge in z-direction in the same way as explained in Sec. 5. The resolutions of all edges and the source term are then merged with each other to get a single resolution for all edges. Consequently, the matrix Q jmax jmin can be recalculated and the simulation can proceed.

Simulation results
The time discretization for each method is realized with the implicit Euler method, where ∆t denotes the step width. The difference between numerical and analytical solution is measured with where t i = i∆t is the time step. The norm · max is used to compareθ with an analytical solution as a function, whereas · ∞ compares the vector u with a vector u ana that contains the coefficients of the analytical solution.

Spectral scaling function method
For the spectral scaling function method (SSM), the 1-D heat equation (2) with dimensionless L = 10, N = 6, λ = 10, c V = 1, ∆t = 1 · 10 −3 and 100 time steps is solved. The initial condition is θ(x, 0) = sin π L x , the source term is q(x, t) = 0, and homogeneous Dirichlet boundary conditions are applied. The error between numerical and analytical solution is shown in Fig. 5. As expected, the error in both norms decreases with decreasing scale.

Quasi-3-D method
The Q3D method is used to simulate the temperature distribution for the benchmark problem defined in [8] and highlighted in red in Fig. 1. The domain consists of three insulated cables. They extend from z = 0 m to z = 1 m and are surrounded and separated by an insulation material (glass fiber). The problem's spatial dimensions  and material data are summarized in table 1 and 2, respectively. In the model, the cables have been homogenized into a bulk model consisting of one material, although in reality they are built up of several strands [21]. An initial condition of θ 0 (x) = 2 K together with Dirichlet boundary conditions at the front (z = 0 m) and back side (z = 1 m) and homogeneous Neumann boundary conditions are applied at the hull.
In order to simulate a quench in the left cable, the time independent source term is applied, where χ q (x) is equal to 1 in the left cable and equal to 0 elsewhere. The remaining values are chosen as q max = 1 · 10 6 W/m 3 , z q = 0.33 m and σ = 0.05 m. For the simulation, scaling functions with N = 6 at scale −5 are used and 200 steps of the implicit Euler method with ∆t = 5 · 10 −5 s are performed. This leads to a total simulation time of 10 ms.
In Fig. 6, the evolution of the temperature inside the three cables is shown. For this plot and for Fig. 7, the temperature was evaluated at z = z q and for each cable in the center of its cross-section. The temperature of the left cable increases immediately, while the other cables start to heat up after a certain time. At the end of the simulation (at t = 10 ms) the left cable has thus the highest temperature. This can be seen in Fig. 7 and Fig. 8 as well. In Fig. 7 one can additionally see how the heat has already been flown from z = 0.33 m in positive and negative z-direction. Fig. 8 shows the temperature distribution at the end of the simulation at z = z q and displays the mesh in the cross-section. One can see that the temperature inside the three cables is nearly constant, whereas in the insulation material high temperature gradients are present. This temperature distribution is the result of the strongly differing values of the thermal conductivity of the cable and insulation (table 2) and of the source term in (14), that is independent of x and y inside the left cable.
The Q3D method with scaling functions (denoted with Q3D sf ) is now compared to a Q3D poly method (from [8]) and to a conventional 3-D FE method. The Q3D poly method uses a polynomial basis with a maximum degree M for the spectral method in longitudinal direction. It decomposes the longitudinal direction in non-equidistant intervals called elements. The problem that is solved can be found in (6) The Q3D sf method is used with N = 6 and the Q3D poly method is used with a maximum polynomial degree M = 4 on each element in z-direction. The FE method uses a tetrahedral mesh with linear, nodal basis functions.
In Fig. 10, the solution of the ARM is compared to those of the SSM. The error of the former method at a minimum scale j min has to be compared to the error of the latter at scale j = j min − 1. Both methods yield the same error, as expected. However, the ARM uses less basis functions than the SSM. Fig. 11 shows the number of used basis functions and the maximum possible number of basis functions for every time step and for the different minimum scale. As one can see, the actual number is below the maximum number of L2 −jmin+1 .

Adaptive Q3D method
For verification, the AQ3D method is used to solve the problem in (6) with dimensionless L = 10, λ = 10, c V = 5, ∆t = 1 · 10 −4 and 10 time steps. The initial condition is taken from (15) and the source term is q(x, t) = 0. Homogeneous Dirichlet boundary conditions are applied on the front and back side, while homogeneous Neumann boundary conditions are applied to the hull. Fig. 12 shows the maximum error between numerical and analytical solution over all edges in z-direction. The error on each edge is calculated with the infinity norm. The number of basis functions is increased by using both finer 2-D meshes in the cross-section and a smaller minimum scale j min in longitudinal direction. Furthermore, the maximum scale is set to j max = −1 and for the spectral method, scaling functions and wavelets with N = 6 are utilized. Lastly, to show the potential reduction of basis functions with the AQ3D method, it is used to solve the benchmark problem with three Rutherford cables defined in [8]. Instead of the cables being extended 1 m in z-direction, they are extended to 10 m for this simulation. Because the center of the source term is still z q = 0.33 m, there are no major differences expected in the temperature distribution between z = 0 m and z = 1 m. The remaining parameters are the same as in Sec. 7.2. Moreover, scaling functions and wavelets with N = 6, j min = −4 and j max = −3 have been chosen. Thus, the results are comparable with the results from Sec. 7.2. The tolerance is τ = 1 · 10 −4 . With the indicated values, the maximum number of basis functions in longitudinal direction is 320. Furthermore, the used mesh in the cross section has 1289 nodes and 2496 elements. This leads to a maximum number of 412 480 basis functions for the AQ3D method.
In Fig. 13, the evolution of the temperature in the three cable is shown for the Q3D and the AQ3D method. Both methods lead to the same evolution. However, the number of basis functions of the AQ3D method is significantly lower than the maximum possible number of basis functions of 412 480 as depicted in Fig. 14. A simulation of this problem with the Q3D method with the same accuracy in the z-direction and the same mesh in the cross-section would require this maximum number of 412 480 basis functions.

Conclusion
A spectral method using Daubechies scaling functions has been formulated and verified. This has been combined with a 2-D FE method with a nodal basis into a quasi-3-D method. With this method, a stack of three simplified Rutherford cables has been simulated successfully. A third method with an adaptive resolution in space and time has been formulated for the 1-D heat equation. This method has been extended to a quasi-3-D method in the same manner as the spectral method. A comparison between the adaptive and associated non-adaptive methods showed a significant reduction of the number of basis functions, preserving the accuracy.