Numerical Modeling of AC Loss in HTS Coated Conductors and Roebel Cable Using T-A Formulation and Comparison With H Formulation

With recent advances in second-generation high temperature superconductors (2G HTS) and cable technologies, various numerical models based on finite-element method (FEM) have been proposed to help interpret measured AC loss and assist cable design. The T-A formulation, implemented in COMSOL, shows great potential for reducing the overall computation costs. In this paper, the performance of the T-A formulation for calculating the AC loss of coated superconductors and cables were assessed and compared against the widely accepted H formulation, with benchmark model of a single REBCO tape in 2D/3D and a 14-strand Roebel cable. Evaluation and comparison on key metrics including the computation time, the number of degrees of freedom and the numerical accuracy were presented, which could provide a reference for researchers in applying the T-A formulation for AC loss calculation.


I. INTRODUCTION
With recent advances in commercialized second-generation high temperature superconductors (2G HTS), proposals of new cable structures and upgrades of previous ones are currently under way, including Roebel assembled coated conductor [1], conductor on round core (CORC R ) [2] and twisted-stacked-tape cable (TSTC) [3]. For instance, Roebel cables, with its compact design and high current-carrying ability, are used in the insert HTS magnet in FRESCA2 [4], [5]. A solenoid demo magnet with CORC R cables retained the critical current over 4 kA in a background field of 14 T, generating a total center magnetic field of 15.86 T [6]. Prototypes of other candidates, such as CroCo [7], VIPER [8] and STAR wires [9] are yielding promising results.
The AC loss in HTS cables increases the refrigeration load to the cryogenic system and affects the efficiency of practical devices [10], [11]. To predict the loss behavior and optimize cable designs, efficient and accurate modeling tools are crucial. The main challenges when modeling HTS machines The associate editor coordinating the review of this manuscript and approving it for publication was Guido Lombardi . and cables are related to the electromagnetic properties [12] and geometrical complexities [10]. The sharp transition from the superconducting state to the normal state and critical currents depending on both the magnitude and direction of the magnetic field are embedded in the constitutive laws [12], [13]. On the other hand, the large aspect ratio of the superconducting layer and 3D structure of cables dramatically increase the size of the problem [14].
Various modeling strategies have been proposed over the years, both analytical [15] and numerical [16], [17]. Amongst them, the finite element method (FEM) with different formulations has been widely accepted and intensively used for its advantage in tackling complex geometries, compatibility with multidisciplinary simulations and transferability in the scientific and industrial communities [10]. In particular for HTS cables, FEM models in 2D with simplified geometries or modified equations, and 3D with full structures were developed, including Roebel [18]- [21], CORC R [22]- [24] and TSTC [25], [26].
In this paper, we focus on the performance of the T-A formulation for calculating the AC losses of HTS coated conductors and cables, with comparison against the H formulation, VOLUME 9, 2021 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ both implemented with the commercial software COMSOL Multiphysics R in 2D and 3D. Main objectives of this study are to test their performances on key metrics including the computation time, number of degrees of freedom (DOFs) and accuracy using benchmark models, and to provide a reference for researchers and engineers in the application of T-A formulation for AC loss calculation. The models, test cases and modeling details are described in section II-III. Results and discussions on single tape benchmark in 2D/3D and a 14-strand Roebel cable are presented in section IV-V. Finally, the main conclusions are summarized in section VI.

II. MODEL DESCRIPTION A. T-A FORMULATION
The T-A formulation was proposed in [27] and [28], and then extended with homogenous and multi-scale concepts in [29]. The current vector potential T T T is defined in the superconducting regions and used as a state variable to calculate the current distribution in superconducting regions, with Faraday's law where ρ is the resistivity of the superconductor. The magnetic vector potential A A A is defined in both superconducting and non-superconducting domains and solved using the A formulation Relative permeability µ r = 1 is assumed in this paper. The surface current is calculated from T T T and imposed to the superconducting boundaries in A A A. The normal component of B is obtained from A A A and serves as a source term for T T T [27]. The thin-strip approximation was applied in the initial version of T-A and the case with finite thickness was presented in [30]. The constrains on the transport current are implemented by imposing the appropriate sets of boundary conditions with variable T T T on the ends or edges of the superconductors. With great potentials in multi-physics coupling and 3D modeling, models of HTS dynamo [31], quench [32], and cables with 3D structure [23] were developed. In our study, we used the original T-A formulation in 2D and 3D, as in [27] and [28], which considers a uniform distribution of current density across the thickness. It should be noted that in cases like the closely packed bifilar stacks [33] and coils [34], where the penetration of the top/bottom faces should be taken into account, a full discretization across the thickness as in the H formulation described below would be preferred.

