H-Formulation FEM Modeling of the Current Distribution in 2G HTS Tapes and Its Experimental Validation Using Hall Probe Mapping

One of the most widespread mathematical formulations applied to simulate the electromagnetic phenomena of coated conductor in the recent literature is the H-formulation. However, the only validation of the model has been indirect by using measurements taken from the applications, as measurements of the energy losses in ac fields, forces developed in levitation systems, or any other parameter related to a specific application. Direct validation of the calculation requires the observation of the local out-of-plane magnetic field over the surface of the sample, and it is only accessible under magneto-optical observations and, in a larger scale and better dynamic range, by the Hall scanning microscopy. We propose here the experimental validation of the H-formulation by comparing the simulated results with measurements made by a Hall probe mapping in a second-generation (2G) tape sample for several dc transported currents at 77 K. This paper presents a methodology to simulate the 2G tape by using only measured data obtained from a sample and its normalized J(B) experimental curves. Some boundary conditions that allow a faster convergence of the problem are investigated. Simulated results of the 2G tape modeled considering only the 1-μm high-temperature superconductor (HTS) layer were compared with others that represent the most important layers of the coated conductor structure in the calculations. The simulated and measured results present a good agreement, proving that this model can calculate precisely the magnetic field and, hence, the current distribution in HTS samples.


Introduction
The recent coated conductors quality improvements allowed to reach a new level of technology development in superconducting large scale applications [1][2]. In order to design devices using superconducting tapes, it is necessary to use computational tools to simulate those materials considering the nonlinearity intrinsically associated to them.
The modeling of high temperature superconductors (HTS) was vastly investigated in the literature, and there are several formulations applied to solve the electromagnetic problem. The formulations more frequently found in the literature applied to solve the problem by the finite element method (FEM) are: the H [3][4][5], the A-V (or A-φ) [6][7][8], the T-Ω [9][10][11] and the E [12] formulations. All those formulations were used to calculate the current distribution inside the superconductor for bulks and tapes. For the case of the coated conductors, one of the models applied frequently in the recent literature is based on the magnetic field strength, also named as H-formulation. One reason for that is the simplicity to implement it in FEM commercial programs, like COMSOL.
Considering the H-formulation, the previous works dedicated to the validation of the model only contrasted the calculated results with measurements provided by a mean value of the current distribution, which represents an indirect comparison. Generally, the results are presented with good agreement [4,[13][14][15][16] for the energy dissipated (power losses) in a superconducting tape transporting an AC current. Also, the Hformulation by FEM simulations was successfully compared with levitation force measurements [17]. Besides, the simulated results for those models were compared with analytical equations [4,18].
The current distribution inside a superconducting tape is directly associated with the external magnetic field produced by it, which can be obtained by mapping its surface using a Hall probe. Using an inversion method based on the Biot-Savart's law, the current distribution in a HTS tape was obtained [19][20].
In this context, the present work has the following contributions: i) To validate the calculated current distribution inside a 2G tape obtained by the Hformulation comparing it with Hall probe mapping measurements; ii) To simulate the 2G tape only using experimental parameters easily obtained from a sample and the manufacturer data; iii) To investigate the influence of relevant parameters in the simulated results, such as consideration of the layers in the modelling, and comparison between different boundary conditions.

The H-formulation
The H-formulation has the advantage of solving directly the differential equations for the magnetic field. It has been extensively presented in the literature to calculate the current distribution in superconductors [3-5, 13-16,18]. By the Faraday's law, for a 2D symmetry presented in Fig. 1, where L is considered infinity long and only the cross section is simulated, it is possible to write the following two partial differential equations (PDEs): By the Ampere's law (neglecting the displacement current) for the symmetry presented in Fig. 1, and the power law, it is possible to obtain a constitutive equation to represent the flux creep and flux flow in the superconducting region, as follow: where E c will be adopted here as 1 µV/cm, and the critical current density J c and n are functions that can be obtained by experimental results or by analytical expression with adjusted parameters.
To solve the problem, equations (1a), (1b) and (2) can be introduced in a software with capacity to resolve generic PDE. In the case of the present work, the problem was solved by the PDE module in COMSOL 4.4, in terms of the magnetic field components (H x and H y ). For the air region, the first term in the right side product of equation (2) was replaced by 1 Ω•m constant, and it must be nonzero in order to avoid a singularity in the system of equations [12]. It is an issue introduced to force a low resistivity air value and obtain the problem convergence.

