Electromagnetic Simulation of No-Insulation Coils Using H – $\phi$ Thin Shell Approximation

When simulating no-insulation high-temperature superconducting pancake coils with the finite element (FE) method, the high aspect ratio of the thin turn-to-turn contact layer (T2TCL) leads to unfavorable meshes in these thin layers as manifested by a high number of degrees of freedom (DoF) or mesh elements of poor quality which decrease the accuracy of the simulation results. To mitigate this issue, we propose to collapse the T2TCL volume into a surface using a thin shell approximation (TSA) for three-dimensional FE analysis. A <inline-formula><tex-math notation="LaTeX">$\vec{H}-\phi$</tex-math></inline-formula> formulation is used and solves for the magnetic field strength <inline-formula><tex-math notation="LaTeX">$\vec{H}$</tex-math></inline-formula> in conducting domains and the magnetic scalar potential <inline-formula><tex-math notation="LaTeX">$\phi$</tex-math></inline-formula> in insulating domains. This formulation avoids spurious currents and reduces the number of DoF in insulating domains. Automatically created thick cuts are used to deal with multiply connected domains. Particular attention is paid to the interpretation of these cuts and the corresponding basis functions in the context of pancake coil geometries. The efficiency of the formulation facilitates the resolution of each turn. In this way, local phenomena such as quench can be captured in a straightforward way. The TSA formulation is verified by comparison against a reference model with volumetrically meshed T2TCL and is shown to be accurate and efficient, significantly reducing the solution time while reducing the effort for creating high-quality meshes. The TSA is implemented in an open-source FE framework and the source code is shared alongside this paper.


I. INTRODUCTION
H IGH-TEMPERATURE superconductors (HTS) are promising for overcoming the magnetic field limits of low-temperature superconductors [1], [2], [3]. However, the low normal zone propagation velocity of HTS [3] makes the detection and protection against quenches a challenge [4]. A no-insulation (NI) winding technique, i.e., without turn-to-turn Manuscript  electrical insulation, was proposed for HTS pancake coils in 2011 [5]. This technique allows for a current by-pass of local hot spots or defects in radial direction [6], [7], [8]. This has enabled the design of so-called self-protected NI coils with high thermal stability (see [9] for an overview). Despite the increased thermal stability of NI coils, quenches can still occur [8], [10] and appropriate analysis techniques are needed to investigate quench detection and protection in NI coils [11].
Most studies of charging and discharging processes as well as thermal analysis of NI coils are performed using equivalent circuit models [6], [8] (again, see [9] for an overview). More comprehensive quench simulations studying also mechanical effects have been carried out by coupling a circuit model with a two-dimensional (2D) magnetic field model, a one-dimensional (1D) heat transfer and a 1D mechanical model [11].
In 2020, a 2D axially symmetric finite element (FE) model of NI coils using a homogenization technique and an anisotropic resistivity has been proposed [12] and extended to model parallelconnected coils forming a magnet [13]. However, quenches are local effects and no-insulation coils commonly exhibit true three-dimensional (3D) geometries (see, e.g., [14]) making the quench problem in NI coils an intrinsically 3D problem requiring 3D FE models. Full 3D resolution models are difficult due to, e.g., the high aspect ratio of the superconducting layer of coated conductors (CC) [15] and the contact or insulation layers between the turns [16], as well as the current flow across the turnto-turn contact layer (T2TCL) [12]. In a classical FE method, the T2TCL volume has to be meshed leading to unfavorable meshes since a high number of degrees of freedom (DoF) is needed to ensure an accurate solution, in particular if the problem is highly-nonlinear as common in the case of HTS materials due to the power law.
To mitigate this issue while still resolving each turn of the NI coil, we propose to collapse the T2TCL volumes to surfaces using a thin shell approximation (TSA) as developed in [17]. The TSA thus alleviates the need to volumetrically mesh the thin layers. 1D Lagrange elements are used to discretize the field behavior across the virtual thickness of the T2TCL [15], [16], [17]. The TSA-based method is verified against reference models with volumetrically meshed TCT2L and is shown to yield accurate solutions with significantly reduced solution time. The code is made available such that it may serve as a community benchmark [18].
The choice of formulation is key to ensuring robustness and efficiency of FE-based simulations involving HTS [19], [20]. field strength H such as the H − φ formulation are preferred over formulations based on the magnetic vector potential A due to the behavior of the involved material properties (if there are no ferromagnetic pieces [21]). The H − φ formulation uses the magnetic field strength as variable in conducting domains and the magnetic scalar potential φ in insulating domains. The latter choice ensures that no spurious current appears in insulating domains without the need for artificially large resistivities as in a pure H formulation solved for the magnetic field strength everywhere [21], [22]. Furthermore, the H − φ formulation also leads to fewer degrees of freedom (DoF) in insulating domains when compared to a pure H formulation. For this reason, the H − φ formulation is used in this work. However, in multiply connected domains, it requires appropriate techniques to ensure uniqueness of the posed problem [23]. In this work, thick cuts [24] are used which can be automatically created and the basis functions related to these cuts allow an evident physical interpretation in the context of pancake coils.
In Section II, the reference FE formulation with volumetrically meshed T2TCL is discussed. It is extended to the TSA formulation in Section III. Section IV summarizes the simulation setup for the verification of the TSA method. The specifics of pancake coil geometries when using a H − φ formulation are highlighted, in particular concerning the basis functions related to the thick cuts. Section V discusses the results of a current-driven powering cycle simulation conducted with the TSA method and compares it against the volumetrically meshed T2TCL reference model in terms of accuracy and efficiency. The findings are summarized in Section VI.