B. H FORMULATION
The H formulation has been widely used and adopted as a standard for new models [10], [35] since it was proposed in [36], [37] in 2D, and [38], [39] in 3D, with the governing equation where H H H is the state variable. The transport current is applied by integrating the current density over the cross-section of the superconductor. Iterative multi-scale and homogenized models were developed for large-scale applications [40], [41]. The Magnetic Field Formulation (mfh) module based on H formulation was integrated in COMSOL Multiphysics R since version 4.3b: from an application perspective, this allowed a more convenient implementation of the H formulation than the original ones based on the Partial Differential Equation (pde) module [42]. In terms of AC loss of HTS cables, the 3D H formulation has been applied to many structures including CORC R [22], Roebel [19] and twisted multifilamentary MgB 2 wires [43]. It is worth noting that the coupled A-H formulation [44], proposed for the modeling of superconducting electrical machines, is also among similar methods. Although not included in this paper, such discussions could be found in [31] and [45].
For the sake of convenience, the different models used in this study depending on their element types are designated as: • T-A(1+1): T-A formulation with linear Lagrange elements for T and linear elements for A; • T-A(1+2): T-A formulation with linear Lagrange elements for T and quadratic elements for A; • T-A(2+3): T-A formulation with quadratic Lagrange elements for T and cubic elements for A; • H(1): H formulation with first-order curl elements; • H(2): H formulation with second-order curl elements. The validation of our models based on the T-A and H formulation is presented in the appendix, where simulations of a single pancake coil are compared against the measured data.

C. TEST CASES
Two sets of geometry were chosen as test cases. One is a standard 4-mm wide coated conductor, and the other is a 14-strand Roebel cable. Details about the test cases are as follows.

1) SINGLE TAPE
The superconducting layer of commercial REBCO coated superconductors is typically 1-2 µm [46]. In 2D models, the superconducting tape is assumed to be infinitely long, which is ensured by the formulations. With the T-A, it is valid to apply the thin-strip approximation to our test case of a single tape, and convenient to impose the transport current constraint by a set of Dirichlet boundary conditions, as proposed in [27] (Fig. 1(a)). While with the H formulation, considering the balance between computation efficiency and accuracy, the thickness of the geometry is artificially scaled up to 10 µm [47] (Fig. 1(b)). Only the superconducting layer is modeled. Meshes used in 2D are also illustrated in Fig. 1, where dmesh is defined as the number of elements on the superconductor in 2D. Although the same control parameters are used to discretize the air domain, the generated meshes are quite different. Using thin-strip approximation, rectangular  elements with high aspect ratio could be avoided, thus the meshes are better structured.
For 3D models, a total length of 2 cm is analyzed. The ideal magnetic insulation boundary conditions are assumed at both ends. The geometries are first divided into smaller parts, which are rectangles in 2D and cuboids in 3D. Here, dmesh is defined as the number of rectangles/cuboids along the cross section, as noted in Fig. 2. They are then converted to triangular/tetrahedral elements by inserting diagonal edges, which is available in the Mesh module in COMSOL Multiphysics R . Meshes used in the 3D single tape models are shown in Fig. 2

.(a) with T-A and (b) with H.
The E − J power-law model was used in all cases, with I c = 140 A, E 0 = 1 µV/cm and n = 21.