Boundary conditions and multilayer modeling
The type of element used, nodal (also named Lagrangian) or edge (also known as rotational), will influence in how the boundary condition should be imposed to the problem and also in how fast the solution is calculated [21]. The application of the edge elements in the H-formulation intrinsically guarantees that the ∇⋅B = 0 condition is satisfied [4], since this condition is present in the solution of each element. Furthermore, it was presented previously that both element types can present the similar results [22].
Several arguments were presented in [4][5] showing the advantages of the edge elements, since it presents a better distribution of the matrix and a faster convergence.
In this work the edge elements will be used. Three cases will be simulated and The boundary condition of cases (i) and (ii) was originally proposed in [18]. It reduces the number of DOF by working with a narrow external air domain, in other words, by having the borders very close to the superconductor. The magnetic field is separated in two components, to be applied in the external border lines: ( , , ) = ( , , ) + ( , , ).
The variables and correspond to the external magnetic field in the x and y directions respectively. In the present work, no external field will be applied, and so = = 0. The variables and correspond to the magnetic field produced by the currents flowing through the superconductor.
For cases (i) and (ii), an intrinsic property of the edge elements is explored, since they have a null component along the perpendicular direction to their edges [4][5]. It simplifies the implementation of the boundary conditions, since equations (3a) and (3b) are respectively applied in x and y border lines (see Fig. 2). By the Biot-Savart's law, the magnetic field produced by the 2G tape currents can be obtained as follow: The simulations presented in this paper are centered in a geometry that corresponds to a 2G tape of 12 mm wide and 110 µm thickness, with the x and y external borders of 6 mm and 1 mm respectively, as shown in Fig. 2. In order to calculate the current density distribution inside the tape, the value of the net current must be inserted in the formulation. In the case of the COMSOL software, it allows this constraint to be enforced by means of Lagrange multipliers, as described in [4], forcing the result of the net current to be equal to the current density integral over the HTS cross section area.

Critical Current Dependence with the Field
There are several ways to model the critical current distribution dependency with the magnetic field. The more general way is the consideration of an analytic expression for J(B) with adjusted parameters [23], in order to fit the simulated and measured results.
Instead of choosing some of the models presented in Table I, we are proposing to use directly the normalized J c (B) curve expression for a 2G tape. This curve can be measured for a given sample.
The I c angular dependence is typically described by the manufactures, by giving two sets of curves, corresponding to B//c and B//ab. In the actual state of the art, the successful introduction of pining centers induces a different angular behavior, but the pattern is well known as its temperature dependence. Both, the temperature dependence with B//c and B//ab and the angular pattern are currently available at the main produces.
Measurements of J c (θ,B,T), are somewhat difficult. They are made typically in terms of I c . One of the points of difficulty is originated by the self-screening effects of the tape, which induce a hysteretic behavior. In order to measure the intrinsic value, characterization is made in small micrometric wide bridges to approach it, neglecting so the screening effects and currents distribution that occurs in the wider samples. Main manufactures provide this kind of measurements at different temperatures and magnetic induction.
For 2G advanced pinning tapes operating at 77 K with magnetic flux density under 100 mT, the observed anisotropy effect is not so severe, allowing an approximating of the problem without considering the angular dependence of the current density. In our case, the I c (B) data were obtained from the manufacturer [33] for an applied DC magnetic field. The critical current was divided by the coated conductor area and normalized for the zero field value. The data is introduced in the program in a lookup where J c0 is the critical current density at zero field and |H µ 0 | is the local magnetic field magnitude. The variable J c is applied in each element to obtain the solution of the problem.
There are two reasons for modelling the problem in this way: it is not necessary to adjust any parameters present in the analytical equation and it was observed a faster convergence than using one of the expressions presented in Table I.
The data were introduced as a table for linear interpolation as shown in Fig. 3. Only the J c (B) component parallel to the c axis was implemented in our model. Obviously the modeling of angular dependence could better represent the physics of the problem, but it was neglected here for simplicity. Fig. 3: Normalized critical current density data obtained from measured results [33] for a 2G tape similar to the studied in this work.

Current Density Distribution in a 2G Tape
All the parameters used in the simulations were based in a commercial 2G tape. The sample tested was produced in 2015 by SuperOx and its dimensions are: 12 mm width, 110 µm thickness and 170 mm long. More information about it and other parameters used for the simulation are presented in Table II.  In all simulated cases, only first order edge elements were used. For cases (i) and (ii), the mesh was created only using quadrilateral mapped elements. For case (iii), the mesh in the areas presented in Fig. 2 is the same of the one made in case (ii), and additionally a free triangular mesh was created externally to this region.
In the simulations, the current was applied in the tape by a sequence of step functions, with an interval of 1 s to reach the next value and keeping the signal constant during 99 s, as presented in Fig. 4. The current was maintained constant during this time in order to allow the flux creep effects pass and to stabilize its distribution along the superconductor. and L2 are compared for the cases (i), (ii) and (iii). Fig. 5 presents the normalized current density calculated for the coarse mesh along L2 (Fig. 5 a) and L1 (Fig. 5 b), for the applied current profile of Fig. 4 during the rise. The results are presented for a time of 99 s after the previous step, when the current distribution is stabilized inside the superconductor. Fig. 6 presents the same results for the current decreased.   In the other hand, the fine mesh presents a profile shape in agreement to the experimental one, but the computation time increases as presented in Table III. From now on, all the simulated results will be presented for the case (ii) and the fine mesh. The formulation used in this work can take into account the flux creep effects observed in experiments with the superconducting material. The capacity of calculating precisely this phenomenon is a sine qua non condition for any model to be implemented in several 2G tape applications. After a step in the applied current (Fig. 4), the flux relaxation will make the current density reduce in magnitude and penetrates more in the tape. Since the flux creep has a logarithmic behavior, after some time the current density variation can be neglected. Fig. 9 presents the normalized current density for several instants along the cut lines L1 and L2. The arrows show the evolution of the current, where each line corresponds to 1 s frame, during 99 s. The same result is presented in Fig. 10 for a 300 A current. It can be noted that the flux creep effect is stronger for 300 A than for 200 A case.   To investigate the effect in other currents, Fig. 12 presents the value of the magnetic induction in the position of the maximal value of the profile (x = -6 mm) for several different currents as a function of the time. Fig. 12 shows that the decay of the magnetic induction for 300 A is higher than for the others currents, indicating that the flux creep is more intense for currents in this value. This figure also presents that after 99 s, the value of the magnetic induction is almost stabilized, and this time is enough to compare the simulated results with the measurements, which will be presented in the next section.