II. H − φ FINITE ELEMENT FORMULATION
In order to simulate the current ramp of a NI coil, the magnetoquasistatic (MQS) equations [25, Section 5.18] are considered with E the electric field in V/m, B the magnetic flux density in T, H the magnetic field strength in A/m, J the current density in A/m 2 and t the time in s. Furthermore, the fields are linked by the material relations with the magnetic permeability μ in H/m and electrical resistivity ρ in Ωm.
The considered computational domain Ω is shown in Fig. 1. It consists of an insulating domain Ω i with σ i = ρ −1 i = 0 S/m and a conducting domain Ω c with σ c = ρ −1 c > 0 S/m. The latter is subdivided into the bare (homogenized) CC Ω c,b and the T2TCL Ω c,cl . The domain is bounded by ∂V = Γ with the unit vector n pointing outwards. We impose assuming that Γ is far enough from the pancake coil. The conducting part Ω c is connected to Γ via current leads in order to impose source currents I src or voltages U src .
In Ω c , edge basis functions are used to represent the magnetic field strength H, while scalar basis functions are used in Ω i to  discretize the scalar magnetic potential φ. In multiply connected domains Ω i , ∇ × H = 0 does not imply that H = ∇φ and since ∇ × ∇φ = 0 for all φ, the latter choice does not allow to represent net source currents. In order to overcome this problem, cuts are introduced in Ω i to make φ discontinuous along those cuts [23]. Cuts as described in [24] are employed whose corresponding basis functions introduce the discontinuity of φ with the associated coefficients representing the net source current.
Following [24], the magnetic field strength is discretized as with N e the edge basis function of edge e, N n the nodal basis function of node n and Ψ C j the edge cohomology basis function, which is dual to the closed loops 1 C j . Examples of these basis functions Ψ C j for the NI coil geometry are shown in Fig. 2.
Their coefficients I C j are equal to the integral of H along C j , i.e., to the current passing through the surface enclosed by C j . Fortunately, it is possible to choose the cohomology basis functions in a way that allows an evident physical interpretation of the coefficients I C j [24,Section 5]. This choice will be discussed in Section IV for our model problem. Furthermore, the thick cuts can be created automatically using an existing opensource implementation in Gmsh [24], [27], even for complicated geometries. Currents I i can be imposed in a strong sense by fixing the corresponding coefficient while voltages V i can only be imposed weakly [28], [29], allowing to simulate an operation of the coil in current-driven (or voltage-driven) mode and to access terminal voltages (currents) easily. Following a Galerkin scheme, the test functions H are chosen in the same space with coefficients h e , φ n and I C j . The corresponding weak H − φ formulation then reads: from an initial solution at time t = 0, find H ∈ H(curl, Ω) s.t.
where H(curl, Ω) is the space of square integrable functions with square integrable weak curl in Ω [28]. Furthermore, (·, ·) Ω denotes the volume integral in Ω of the scalar product of the two arguments, while ·, · Γ denotes the surface integral on Γ of the two arguments. The second integral is only computed in Ω c in contrast to a pure H formulation which uses a large, but finite resistivity causing the integral to be computed in Ω.

