Nonlinear Magneto-Quasistatic Simulation of Superconducting Tapes With a-ψ Algebraic Formulation

The study of superconducting tapes requires specialized formulations to account for the high aspect ratio of the tapes and the extremely high conductivity of the thin superconducting layer. This article presents a new formulation coupling the line integral of the magnetic vector potential <inline-formula> <tex-math notation="LaTeX">$a$ </tex-math></inline-formula> defined in the whole domain to the stream function <inline-formula> <tex-math notation="LaTeX">$\psi $ </tex-math></inline-formula> defined on the nodes of the superconducting layer. Results on two benchmarks show that the proposed <inline-formula> <tex-math notation="LaTeX">$a-\psi $ </tex-math></inline-formula> formulation is effective in solving problems where the currents are limited to thin layers also in case of a strong material nonlinearity.


I. INTRODUCTION
H IGH-TEMPERATURE superconducting (HTS) materials represent a new frontier in high-energy applications: operating at elevated temperatures, when compared with low-temperature superconductors, they are able to withstand stronger magnetic fields and guarantee high transport current with minimal loss of energy [1].In this context, ReBCO coated conductors are a promising choice which could fulfill both technical and economical issues.ReBCO refers to composite material based on rare-earth elements, available in thin tape form in a layered structure, enveloped by a copper support layer.In this work, Yttrium-based tapes, namely, YBCO tapes, are considered (Fig. 1).The high tolerance of HTS tapes to the tensile strength and compressive strain due to the presence of the hastelloy substrate allows the development of the conductor on round core (CORC) technology: the tape is wounded around a central copper core, guaranteeing mechanical and electrical isotropy [2], which makes them suitable for high magnetic field applications.Due to the high aspect ratio of the tape, between the extremely thin thickness (from 40 to 100 µm) and the width, they can suitably be represented with infinitely thin sheets [3].Instead, when 3-D models of the tapes are considered, objects are discretized with volume elements (usually tetrahedra or hexahedra), with two major drawbacks.
1) A large number of elements are required for the discretization, increasing the number of unknowns.
2) The elements are characterized by a large aspect ratio, affecting the efficiency of the numerical solution [4].This article presents a new formulation based on the cell method that couples a 2-D discretization of the supercon-  ducting layer with a 3-D discretization of the surrounding environment.Particular attention is devoted to the treatment of the electrical resistivity, to avoid singularities of the final system when it reaches the typical infinitesimal values of superconductors.

II. MATERIAL MODELING
The YBCO tape under study is a second-generation coated conductor made by different layers as shown in Fig. 1.The tape width is 4 mm, whereas the thicknesses of the different layers are reported in Table I along with the electrical conductivities at 77 K [5].The tape conductivity is homogenized considering that all the layers are connected in parallel.The linear conductivities, due to the copper (Cu) stabilizer, the silver (Ag) layer, and the Hastelloy (Ha) substrate, are averaged using their thicknesses δ as weights This value is then scaled to reference the conductivity to the thickness of the superconductive (sc) YBCO layer Finally, the linear contribution is averaged with the nonlinear conductivity of YBCO, assuming the conventional power law with constant critical current J c [6] where E 0 = 0.1 mV/m, J c = 2.85 × 10 10 A/m 2 , and n = 30.5 [6].By exploiting the electrical parallel between the superconducting YBCO layer and the equivalent linear materials, it is possible to apply the current divider rule.Due to the scaling (2), the equivalent cross section of the two layers is identical, so the current divider rule holds in terms of current densities where J is the total equivalent current density in the tape.Equation ( 4) can be solved for different values of J obtaining the homogenized characteristic σ = σ sc (J ) + σ ′ lin of the tape.