Model Validation
This section presents the Hall scanning magnetometer [36][37] results on a 2G tape sample and compare them with the FEM simulations. A Hall probe sensor placed in a flexible strip was adjusted to gently touch the sample. In this condition, the active part in the sensor (active area of 300 x 300 µm 2 ) has a distance of 400 µm from the sample surface. And so, the distance from the HTS layer and the active part of the Hall sensor is approximately 450 µm. For this reason, the simulated results will be presented for y = 450 µm (Fig. 1). To make the mapping of the magnetic induction B y (x,z), the Hall sensor was placed in a xyz system. All the measurements were made with the tape and the sensor in a LN2 bath (77K). The Hall mapping system precision is 0.01% for full scale range of 0.5 T. Fig. 13: Experimental rig designed to measure the external field produced by a transport current in a 2G tape sample.
The magnetic induction scanning was made for the current profiles presented in Fig. 4.
After each current was applied, an interval of time higher than 120 s was waited before the map was started. However, it is not a problem to compare the simulated and measured results, since the magnetic induction practically does not change after a longer time, as presented in Fig. 12, due to the logarithmic behavior of the flux creep.
To determine the critical current a nanovoltmeter was used to measure the voltage between two points separated 10 cm in the sample. The V(I) curve was measured in the sample as presented in Fig. 14, in order to obtain the critical current (I c ) and the n exponent of equation (2). The determined I c is 447 A (1 µV/cm) and the n exponent is 32, after fitting the measured data. We can observe that the results of Fig. 15, for increasing currents up to 300 A, present a better agreement than the results of Figs. 16 (currents closer to I c ) and 17 (decreasing currents). For the current rise ( Fig. 15 and Fig. 16 a and   The sample inhomogeneity can be estimated by a magnetic induction scanning map of an extension area. Figs. 18 and 19 show a B y (x,z) map for 400 A (during the current rise) and 0 A (remnant field), respectively. In both figures it is possible to observe that the magnetic field profile is not exactly the same for two different z positions. The maximum fluctuation observed in the measured peak value of the 400 A map (Fig. 18) was about 6 %. For the remnant state (Fig. 19) this oscillation along the measured peak values was about 16 %. It means that if an adjustment in the simulation parameters were made to reproduce a determined experimental result for one cross section along z, the same adjust would not fit exactly to another line.

Conclusion
This work presented a validation of the H-formulation applied to calculate the current distribution in a 2G tape. The simulations were made only using parameters obtained from the experimental results. Instead of using any analytical expression with adjusted parameters for the critical current dependency with the field, the normalized J c (B) data, obtained from the I c (B) curves supplied by the manufacturer, was applied directly in the computation code. The J c angular dependency in the 2G tape was neglected here.
The critical current I c and the n exponent were obtained from the I(V) curve measured in a sample, and these parameters were introduced into the simulation code. By the applied methodology it is not necessary to determine any constant to adjust the simulated data in order to fit with the experimental results with the calculations.
Two boundary conditions were compared: a boundary far from the HTS with zero field condition and a boundary close to the HTS with an imposed field in the border determined by the Biot-Savart's law (integral equation-FEM coupling). The calculated current distribution inside the 2G tape was the same for both cases, but it was possible to reduce drastically the computational time in the last one.
Simulations in which the 2G tape was only modelled by the 1 µm HTS layer were compared with others where the most important layers were included, such as the substrate, the silver and the copper covering. The current distribution inside the HTS layer was the same for all the studied cases, which considered a value up to 1.03 I c . The results showed that if the transported current is not much higher than I c , only the superconducting layer should be modelled, allowing a faster computational solution.
The model was able to calculate the external field produced by the 2G tape in agreement with the experimental result, presenting a maximum discrepancy about 10 % and 30 %, for the currents rise and drop, respectively. Furthermore, fluctuations were observed along the magnetic field mapping, and they can be attributed to inhomogeneities in the sample. The maximum fluctuation observed in the B y (x,z) peak value of the 400 A and the remnant state were approximately 6% and 16 %, respectively. The fluctuation can partially explain the observed discrepancy between measurements and simulations, since if the simulation parameters were adjusted to reproduce a determined experimental profile B y (x) of the map, it would not fit exactly to another line.