III. THIN SHELL APPROXIMATION
We propose to use a modification of the TSA in [17] by i) allowing the thin layer to be located in-between conducting materials and ii) removing the need for modifications on the mesh level as explained later in this section. The T2TCL volume Ω c,cl is collapsed into a surface Γ c,cl as shown in Fig. 3, which removes the need for a volumetric mesh representation of Ω c,cl .
Since Ω c,cl is conducting, the tangential magnetic field needs to be discontinuous in order to be able to represent surface current densities. To this end, two different fields H + and H − (as well as test functions H + and H − ) with support restricted to the positive and negative side of Γ c,cl , respectively, are introduced. They represent H on both sides of Γ c,cl .
The discontinuity of the tangential magnetic field strength can be introduced by two equivalent approaches: i) using a duplication of mesh DoF as in [17] or by ii) using dedicated FE basis functions with support restricted to one side of Γ c,cl [30]. In this work, we follow the latter approach as it allows to couple the thermal TSA of [16] and this work without needing two different meshes for the thermal and MQS system.
An interface condition on Γ c,cl is enforced by considering an additional surface contribution in the formulation (5), i.e., where [ H ] = H + | Γ c,cl − H − | Γ c,cl is the jump of H across Γ c,cl and w is the unit vector normal to Γ c,cl . The surface term is used for an internal discretization of the MQS problem in the internal domainΩ c,cl which represents the volumetric T2TCL Ω c,cl as shown in Fig. 4. To this end, a local coordinate system ( u, v, w) is used in the following with ( u, v) oriented along the tangential direction of Γ c,cl . To build the mesh for the internal FE discretization,Ω c,cl is split into N auxiliary layersΩ The latter formulation is still comprised of volume integrals in Ω (k) c,cl . In order to decompose them into integrals over Γ  c,cl is assumed as Herein, H j t denotes the tangential magnetic field strength on Γ (j) c,cl and enforcing H 0 t = H t,− and H N t = H t,+ links the internal and external problem (see also Fig. 4). Furthermore, Ξ j (w) denotes but higher order functions can be used, too. Using the same ansatz for the test function H leads to a decomposition of the internal problem (7) into surface integrals over Γ (k) c,cl and 1D FE matrices in [w k−1 , w k ], see [17] for the final formulation.

IV. MODEL SETUP AND IMPLEMENTATION
The H − φ formulation with and without the TSA discussed in the previous sections are implemented in the open-source FE framework GetDP 3.4 [31] using the open-source software Gmsh 4.10.5 [27] for meshing. In particular, the Python application programming interface of Gmsh [32] is used to create the pancake coil geometry and mesh. An implicit Euler scheme with adaptive time step [33] is used for time integration and the FE basis used is of lowest order although higher order basis functions could be used. Since the T2TCL consists only of one material, N = 1 is used for the TSA model. Furthermore, Gmsh is also used to automatically create the necessary cuts [24]. The created cuts are shown in Fig. 2, with cut Ψ C 1 located in-between the two current leads and Ψ C 2 in the bore of the NI coil. These cuts allow an evident physical interpretation of the coefficients I C j . The input source current flowing into the current lead is represented by I C 1 , i.e., I C 1 = I src and I C 2 is equal to the total current flowing in azimuthal direction, i.e., I C 2 = I θ,total . To this end, I C 1 is imposed since the source current is known beforehand, while I C 2 is an output of the simulation, since I θ,total is not known beforehand due to the radial currents flowing across the T2TCL. For an insulated pancake coil, I C 2 would be known as well due to the absence of radial currents through the T2TCL.
The coil and simulation parameters are summarized in Table I. The equivalent resistivity ρ c,b reads with f HTS the volume fraction of superconducting material in the CC, f NC = 1 − f HTS the fraction of non-superconducting material (here, a mix of copper, silver and Hastelloy ® ) and ρ NC the homogenized resistivity of the latter. The resistivity of the HTS is given by the power law, i.e., Herein, E c = 1 × 10 −4 V/m, λ the fraction of current flowing in the HTS layer and all other parameters as defined in Table I. As ρ HTS is itself a function of λ, a non-linear root-finding problem needs to be solved to find λ as for example detailed in [20, Section III.D]. The idea is to equate the local electric fields across the HTS and non-HTS layers, solve this equation for λ using the bisection method and compute ρ c,b . A quasi Newton-Raphson iteration is used for linearization of the non-linear system; the Jacobian of ρ c,b is approximated by neglecting the Jacobians of λ, ρ NC and J c . The average azimuthal and radial current are computed by Alternatively, they could be computed from the field solution by averaging the integrated radial (or azimuthal) component of the current density J. The source voltage between the current leads U src is computed as a post-processing quantity using global quantities following the approach of [28].