2) ROEBEL CABLE
A 14-strand Roebel cable fabricated by the research group at Karlsruhe Institute of Technology (KIT) is chosen. Previously, comprehensive studies with AC loss measurements [48] and numerical analyses [49], [50] were conducted. A full 3D model with the H formulation and periodic boundary conditions can be found in [19]. Key parameters used in this paper are the same as ones in [19] and summarized in Table 1. For simplicity and consistency with previous publications, we used a constant J c with total critical current I c = 465 A, as in [19] and [49]. Fig. 3(a) and (b) shows the geometries and meshes (dmesh = 20) used in this study.
All the simulations in this paper were carried out on a PC with Intel (R) Core (TM) i7-7700K CPU @ 4.20 GHz processor and 32 GB of memory, and solved using COMSOL Multiphysics R 5.4 [47], [51]. The H formulation was implemented with the embedded General Partial Differential Equation (PDE) module and the T-A with the Magnetic Fields (mf) module coupled with PDE module in lower dimensions. For all studies, 101 steps over one period (0.02 s) were saved for the post-process. Detailed configurations and justification for the parameters in each model are explained in the next section.

III. MODEL DETAILS A. ELEMENT ORDER
As presented in several reports [37], [39], the common practice when modeling with H formulation is to use edge elements (or 'curl elements' in COMSOL R ). In our study, both first and second order curl elements were considered in 2D, and first-order curl elements were tested in the case of 3D.
With the T-A formulation, it was reported in [29] and [31] that with first or second order Lagrange elements for both T VOLUME 9, 2021  [19]. In consistence with strip approximation, one layer of elements is considered across the thickness of the strand. Meshes with dmesh = 20 are plotted, respectively. The size of the air domain shown in this figure is not to scale. and A, spurious oscillations could appear in the current density J profile at certain study steps. However, because such oscillations mostly occur within the subcritical region, where the current density is lower than the critical current density J c , the influence on the resulting AC losses is insignificant. This pattern of oscillation was found in our study as well. For instance, in Fig. 4(a) with transport current with amplitude of I 0 =0.9 I c and 100 elements across the tape width at 0.02 s with first-order Lagrange elements, each neighboring point alternates above and below the expected value. Plot after 'smoothing' [52], which is a post-processing option provided by COMSOL R to average any two values with the same coordinates, is also shown.
In fact, as stated in the governing equations of T-A formulation, the current density J takes the first derivative of the state variable T and the second derivative of state variable A. To match the levels of numerical accuracy in coupled equations, the element order used for discretizing the A-field should always be one order higher than that in T . As demonstrated in Fig. 4(b), where paired element orders (linear elements for T , quadratic for A and quadratic elements for T , cubic elements for A) were applied, there was no oscillation observed. Also, the smoothed curve in Fig. 4(a) followed the trend quite closely. In each test case, choices of element orders were made accordingly.

B. TOLERANCE
A tolerance criterion is set up in the adaptive timestepping solver (DASPK), which is the default solver for time-dependent studies in COMSOL Multiphysics R , to determine the accuracy of a study step and terminate the iteration, defined as [52] where M is the number of fields, N j is the number of DOFs in field j, A i is the absolute tolerance (atol) for DOF i, E i is the estimated absolute error, U i is the solution, and R is the relative tolerance (rtol). Despite the various methodologies to adjust the tolerance factors, for convenience, atol is set as a scaled factor of 0.01 and 0.1 for 2D and 3D, respectively, throughout the following studies. Fig. 5(a) shows the calculated AC losses with 0.9 I c in a single tape, with regard to tolerance rtol and meshing quality. In each case, a certain threshold of tolerance factor is required to converge to a satisfactory solution. Generally, a larger tolerance may lead to errors in the current distribution and subsequently errors in AC loss, and a smaller tolerance could guarantee a better convergence but would take longer computation time. With finer meshes, a stricter tolerance factor is needed. In theory, tolerance factors should be checked for each model for optimal performance. However, significant errors in calculated AC loss almost always correspond to clear spikes in instantaneous losses. Fig. 5(b) shows the instantaneous losses with rtol = 5e-4, 1e-4 and 1e-5, and 100 meshing elements, where irregular spikes can be easily recognized. Monitoring the instantaneous losses in each time step is a practical way to check the calculation, as recommended in [51]. Spikes in instantaneous losses may suggest that the tolerance configuration should be tightened. An alternative is to change the update frequency of the Jacobian matrix in Newton iteration to 'once per time step' or 'once every iteration'. This could improve the stability of the solver and avoid applying extremely tight tolerance factors throughout the calculation. Still, longer computation time should be expected.
In practice, the determination of tolerance settings may take some trials and errors, depending on each individual model and accuracy requirement. In the following sections, the relative tolerance in 2D and 3D models are 1e-5 and 1e-3, respectively, and Minimal update frequency is used, unless otherwise stated.