III. DOMAIN DISCRETIZATION AND INTEGRAL VARIABLES
Following the principles of the cell method [7], the domain under study is subdivided into two complementary regions: the conductive domain, with negligible thickness compared with other dimensions, represented by 2-D layers, and the surrounding 3-D nonconductive domain.The latter is discretized with a tetrahedral mesh G, while the former is discretized with a surface triangular mesh G s , created with the constraint that the faces of G s are also faces of G. From the tetrahedral mesh, a barycentric dual mesh G is derived [7], where dual nodes are located in correspondence of the tetrahedra barycenters, dual edges connect adjacent nodes and pass through the primal face mid-points, and dual faces hinge around primal edges and are created from a closed loop of dual edges, Fig. 2(a).Similarly, a surface dual mesh Gs is built starting from G s , as shown in Fig. 2(b).
A. a Formulation in the 3-D Domain In the 3-D subdomain, the electromagnetic variables associated with the spatial entities of G and G are [see also Fig. 2(a)]: the line integral of the magnetic vector potential a along primal edges, the magnetic flux b through primal faces, the magneto-motive force h along the dual edges, and the electric current ĩ through the dual faces (the bold notation is used for the vectors that collect these quantities).According to [8], the following equations hold: where C and C are the face-to-edge incidence matrices defined on G and G, respectively, and M ν is the constitutive reluctance matrix.Using the identity C T = C, the previous equations can be combined in Standard tree gauging can be applied to guarantee the uniqueness of the solution [9].

B. ψ Formulation in the 2-D Domain
On the spatial entities of the surface mesh G s and Gs , the following quantities are defined [see also Fig. 2(b)]: the stream function ψ on primal nodes, the electric current i through primal edges, the electric voltage ẽ along dual edges, and the magnetic flux b through dual faces.As detailed in [10], after defining the classical constitutive resistance matrix M ρ and using the edge-to-node G s and dual face-to-edge Cs surface incidence matrices, it is possible to write the current field equations In (10), the nonlinear Ohm's law is linearized using the fixed-point scheme with residual R and fixed point resistivity ρ 0 .Using the identity Cs = G T s , the previous equations are combined in

C. Coupling Terms
It is worth noting that the assignment of physical variables to the space elements defined in Section III-B does not follow the conventional association rules [7].For example, the stream function is associated with primal nodes instead of the conventional assignment to dual nodes.In this way, the total current in a conductor can be defined as the difference in stream functions values using (9).The coupling (8) and (12) requires to express the current through dual faces ĩ in terms of ψ and the magnetic flux through the dual faces b in terms of a.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply. 1) ψ − ĩ Relation: To map the vector of stream functions ψ defined on the primal nodes to the corresponding vector of stream functions defined on the dual nodes ψ, two steps are necessary.First, the dual grid is augmented by adding the mid-edge nodes to the list of dual nodes [11].In this way, ψ can be easily interpolated using the face-to-node incidence matrix T and the edge-to-node incidence matrix G + s , where all the entries are taken as positive, as shown in Fig. 3(a) The current vector through the dual edges of the surface mesh is then obtained using the discrete gradient on the augmented dual grid Gs,aug where the matrix P projects the currents defined on the 2-D surface mesh into the corresponding currents defined on the 3-D mesh.From a geometrical point of view, P is the matrix that maps the edges of the 2-D mesh into the edges of the 3-D mesh.
2) a − b Relation: The mapping between the magnetic vector potential circulations a to the magnetic flux densities b is obtained in a similar fashion.Initially, the flux densities through the faces of the 3-D mesh that belong also to the 2-D mesh are selected using the selection matrix S. Then the flux through a triangular face is divided in three contributions associated with the dual faces, as shown in Fig. 3(b), using again the face-to-node incidence matrix