V. VERIFICATION: POWERING CYCLE SIMULATION
For verification, a powering cycle of the NI coil has been simulated with three models i) ref: a fine-meshed reference model (see Fig. 2), ii) vol: a coarser-meshed version of ref and iii) TSA: a coarsely meshed TSA model (see Fig. 9 for the TSA mesh). The first two have a volumetrically meshed T2TCL. The mesh sizes of the latter two models have been chosen as coarse as possible under the condition that absolute errors in comparison with the ref model 2 are below 5% of the maximum over time of the absolute value of the source voltage U src and central axial magnetic flux density B z,c . This choice allows to compare the TSA and vol model in terms of efficiency. Let us note that the high aspect ratio of mesh elements for coarsely meshed T2TCL leads to convergence issues for the iterative scheme needed to solve the non-linear system. Such problems have not been observed for the TSA model. This is the reason that a coarse mesh can still be used for the TSA while keeping an accurate solution, which is not the case for the vol model. Fig. 5 depicts the applied source current and the computed central magnetic flux density which shows the typical delay with respect to the source current due to the radial currents. Furthermore, the exponential decay of the field after turning the source current off can be seen. All models are in good agreement. Fig. 6 shows the average azimuthal and radial currents for all models. As expected, I θ shows the same qualitative behavior as   B z,c as it is the source of the central axial field component. Fig. 7 shows the source voltage which follows the same qualitative behavior as I r as the latter leads to the main resistive voltage contribution. As seen in both figures, all models are in good agreement which is also supported by the absolute errors with respect to the ref model shown in Fig. 8, which are below 5% of the maximum of the ref model. The magnetic flux density in the r − z cut plane can be seen in Fig. 9, which also shows the mesh of the TSA model.
Having established the accuracy of the TSA model, the number of DoF and run time of the simulation are shown in Table II. The TSA model has 2.2 times fewer DoF and 1.9 times shorter total solution time compared to the vol model. As discussed in [16], the more the T2TCL thickness approaches zero (for directly touching turns), the more the solution time is reduced. In such a case, the TSA becomes the only practical method expected to lead to a (correct) solution. Furthermore, the effect increases with the number of turns. Last, the construction of a high-quality   volumetric mesh of the T2TCL commonly requires user-defined mesh controls which are not needed for the TSA [16]. In summary, the TSA provides an accurate and efficient alternative.

VI. CONCLUSION
This paper discusses the H − φ TSA for T2TCL representation in the FE simulation of the current cycle of HTS NI coils in 3D. The physical interpretation of the required cuts has been highlighted. The formulation has been implemented in the open-source FE framework Gmsh/GetDP, which also allows for the automatic creation of the cuts. A reference implementation has been published alongside this paper. The numerical simulation results for a current-driven powering cycle of a model NI coil have been compared against a reference implementation with volumetrically meshed T2TCL showing excellent qualitative and quantitative agreement. The TSA leads to a significant reduction of number of DoF and run time required for desired accuracy. Its efficiency facilitates the mesh resolution of each turn which is important for accurately capturing local phenomena such as quench. Furthermore, the TSA drastically simplifies the creation of high-quality meshes.