IV. TEST CASE 1: SINGLE TAPE
Several criteria are established to assess the performance of T-A and H formulation, including the total computation time, number of DOFs and solver time per thousand DOFs. For accuracy, solutions with the classic Norris's analytical model [53] are adopted as a reference, expressed as where I 0 is the peak transport current and I c is the critical current. Nonetheless, this is a solution to superconductors with the critical-state model (CSM) instead of power law as used in FEM. Here we used the relative difference in AC losses when doubling the number of elements across the tape width to represent its sensitivity to meshing quality, which is defined as where m is a positive integer, Q m is the AC loss in J/m/cycle, dmesh is the number of elements on the superconductor in 2D (Fig. 1), or correspondingly a control parameter in 3D as explained in Fig. 2. The relative error r defined above, together with the absolute value of AC losses, should give a comprehensive picture of the accuracy of each model.

A. COMPARISON OF 2D MODELS
In 2D, all cases except for T-A(1+1) were tested. Fig. 6 shows the total number of DOFs with dmesh. DOFs of T-A(1+2) are comparable to or larger than that of H(1), in spite of the simplified geometries. This is because in T-A(1+2), second order Lagrange elements are assigned to elements in the air domain, which needs to be sufficiently large in FEM. However, in H(1), the variable H possesses two components, which are discretized with first-order curl elements. With higher-order discretization, less DOFs are generated with T-A(2+3) comparing to H(2). Calculated AC losses with low (0.2 I c ), moderate (0.6 I c ) and high (0.9 I c ) transport currents are summarized in Fig. 7, where the dash-dotted lines are derived from Norris's analytical model. Overall, the computed AC losses converge with increasing number of elements in all four cases. With lower transport currents, a finer mesh is required. This is consistent with the pattern of current distribution in a superconducting strip with transport current, which indicates possible improvements of the FEM models by adjusting the distribution of meshing nodes without increasing the total number of DOFs. In fact, such principles are common practices in FEM modeling. For instance, in both homogeneous T-A VOLUME 9, 2021 formulation [29] and H formulation [40], carefully arranged meshes are recommended for better performance.
With respect to the absolute values, four FEM models with the power-law model converge to similar values when discretized with sufficient elements. The discrepancy between numerical models and Norris's analytical model could largely be attributed to the different constitutive laws, which is the E-J power-law model with n-index of 21 in FEM comparing to CSM. This could be further mitigated by increasing the n-index to imitate the immediate transition into the normal state as in the CSM [54], [55].
When increasing dmesh, the trends of convergence with corresponding element orders are close, for instance, H(1) with T-A(1+2), and H(2) with T-A(2+3). The difference between T-A and H formulation with dense meshes is minimal and mainly results from the assumptions to the thickness of the superconducting strip. With thin-strip approximation, the assumed thickness in T-A is effectively a scaling constant with no real impact on the output. But with H formulation, AC loss is dependent on the thickness of the geometry. As suggested in [47], with middle to high transport currents, a thickness of no greater than 20-30 µm is preferred. 10 µm is used in our study, and the error related to tape thickness should be minimal in most cases. Fig. 8 shows r , as defined in (7), and solver time with transport currents 0.2 I c , 0.6 I c and 0.9 I c . The jumps of computation time with 0.9 I c , dmesh ≥ 200 and 0.6 I c , dmesh ≥ 320 are due to change in the update frequency of Jacobian. The basic trade-off remains that assigning more elements at the superconducting strip leads to better accuracy as well as longer computation time. With higher transport current and higher element order, fewer uniformly distributed meshing elements are required to converge to a stable value.
To take full advantage of the current-carrying ability of HTS tapes, higher operating current is preferred in most applications [47]. In practice, it is possible to lower the number of degrees of freedom while achieving a reasonable level of convergence. Specifically, with T-A in Fig. 8(b), if a relative error threshold of 0.5% with 0.6 I c was set, a minimum number of 40-50 elements would suffice for T-A(1+2), and 10-16 for T-A(2+3).
In terms of computation time, with dmesh ≤ 100 in H(1) and dmesh ≤ 50 in H(2), the computation time of H(1) and H(2) is comparable to that of T-A(1+2) and T-A(2+3). But as dmesh increases, it rises sharply with H(1) and H(2) while remains relatively stable with T-A. Considering that the number of DOFs of H(1) is similar to or even less than that of T-A(1+2) (Fig. 6), the improvement of efficiency with larger dmesh is mostly attributed to different formulations. With H(2) and T-A(2+3), it should be both the reduction of DOFs and difference in formulations. This also indicates that T-A formulation could show greater advantage in saving total solver time with large-scale models in 2D.
Considering the size of the model, the number of DOFs with H(1), T-A(1+1) and T-A(1+2) are plotted in Fig. 9.