D. Final Equations
Combining the equations of the 3-D and 2-D domains ( 8) and ( 12) with the coupling ( 14) and (15), the final system becomes where The source term is set by suitable Dirichlet boundary conditions on a (external field) or ψ (transport current).The nonlinear transient problem is solved using the θ -method (here θ = 0.5) with a fixed time step.At each time step, a fixed-point nonlinear problem is solved.This coupled scheme exploits the initial factorization of the solution matrix that is then used either for the nonlinear iterations or for the time stepping scheme.
The choice of the initial conductivity σ 0 = (1/ρ 0 ) is crucial for the convergence of the nonlinear iterations.A preliminary study, not reported here for brevity, showed that using J = 1.2 J c in (4) provides a good convergence in all the benchmarks.Another key strategy for the convergence is the under-relaxation of the update of the residual term R k .In fact, due to the large variation in the conductivity values of YBCO, the residual term is subject to large variations, especially when the solution is far from the convergence at the initial iterations.For this reason, a dumping factor α = 10 −4 is used to update the residual IV. RESULTS

A. Verification With a Nonlinear Static Benchmark
The first benchmark is used to test the behavior of the proposed formulation with the high nonlinearity of the conductivity (3).It consists of a circular ring made by YBCO tape having the inner and outer radii equal to R i = 50 mm and R o = 54 mm, respectively.A transport current I 0 ranging from 0.1 to 1.1 times the critical current I c is imposed.The mesh consists of ∼207k tetrahedra (surrounding air) with average quality index 0.7430, and ∼24k triangles (superconducting tape) with average quality index 0.9963.The average number of nonlinear iterations to reach a relative error on the solution vector below 10 −6 is 129 for a total simulation time of 28.1 s (factorization: 5.6 s, nonlinear iterations: 22.5 s).Timings refer to a pure MATLAB implementation of the algorithm, running on a workstation equipped with two Intel Xeon Gold 6154 18-core at 3.0 GHz and 512 GB of RAM.
The reference solution of this benchmark can be obtained by subdividing the ring in N t flux tubes (see Fig. 4), where the nonlinear conductance of each tube is The current in each flux tube i k is calculated by solving the system of nonlinear equations Fig. 5 shows comparison of the power losses calculated with the a − ψ formulation with those calculated with (19) using Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.N t = 20.A perfect agreement can be observed, showing the good convergence of the nonlinear scheme also when the transport current is in the order of the critical one.

B. Validation With a Nonlinear Transient Benchmark
The second benchmark is a single-layer, three-tape CORC cable.The three tapes are wounded around a former (not included in the simulation) having 4.76 mm diameter with a pitch of 40 mm [12].The mesh, shown in Fig. 6, consists of ∼67k tetrahedra with average quality index 0.7631, and ∼7k triangles with average quality index 0.9976.The external excitation is a 50 Hz uniform magnetic flux density, orthogonal to the cable axis, with intensity ranging from 10 to 90 mT.The transient simulation is run for three periods, with a fixed step size of 0.1 ms.The numerical results are compared with the experimental measurements reported in [12] and the simulation results obtained in [13] using an H -formulation in terms of ac losses per unit length.Losses are computed in the last of the three periods (60 ms) simulated.The average simulation time for each solution is 1.13 h.The results of Fig. 7 show a good agreement for all the values of applied magnetic field up to 60 mT.Beyond this value, heating effects that lead to the reduction in the value of J c are reported in [12], with a consequent reduction in the measured losses.
V. CONCLUSION In this article the a − ψ formulation in terms of algebraic quantities is presented.The coupling between the 3-D and 2-D equations makes it possible to keep the mesh quality also in case where the aspect ratio is critical.The formulation is verified with a semi-analytical problem and validated versus measurements available in the literature.The results show that the method is effective and efficient in solving problems with superconductive tapes also in case of strong nonlinearities.

Fig. 4 .
Fig. 4. Geometry of the static benchmark with current flux tubes' subdivision (not to scale).

Fig. 5 .
Fig. 5. Power losses for the circular ring static benchmark.

TABLE I THICKNESS
AND CONDUCTIVITY OF THE DIFFERENT LAYERS OF THE YBCO TAPE UNDER STUDY