Number of DOFs of T-A(1+1)
is about 30% less than that of H(1), and the difference diminishes with increasing dmesh. However, the number of DOFs with T-A(1+2) is 3 to 4 times that of H(1). This is mainly because in 3D, there are three components in variable A. Comparing first-order curl elements with second-order Lagrange elements, the latter could be more costly.
AC losses with respect to dmesh in 3D in Fig. 10 share similar characteristics as in 2D in Fig. 7. The H formulation has greater errors when the meshing quality is low, but as dmesh increases, the accuracy of all three models are close, with the exception of I 0 = 0.2 I c . Since it is neither commonly required nor practical to achieve high computation accuracy in 3D with low applied current, the case with 0.2 I c is ignored in the following discussions. Fig. 11 shows the normalized current density J x /J c profile with dmesh = 32 and I 0 /I c = 0.9 at t = 0.02 s in T-A(1+1), T-A(1+2) and H(1). In Fig. 11(a) with T-A(1+1), 3D oscillation pattern appears in neighboring domain elements. Similar to that in Fig. 4 in section III-A, the normalized current density profiles in the critical regions remain stable while oscillations mostly occur in the subcritical region. As for the AC losses, the difference between the element orders of variable A does lead to larger errors in T-A(1+1) comparing to T-A(1+2) with the same geometric meshing.
As far as efficiency is concerned, the total computation time and relative error r are plotted in Fig. 12. It is shown that T-A(1+1), with the least DOFs and shortest solver time amongst all three, requires better resolution of meshes to compensate for the lower-order shape functions comparing to T-A(1+2). Contrary to 2D, H(1) in 3D is less time-consuming than T-A(1+2) with the same dmesh. T-A(1+2) is both of the highest accuracy and requires the longest computation time with this model. However, as demonstrated in Fig. 6 in 2D and Fig. 9 in 3D, models with second-order elements of variable A generates much more DOFs in 3D but not in 2D. To take this factor into account, the average time per thousand DOFs is compared in Fig. 13. In a double log plot, the slope of H(1) is about VOLUME 9, 2021 twice that of T-A(1+1) and 1.5 times that of T-A(1+2). This indicates that H(1) has advantages in total solver time and DOFs, as well as comparable time per DOF with models of relatively small scale. But with larger superconducting domains, the average processing time in H(1) is considerably longer as condition of the matrix worsens. In this case, the total computation time would depend on the tradeoff between the increase in the number of DOFs and the average time per DOF.

V. TEST CASE 2: ROEBEL CABLE
In this section, we take the example of AC losses with transport currents of the Roebel cable described in section II-C-2. Performance of T-A(1+1), T-A(1+2) and H(1) with meshes in Fig. 3 are summarized in Table 2.
In consistence with conclusions from the previous section, H(1) has about 1/4 the DOFs of T-A(1+2) but longer average time per thousand DOFs, resulting in almost the same total computation time as T-A(1+2). T-A(1+1), although with 'oscillation' patterns, requires the least DOFs as well as the least time. The AC losses calculated with three models are reasonably close. Using H(1) as a reference, the relative difference is within 5%. Overall, in some cases, despite the large matrix size, T-A(1+2) could be a good choice for full 3D modeling of HTS cables with a balance of satisfactory accuracy and acceptable computation time. However, more detailed studies are required to predict whether T-A(1+2) or H(1) would have the advantage in total solver time, especially with a new model or a different set of computer configurations. On the other hand, the application of T-A(1+2) may be limited by the computer memory, due to the drastic increase  in the number of DOFs. Then T-A(1+1) with a finer mesh could be a good alternative with reasonable numerical errors. Fig. 14 shows the normalized current density on the crossover region with peak current. In Fig. 14 (c) with H(1), the outer surface of the crossover region carries larger currents than the inner side, which is similar to the observation in [19]. This difference across the thickness of the strand is suppressed with thin-strip approximation, where a uniform distribution of current density is assumed, as in Fig. 14(a) and (b). Although the overall current distribution and the total AC loss are similar in this case, this effect could be more significant in others, for instance, when Roebel cable is subjected to an external field parallel to the strand surface. Comparing to thin-strip approximation, H formulation in full 3D retains the ability to characterize the field penetration from the wide face of the strand, and as a result, may affect the overall AC losses.
The full AC loss versus normalized applied current is plotted in Fig. 15, with solutions from T-A(1+2), dmesh = 64 in 2D and T-A(1+1)/T-A(1+2), dmesh = 20 in 3D. The experimental data and simulation results with critical-state model were obtained from [48] and [19]. At low current level, losses with T-A(1+2) are consistent with those predicted by the H formulation in both 2D and 3D, and higher than the measured data, while losses from 3D T-A(1+1) are even higher. It is highly probable that calculation errors relating to poor meshing quality could be non-negligible at this transport current level. However, it is hard to confirm without further refinement. With high transport currents, all simulation models yield similar results. The measured AC losses are higher than the calculated ones, which may result from the coupling losses amongst strands in experiments or the anisotropy of critical currents. The computation time of 3D T-A(1+2) presented in Fig. 15 varies from 4.8 to 12.7 hours (0.22 I c to 0.73 I c ), with spikes over 0.77 I c , comparing to 1.7 -3.7 hours (0.26 I c -0.86 I c ) with T-A(1+1).

VI. CONCLUSION
In this work, we presented an assessment of the T-A formulation for calculating the AC loss of HTS coated conductors, compared against the H formulation, concerning number of DOFs, computation time and accuracy. With sufficiently fine meshes, all models showed good agreements with analytical solutions as well as with each other.
In the 2D single tape model with coarse mesh, the total solver time and number of DOFs of T-A and H with equivalent element orders are comparable. With further refinement of the mesh, T-A(1+2) requires much less computation time than H(1) even with a greater number of DOFs. The accuracy was directly related to the element orders. In the 3D model of a single tape, quadratic elements of variable A in 3D T-A(1+2) lead to a drastic increase of DOFs and, subsequently, total solver time. Whereas T-A(1+1) exhibits good accuracy of AC loss with refined meshes. However, compared to H(1), the average time per thousand DOFs of T-A(1+1) and T-A(1+2) stays low.
As the complexity of the model increases, like the full 3D 14-strand Roebel cable, T-A(1+2) remains advantageous in terms of time per thousand DOFs but requires much greater computer memory, when compared with H(1). Alternatively, T-A(1+1) could reduce the number of DOFs and the total solver time, but extra attention should be paid to the oscillation patterns and meshing quality.
It should be pointed out that, throughout this study, only the T-A formulation with thin-strip approximation was considered.

APPENDIX
To validate the numerical models in this study, the measured transport AC losses of a 24-turn pancake coil in liquid nitrogen presented in [56] were compared against simulations using H(1) and T-A(1+2). Details about the experiments were presented in [57]. The anisotropic field dependence of J c was modeled in accordance with the method in [57] and [56]. Power-law with E 0 = 1 µV/cm and n = 18.3 was used, as in [56]. Fig. 16 shows a good agreement both between the simulated results and with the measured data. The maximum relative differences with respect to the measurements below the critical current are about 20%.