Advanced electromagnetic modeling of large-scale high-temperature superconductor systems based on H and T-A formulations

The H formulation uses the magnetic field strength ${\mathbf{H}}$ as dependent variable. Within a bounded universe the different materials are represented by different subdomains. Each subdomain has different properties, i.e. resistivity $\rho$ and permeability $\mu$. These values define the constitutive relations ${\mathbf{E}} = \rho {\mathbf{J}}$ and ${\mathbf{B}} = \mu {\mathbf{H}}$, where ${\mathbf{E}}$ and ${\mathbf{H}}$ are the electric and magnetic field strength, respectively, and ${\mathbf{B}}{\kern1.5pt}$ and ${\mathbf{J}}$ are the magnetic flux and current density, respectively.

To derive the governing equation of the H formulation, Ampère's law is written neglecting the displacement current. Then, Ampère's and Faraday's laws are given by

$$\begin{equation}\nabla \times {\mathbf{H}} = {\mathbf{J}}\end{equation} \tag{ 1 }$$

$$\begin{equation}\nabla \times {\mathbf{E}} = - \frac{{\partial {\mathbf{B}}}}{{\partial t}}.\end{equation} \tag{ 2 }$$

If we consider just linear magnetic materials ($\mu = {\text{const}}$), equations (1) and (2) yield the governing equation of the H formulation

$$\begin{equation}\nabla \times \left( {\rho \nabla \times {\mathbf{H}}} \right) = - \mu \frac{{\partial {\mathbf{H}}}}{{\partial t}}.\end{equation} \tag{ 3 }$$

Gauss's law $\nabla \cdot {\mathbf{B}} = 0$ is fulfilled by means of the election of the initial conditions, ${\mathbf{H}}\left( {t = 0} \right) = 0$, as explained in [15, 42].

Figure 2 depicts a two-dimensional (2D) planar model, where the x–y plane contains the cross-section of the superconductors. The subdomains ${{{\Omega }}_{{\text{sc}}}}$, ${{{\Omega }}_{\text{n}}}$, and ${{{\Omega }}_{{\text{sm}}}}$ represent the superconductor, normal conductor and surrounding mediums, respectively. The surrounding medium subdomain includes the insulating materials and the cryogenic liquid. In the 2D planar model, ${\mathbf{H}}$ has two non-zero components $({H_x},{\kern1.5pt}{H_y})$, while${\kern1.5pt}{\mathbf{E}}$ and ${\mathbf{J}}$ have just one non-zero component each $\left( {{E_z},{\kern1.5pt}{J_z}} \right)$.

Zoom In Zoom Out Reset image size

The transport currents ${I_k}$ in each conductor are imposed by means of integral constraints. One constraint is required for each conductor, as explained in [15, 42].

The selection of the elements used in the FEM discretization plays also an important role on the accuracy and computational speed of the numerical model. In the case of the H formulation, several arguments are presented in [15, 42, 52] showing the advantages of the first-order edge elements over other kind of elements.

The T-A formulation, as described in [47, 48], relies on the primary assumption that the thin superconducting layers of the HTS tapes can be modeled as one dimensional (1D) objects in a 2D model. The infinitely thin approximation is meaningful when dealing with superconductors wires having large aspect ratio (width/thickness), like the 2G HTS tapes, where this ratio is in the range of ${10^4}$ [1]. The T-A formulation requires the implementation of both the T and the A formulations, and both state variables ${\mathbf{T}}$ and ${\mathbf{A}}$, current and magnetic vector potentials, are evaluated.

Here, a 2D planar geometry is assumed, as illustrated in figure 3. The bounded universe is made of 1D superconducting layers and the surrounding medium ${{{\Omega }}_{{\text{sm}}}}$. The normal conductor layers of the HTS tapes are not considered. It is assumed that the current only flows through the superconducting layers, and the surrounding medium is assumed to be non-conductive. The current vector potential ${\mathbf{T}}$ is exclusively defined over the superconducting layers, while the magnetic vector potential ${\mathbf{A}}$ is defined over the entire bounded universe.

Zoom In Zoom Out Reset image size

Faraday's law (equation (2)) and the definition of the current vector potential $({\mathbf{J}} = \nabla \times {\mathbf{T}}$) are used to derive the governing equation of the T formulation

$$\begin{equation}\nabla \times \rho \nabla \times {\mathbf{T}} = - \frac{{\partial {\mathbf{B}}}}{{\partial t}}.\end{equation} \tag{ 4 }$$

In the 2D case depicted in figure 3, as long as the thickness of the superconducting layer can be neglected, ${\mathbf{J}}$ and ${\mathbf{T}}$ have only one non-zero component $\left( {{J_z},{\kern1.5pt}{T_y}} \right)$, and equation (4) is simplified as follows:

$$\begin{equation}\frac{\partial }{{\partial x}}\left( {\rho \frac{{\partial {T_y}}}{{\partial x}}} \right) = \frac{{\partial {B_y}}}{{\partial t}}.\end{equation} \tag{ 5 }$$

The transport currents in each tape ${I_k}$ are imposed by setting the boundary conditions for ${T_y}$. The values of ${T_y}$ at the edges of the 1D layer (${T_1}$ and ${T_2}$) must fulfill the following relation

$$\begin{equation}{I_k} = \left( {{T_1} - {T_2}} \right)\delta \end{equation} \tag{ 6 }$$

where $\delta$ is the real thickness of the HTS layer. Usually, ${T_1}$ or ${T_2}$ is set to zero, and the other value is computed by means of equation (6).

The component of B perpendicular to the superconducting layer ${B_y}$, required to compute ${T_y}$ in equation (5), is obtained by calculating ${\mathbf{A}}$. Ampère's law (equation (1)) and the definition of the magnetic vector potential ${\mathbf{A}}$ $\left( {{\mathbf{B}} = \nabla \times {\mathbf{A}}} \right)$ are used to derive the governing equation of the A formulation

$$\begin{equation}\nabla \times \nabla \times {\mathbf{A}} = \mu {\mathbf{J}}.\end{equation} \tag{ 7 }$$

In 2D cases, ${A_z}$ is the only non-zero component of ${\mathbf{A}}$, therefore equation (7) is simplified to ${\nabla ^2}{A_z}{ } = 0$. At first glance, equation (7) should be simplified to ${\nabla ^2}{A_z}{ } = - \mu {J_z}$. However, as the current flows only through the 1D superconducting layers, ${J_z}$ is zero all over the bounded universe, and the current is imposed by means of boundary conditions.

Thus, in order to couple ${J_z} = \partial {T_y}/\partial x$ with the A formulation, the surface current density ${\mathbf{K}}$ is introduced as

$$\begin{equation}\mathbf{{{K}}} = \delta \mathbf{{{J}}}.\end{equation} \tag{ 8 }$$

In the 2D case depicted in figure 3, ${\mathbf{J}} = \left( {0,{\kern1.5pt}0,{\kern1.5pt}{J_z}} \right)$.

As previously mentioned, ${\mathbf{K}}$ is then imposed into the A formulation as an external surface current density by means of boundary conditions of the form

$$\begin{equation}\boldsymbol{{{\hat n}}} \times \left( {{{\mathbf{H}}_1} - {{\mathbf{H}}_2}} \right) = {\mathbf{K}}\end{equation} \tag{ 9 }$$

where $\boldsymbol{{{\hat n}}}$ is the unit vector normal to the tape, and ${{\mathbf{H}}_1}$ and ${{\mathbf{H}}_2}$ are the magnetic field strength vectors above and below the HTS layer, respectively.

As in the case of the H formulation, the selection of the elements used in the FEM matters. Two kind of elements are required, Lagrange second-order elements are used to approximate ${\mathbf{A}}$ and Lagrange first-order elements for ${\mathbf{T}}$. These specific choices avoid spurious solutions, as justified in [50].

The case study used in this work is the same racetrack coil presented in [13, 50]. This coil has ten pancakes, each consisting of 200 turns, bringing a total number of turns equal to 2000. The transport current in each turn is the same, because they are connected in series. The geometric parameters of the coil are recalled in table 1. A unit cell is defined as the rectangular region occupied by a tape and its immediate surrounding medium. It is considered that the winding of the coil is even, thus all the unit cells have the same dimensions. The symmetry of the coil allows modeling just one quarter of the coil's cross-section. Therefore, it is possible to consider only five stacks, each consisting of 100 turns in a planar 2D geometry. The coil, its cross-section and the modeled section are depicted in figure 4.

Zoom In Zoom Out Reset image size

Table 1. Case study geometric parameters.

Parameter	Value
Number of pancakes	10
Turns per pancake	200
Unit cell width	4.45 mm
Unit cell thickness	293 $\mu {\text{m}}$
HTS layer width	4 mm
HTS layer thickness $\delta$	1 $\mu {\text{m}}$

The constitutive $E - J$ relation of the HTS material is modeled by the power-law [53], therefore the resistivity of the superconducting subdomains ${{{\Omega }}_{{\text{sc}}}}$ is given by

$$\begin{equation}{\rho _{{\text{HTS}}}} = \frac{{{E_{\text{c}}}}}{{{J_{\text{c}}}}}{\left| {\frac{{\mathbf{J}}}{{{J_{\text{c}}}}}} \right|^{n - 1}}.\end{equation} \tag{ 10 }$$

The Kim-like model [14, 54] is used to describe the anisotropic behavior of the HTS tapes, so that ${J_{\text{c}}}$ is defined by

$$\begin{equation}{J_{\text{c}}}\left( {\mathbf{B}} \right) = \frac{{{J_{{\text{c}}0}}}}{{{{\left( {1 + \frac{{\sqrt {{k^2}B_\parallel ^2 + B_ \bot ^2} }}{{{B_0}}}} \right)}^\alpha }}}\end{equation} \tag{ 11 }$$

where ${B_ \bot }$ and ${B_\parallel }$ are the magnetic flux density components perpendicular and parallel to the wide surface of the tape, respectively. The parameters of equations (10) and (11) are summarized in table 2.

Table 2. Case study electromagnetic parameters.

Parameter	Value
${E_{\text{c}}}$	$1 \times {10^{ - 4}}{\kern1.5pt}{\text{V}}{{\text{m}}^{ - 1}}$
$n$	$38$
${J_{{\text{c}}0}}$	$2.8 \times {10^{10}}{\text{ A}}{{\text{m}}^{ - 2}}$
${B_0}$	$0.04265{\text{ T}}$
$k$	$0.29515{\kern1.5pt}$
$\alpha$	$0.7$

The H formulation model that considers in detail each individual tape, presented in [13, 50], is used in this article to validate the rest of the models, and it is hereinafter called reference model.

The HTS tapes are composed by one layer of superconductor and different layers of normal conductors e.g. copper, silver, substrate [1]. The resistivity of the superconducting layer is several orders of magnitude lower than the resistivity of the other normal conductor layers forming part of the HTS tapes [15]. Therefore, the reference model does not include any normal conductor subdomain ${{{\Omega }}_{\text{n}}}$. The resistivity of the surrounding medium subdomain ${{{\Omega }}_{{\text{sm}}}}$ is considered to be ${\rho _{{\text{sm}}}} = 1~{{ \Omega \text{m}}}$ [15]. No magnetic materials are considered, then the permeabilities of the superconducting ${{{\Omega }}_{{\text{sc}}}}$ and surrounding medium ${{{\Omega }}_{{\text{sm}}}}$ subdomains are equal to the permeability of vacuum ${\mu _0}$.

Figure 4 depicts the geometry of the reference model including the numbering of pancakes and tapes. The mesh of the unit cells is structured with one element along the HTS layers' thickness and 100 elements along their width. An increasing number of elements towards the edges of the tapes allows increasing the accuracy of the ${\mathbf{J}}$ distribution in the regions where the magnetic field penetrates the tapes [13].

An assessment of the number of elements along the tapes' width was presented in [50]. The results demonstrated that, for the test conditions (transport current of 11 ${\text{A}}$, and 50${\text{ Hz}})$, the compromise between accuracy and computation time is fulfilled with 60 elements. The H formulation model that considers in detail each tape of the system and uses 60 elements along the tapes' width is hereinafter called H full model. The mesh of the unit cells of the H full model is also presented in figure 4. Accordingly, the difference between the reference and the H full models is the distribution and the number of elements along the tapes' width, 100 elements with a non-uniform distribution for the reference model, and 60 elements with a uniform distribution for the H full model. In addition, throughout the rest of this work, it is assumed that all models have 60 elements along the tapes' width, with the exception of the reference model. The presentation of the H full model with 60 elements is relevant to compare the computation times and decide which reductions are due to the number of elements and which reductions are due to the choice of the strategy.

The T-A full model uses the T-A formulation, and as well as the reference and H full models, considers in detail all the tapes. Specifically, in the case of the T-A formulation models, 'considers in detail all the tapes' means that the current vector potential ${\mathbf{T}}$ is computed along every single-tape. The mesh of the unit cells is also structured as shown in figure 4. In this case, the HTS layers have no thickness and the mesh is made of 1D elements uniformly distributed along the HTS layers' width.

The reference, H full and T-A full models were simulated for one cycle of a sinusoidal transport current with an amplitude of 11 ${\text{A}}$, and a frequency of 50${\text{ Hz}}$. The value 11 ${\text{A}}$ was chosen because at the peak of the cycle the tape 1 of pancake 5 is completely penetrated by the current density. The simulation results are compiled in figure 5, in tabular format. The first column contains the normalized current density ${{\mathbf{J}}_{\text{n}}} = {\mathbf{J}}/{J_{{\text{c }}}}$. The magnitude of the magnetic flux density $\left| {\mathbf{B}} \right|$ is presented in the second column. The plots of these first two columns show the results at the first negative peak of the transport current $t = 15{\text{ ms}}$. It can be seen that both the ${{\mathbf{J}}_{\text{n}}}$ and $\left| {\mathbf{B}} \right|$ plots are indistinguishable to the naked eye between the different formulations.

Zoom In Zoom Out Reset image size

The third column of figure 5 contains the average losses plots. The x-axis in these plots represents the tapes' number. There are five curves in each plot, one for each pancake. The numbering of the tapes and pancakes follows the order presented in figure 4. The losses estimated with the H and T-A full models are very similar to those estimated with the reference model, but there are visible differences, particularly in the first two pancakes. Due to the higher current penetration, the losses in pancake 5 are almost three orders of magnitude larger than the losses in pancake 1. Although there are variations in the losses at the end of the pancakes, the losses in a given pancake remain within the same order of magnitude.

The quantitative comparison of the models is carried out by calculating the relative error of the average losses, the coefficient of determination ${R^2}$ of the ${\mathbf{J}}$ distributions, and the normalized computation time. These data are listed in the fourth column of figure 5.

The average power losses are obtained using data of the second half of the cycle, as follows:

$$\begin{equation}{P_{av}} = \frac{2}{\tau }\mathop \int \limits_{\tau /2}^\tau \mathop \iint \limits_{{{{\Omega }}_{sc}}}^{} {\mathbf{E}} \cdot {\mathbf{J}}{ }\ dsdt\end{equation} \tag{ 12 }$$

where $\tau$ is the period of the sinusoidal cycle, and ${{{\Omega }}_{{\text{sc}}}}$ are the superconducting subdomains. In the 2D models addresses in this paper, the integral over ${{{\Omega }}_{{\text{sc}}}}$ is a surface integral, then the double integral symbol is preferred. The relative error of the average losses, expressed in percent, is defined as

$$\begin{equation}{\text{er}}_{\text{P}} = \frac{({P_{{\text{M}}\_{\text{av}}}} - {P_{{\text{R}}\_{\text{av}}}})}{P_{{\text{R}}\_{\text{av}}}} \times 100\% \end{equation} \tag{ 13 }$$

where ${P_{{\text{R}}\_{\text{av}}}}$ and ${P_{{\text{M}}\_{\text{av}}}}$ are the average losses computed with the reference model and the model that is being compared, respectively. The previous definition allows for positive and negative errors to be computed. Therefore, it will be possible to know if a model overestimates (${\text{e}}{{\text{r}}_{\text{P}}} > 0$) or underestimates (${\text{e}}{{\text{r}}_{\text{P}}} < 0$) the losses. For the test conditions, ${P_{{\text{R}}\_{\text{av}}}} = 127.24{\kern1.5pt}{\text{W}}{\kern1.5pt}{{\text{m}}^{ - 1}}$.

Unlike the average losses which are scalars, the ${\mathbf{J}}$ distributions are multivariable functions. The coefficient of determination is a widely used metric to evaluate the goodness of the fit [55], here it is used to compare the ${\mathbf{J}}$ distributions, and is defined as

$$\begin{equation}{R^2} = 1 - \frac{{\mathop \smallint \nolimits_0^\tau \mathop \iint \nolimits_{{{{\Omega }}_{sc}}}^{} {{\left( {{{\mathbf{J}}_R} - {{\mathbf{J}}_M}} \right)}^2}dsdt{\kern1.5pt}}}{{\mathop \smallint \nolimits_0^\tau \mathop \iint \nolimits_{{{{\Omega }}_{sc}}}^{} {{\left( {{{\mathbf{J}}_R} - \overline {{{\mathbf{J}}_R}} } \right)}^2}dsdt{\kern1.5pt}}}\end{equation} \tag{ 14 }$$

where ${{\mathbf{J}}_{\text{R}}}$ and ${{\mathbf{J}}_{\text{M}}}$ are the ${\mathbf{J}}$ distributions computed with the reference and tested models, respectively. $\overline {{{\mathbf{J}}_{\text{R}}}}$ is the mean value of ${{\mathbf{J}}_{\text{R}}}$. It must be remembered that ${R^2} = 1$ means a perfect matching between ${{\mathbf{J}}_{\text{R}}}$ and ${{\mathbf{J}}_{\text{M}}}$. The coefficient ${R^2}$ has an advantage over the error ${\text{e}}{{\text{r}}_{\text{P}}}$. The averaging nature of ${\text{e}}{{\text{r}}_{\text{P}}}$ tends to hide local and instantaneous errors, e.g. an instantaneous excess in the losses may be compensated by another instantaneous deficit; on the contrary, these same errors have a cumulative effect in ${R^2}$.

The normalized computation time is defined as

$$\begin{equation}\overline {{\text{ct}}} = \frac{{{\text{c}}{{\text{t}}_{\text{M}}}}}{{{\text{c}}{{\text{t}}_{\text{R}}}}}{\kern1.5pt} \times {\kern1.5pt}100{{{\kern1.5pt}\% }}\end{equation} \tag{ 15 }$$

where ${\text{c}}{{\text{t}}_{\text{R}}}$ and ${\text{c}}{{\text{t}}_{\text{M}}}$ are the computation times required by the reference model and the model that is being compared, respectively. For the test conditions, ${\text{c}}{{\text{t}}_{\text{R}}} = 31{\text{ h }}32{\text{ min}}$.

The results of the fourth column in figure 5 show that the accuracy of the H and T-A full models is satisfactory compared to the reference model. The errors ${\text{e}}{{\text{r}}_{\text{P}}}$ are lower than $1.7{\kern1.5pt}\%$, and coefficients ${R^2}$ are larger than 0.98. The computation time are 17 h 36 min ($\overline {{\text{ct}}} = 75.8{ }\%$) and 3 h 14 min ($\overline {{\text{ct}}} = 10.25{ }\%$), for the H and T-A full models, respectively. These values, more specifically the normalized computation times, demonstrate that the T-A formulation allows building more efficient full models in term of computation time keeping a fairly good accuracy.

The homogenization method consists in modeling the stacks of HTS tapes as homogeneous bulks. When the stack is transformed into an anisotropic bulk its geometric features are 'washed out'. This process is depicted in figure 6. The model should include additional features that allow the electromagnetic behavior of the homogeneous bulk to resemble that of the original stack.

Zoom In Zoom Out Reset image size

The non-superconducting materials forming part of the stack have resistivity values several orders of magnitude larger than those of the HTS material, hence only the HTS material resistivity is considered in the homogenization process. The resistivity of the bulk is derived from equations (10) and (11). But the ${J_{\text{c}}}$ value must be replaced by a homogenized critical current density ${J_{{\text{ch}}}}$, defined as

$$\begin{equation}{J_{{\text{ch}}}} = {f_{{\text{HTS}}}}{ }\ {J_{\text{c}}}\end{equation} \tag{ 16 }$$

where

$$\begin{equation}{f_{{\text{HTS}}}} = \frac{\delta }{\Delta }\end{equation} \tag{ 17 }$$

$\delta$ is the real thickness of the HTS layer, and ${{\Delta }}$ is the thickness of the unit cell. Thus, the superconducting properties of the HTS tapes are diluted in the cross-section of the bulk.

In the H full model, it is necessary to add one integral constraint per tape to impose the desired transport current. In the H homogenous model, the bulk subdomains ${{{\Omega }}_{\text{h}}}$ are further subdivided into bulk's subsets ${{{\Omega }}_{{\text{h}}\_{\text{i}}}}$, as depicted in figure 6. Then, one constraint is necessary for each bulk's subset, instead of one constraint per tape. The transport current imposed in a given subset is ${I_{{\text{tr}}}} = h{I_k}$, where ${I_k}$ is the transport current in the original tapes and $h$ is the number of tapes covered by the subset. The losses are computed in each bulk's subset. The losses in each subset are divided by the number of tapes included in each subset, then it is possible to approximate the losses along the stack. A detailed description of the H homogenous strategy can be consulted in [15].

The H homogenous model of the case study considers five bulks, one for each pancake. Each bulk is subdivided into six subsets, as shown in figure 7. The subsets in the upper part of the pancakes have a larger aspect ratio than the ones closer to the symmetry plane. Such kind of distributions has been recommended in [13, 15]. The mesh of the bulks is structured and considers one element along the subset's thickness and 60 elements along the tapes' width, as depicted in figure 7.

Zoom In Zoom Out Reset image size

Unlike the homogenization and the multi-scaling methods, the densification method is an original contribution of this article. The densification method addresses the analysis of stacks by means of a reduced number of tapes, called densified tapes. The densified tapes merge a given number of tapes into a single tape. In the densification process, the densified tapes preserve their original geometry and concentrate the transport current of their surrounding tapes, while the surrounding tapes are erased.

The idea of the densification, as in the case of the homogenization, is to build models with a smaller number of elements. Nevertheless, the reduced number of elements should not compromise the accuracy of the models. In the densification the number of elements is reduced by means of a reduced number of densified tapes. In the homogenization the electromagnetic behavior of the original stack is preserved by means of the distribution of the transport current all over the bulks, here this requirement is met by means of the concentration of the transport current in the densified tapes.

As in the previous models, the densified model does not include the normal conductors forming part of the HTS tapes. The resistivity of superconducting subdomains of the densified tapes is derived from equations (10) and (11), where the ${J_{\text{c}}}$ value is replaced by a densified critical current density ${J_{{\text{cd}}}}$ defined as

$$\begin{equation}{J_{{\text{cd}}}} = d{J_{\text{c}}}\end{equation} \tag{ 18 }$$

where $d$ is the number of tapes merged into a single densified tape.

The transport current in the densified tapes is ${I_{{\text{tr}}}} = d{I_k}$, where ${I_k}$ is the transport current in the original non-densified tapes. It is necessary to add one integral constraint per densified tape with the proper transport current ${I_{{\text{tr}}}}$. Then, as in the case of the H homogenous model, the number of constraints is reduced.

The densification process is depicted in figure 8. In this example, a given densified tape is built out of three tapes, labeled $r - 1$, $r$ and $r + 1$, therefore $d = 3$. The densified tape is located at the position of the original tape $r$. It is not necessary for $d$ to be an integer. The parameter $d$ may be equal to other real positive number. For instance, a stack made of five tapes can be modeled by means of two densified tapes. In this case, the densified tapes may merge three and two tapes, respectively. In another possible scenario, the densified tapes may merge 2.5 tapes, then the parameter $d$ should be $d = 2.5$ for both densified tapes.

Zoom In Zoom Out Reset image size

Once the ${\mathbf{J}}$ distribution is computed, the losses can be calculated in the densified tapes. The losses in the densified tapes are divided by their corresponding $d$, and these values are used to interpolate the losses in each tape of the original stack.

The accuracy of the densified models may be degraded due to the nature of the densified tapes, larger self-fields and larger distances between tapes than in the full models. Therefore, the number and the position of the densified tapes must be carefully assessed. In order to build a successful H densified model of the case study, different sets of densified tapes were tried. According to our heuristic criterion, the compromise between accuracy and computation time is fulfilled with a set of 31 densified tapes per pancake. The geometry of the model and the position of the densified tapes are depicted in figure 9. The first 21 densified tapes merge four tapes each ($d = 4$). For the following four tapes, the parameter $d$ takes the values {3, 3, 2, 2}, respectively. Finally, for the upper six tapes $d = 1$, meaning that these six tapes are not densified. The distribution of densified tapes is denser at the upper part of the pancake than at the bottom to achieve the required accuracy in the regions where larger variations in the ${\mathbf{J}}$ distributions arise.

Zoom In Zoom Out Reset image size

The idea of the multi-scaling method is to break up the model into several smaller models. In this way, it is possible to reduce the size of the problem by analyzing in detail a subset of significant tapes called analyzed tapes.

The multi-scale models, as described in [13], are formed by two 2D submodels. The first submodel is an A formulation magnetostatic model of the full coil including all the tapes with their actual geometry. This submodel, called coil submodel, does not consider any superconducting properties, hence the results depend on a predefined ${{\mathbf{J}}_0}$ distribution. The second submodel, called single-tape submodel, is an H formulation model of a unit cell containing just one tape. The single-tape submodel does not consider the normal conductor layers of the HTS tape, and the HTS layer is considered with its actual thickness. Both submodels are depicted in figure 10.

Zoom In Zoom Out Reset image size

The computational process is carried out in two steps. The first step is to use the coil submodel to estimate the background magnetic field strength ${\mathbf{H}}$ all across the bounded universe. Then, the ${\mathbf{H}}$ field along the boundary of the unit cells of the analyzed tapes is exported to the single-tape submodel as a time-dependent Dirichlet boundary conditions. The second step of the computational process is the use of the single-tape submodel. In this second step, the losses in all the analyzed tapes are calculated. Finally, the losses in the non-analyzed tapes are obtained by interpolation.

Breaking up the model into several smaller models not only reduces the computational burden, but also allows the parallelization of the problem, further reducing the computation time. A detailed description of the H multi-scale strategy can be consulted in [13].

The H multi-scale model of the case study uses 6 analyzed tapes in each pancake, 30 analyzed tapes in total. The distribution of the analyzed tapes is selected to be analogous to the distribution bulk's subsets in the H homogenous model. This means that the positions of the analyzed tapes correspond to the center of each bulk's subset. The set of analyzed tapes in each pancake is {25, 66, 88, 96, 99, 100}. The position of the analyzed tapes is shown in figure 11(a). The distribution of analyzed tapes also respects the directives proposed in [13, 49]. The mesh of the unit cells is structured and considers one element along the HTS layers' thickness and 60 elements along their width.

Zoom In Zoom Out Reset image size

The accuracy of the multi-scale models depends on the accuracy on the background magnetic field, which in turn depends on the predefined ${{\mathbf{J}}_0}$ distribution. The lack of knowledge of the predefined ${{\mathbf{J}}_0}$ distribution is the main limitation of the H multi-scale strategy [13]. To address this issue and preserve the capability to analyze the system by means of a reduced set of analyzed tapes, the iterative multi-scale strategy was proposed in [49]. The iterative multi-scale strategy is the iterative implementation of the multi-scale strategy. The iterative multi-scale strategy allows obtaining a new and more accurate dynamic solution at each iteration.

In the multi-scale models, the background magnetic field is transferred from the coil submodel to the single-tape submodel, whereas, in the iterative strategy, in addition to the background magnetic field, the current density is passed back from the single-tape submodel to the coil submodel. At the beginning of the procedure, the ${\mathbf{J}}$ distribution in every tape is supposed to be uniform, then the coil submodel is used to estimate the background magnetic field. The ${\mathbf{H}}$ field along the boundary of the analyzed tapes is exported as time-dependent boundary condition to the single-tape submodel. Now, the single-tape submodel is not only used to compute the losses but also the current density. An interpolation method is used to estimate the ${\mathbf{J}}$ distributions in the non-analyzed tapes. The new ${\mathbf{J}}$ distribution for all the tapes is exported to the coil submodel and a new background magnetic field is computed. The process is repeated to obtain more accurate estimations.

To exit from the iterative loop, the ${\mathbf{J}}$ distribution for all the tapes of the current iteration is compared with the distribution of the previous iteration until a convergence criterion is met. The error of the ${\mathbf{J}}$ distributions at the iteration $k$ is defined as

$$\begin{equation}e{J_k} = \frac{{\sqrt {\mathop \smallint \nolimits_0^\tau \mathop \iint \nolimits_{{{{\Omega }}_{{\text{sc}}}}}^{} {{\left( {{{\mathbf{J}}_{k - 1}} - {{\mathbf{J}}_k}} \right)}^2}{ }dsdt} }}{{\sqrt {\mathop \smallint \nolimits_0^\tau \mathop \iint \nolimits_{{{{\Omega }}_{{\text{sc}}}}}^{} {{\left( {{{\mathbf{J}}_{k - 1}}} \right)}^2}{ }dsdt} }}\end{equation} \tag{ 19 }$$

where ${{\mathbf{J}}_{k - 1}}$ and ${{\mathbf{J}}_k}$ are the ${\mathbf{J}}$ distributions for the iteration $k - 1$ and $k$, respectively. If the error $e{J_k}$ is smaller than a user-predefined criterion $\varepsilon$, then the process is completed. A detailed description of the H iterative multi-scale strategy can be consulted in [49].

The H iterative multi-scale model uses the same set of analyzed tapes used by the H multi-scale model, therefore figure 11(a) also represents the geometry of coil submodel of the H iterative multi-scale model.

The iterative multi-scale strategy requires the interpolation of the ${\mathbf{J}}$ distributions, the linear and the inverse cumulative density function (ICDF) interpolation methods were used. The ICDF interpolation method [56] was adapted to interpolate ${\mathbf{J}}$ distributions in [49]. This method produces more realistic current density distributions, avoiding some issues produced by the conventional linear interpolation. The implementation of ICDF interpolation, as presented in [49], requires the computation of integrals, derivatives and inverse functions. In the H iterative multi-scale model, this method is implemented in a MATLAB script. The following simultaneous multi-scale models presented below were implemented in a single COMSOL model, and for convenience just the simpler linear interpolation was used. This is not a major drawback because the ICDF interpolation makes only a marginal contribution to the accuracy of the model [49].

As described above, the H multi-scale and the H iterative multi-scale strategies use two different submodels. The computation of the ${\mathbf{J}}$ distributions in the analyzed tapes using the single-tape submodel can only be performed after the computation of the background magnetic field using the coil submodel. Therefore, the computation of the ${\mathbf{J}}$ distributions and the background field is not carried out simultaneously.

In this section, a new strategy called simultaneous multi-scale is proposed. The simultaneous multi-scale strategy allows simultaneously solving the ${\mathbf{J}}$ distribution and the background magnetic field. The strategy relies on the possibility to include an additional contribution to the Ampère's law (equation (1)). This summand allows imposing an external current density ${{\mathbf{J}}_{\text{e}}}$ in the superconducting subdomains ${{{\Omega }}_{{\text{sc}}}}$ of the non-analyzed tapes, as follows

$$\begin{equation}\nabla \times {\mathbf{H}} = {\mathbf{J}} + {{\mathbf{J}}_{\text{e}}}.\end{equation} \tag{ 20 }$$

Faraday's law (equation (2)) and the constitutive relations of the materials are used to derive the governing equation of the H formulation, which is expressed as

$$\begin{equation}\nabla \times \left( {\rho \left( {\nabla \times {\mathbf{H}} - {{\mathbf{J}}_{\text{e}}}} \right)} \right) = - \mu \frac{{\partial {\mathbf{H}}}}{{\partial t}}.\end{equation} \tag{ 21 }$$

The external ${{\mathbf{J}}_{\text{e}}}$ in the superconducting subdomains ${{{\Omega }}_{{\text{sc}}}}$ of the analyzed tapes and in the surrounding medium subdomain ${{{\Omega }}_{{\text{sm}}}}$ is zero. The external ${{\mathbf{J}}_{\text{e}}}$ in the superconducting subdomains ${{{\Omega }}_{{\text{sc}}}}$ of the non-analyzed tapes is approximated by interpolating the ${\mathbf{J}}$ distributions of the analyzed tapes.

The resistivity in the superconducting subdomains ${{{\Omega }}_{{\text{sc}}}}$ of the analyzed tapes is defined by equations (10) and (11). The resistivity of the superconducting subdomains ${{{\Omega }}_{{\text{sc}}}}$ of the non-analyzed tapes is considered to be the resistivity of the surrounding medium, ${\rho _{{\text{sm}}}} = 1{{{\kern1.5pt}\Omega \text{m}}}$. This value is orders of magnitude larger than the resistivity of the superconducting subdomains, therefore the induced current density in the non-analyzed tapes has a negligible impact when compared with the external ${{\mathbf{J}}_{\text{e}}}$.

The H simultaneous multi-scale model of the case study considers the same set of 30 analyzed tapes of the H multi-scale model. The non-analyzed tapes in the H simultaneous multi-scale model preserve their original geometry. Hence, figure 11(a) also represents the geometry of the H simultaneous multi-scale model. It is possible to reduce the number of degrees of freedom (DOF) and the computational burden of the H simultaneous multi-scale model by means of the homogenization or the densification of the non-analyzed tapes. Therefore, two additional models are presented here: the H simultaneous multi-scale homogenous model and the H simultaneous multi-scale densified model.

In the H simultaneous multi-scale homogenous and H simultaneous multi-scale densified models not all the non-analyzed tapes are homogenized or densified. The non-analyzed tapes adjacent to the analyzed tapes keep their original shape. These non-homogenous/non-analyzed or non-densified/non-analyzed tapes are used to move the distortions in the magnetic field produced by the homogeneous or densified tapes away from the analyzed tapes. The geometries of the three H simultaneous multi-scale models are presented in figure 11.

To compare and validate the strategies described in this section, these models too were simulated for a sinusoidal transport current (11 A, 50 Hz). The results of the H homogenous, H densified, H multi-scale and H iterative multi-scale models are presented in figure 12. The results of the three H simultaneous multi-scale models are presented in figure 13. These figures have the format of figure 5. The first two columns show ${{\mathbf{J}}_{\text{n}}}$ and $\left| {\mathbf{B}} \right|$, both at $t = 15{\text{ ms}}$, the average losses are displayed in the third column, and the last column contains the quantitative data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first row in figure 12 shows the results of the H homogenous model. This model successfully reproduces the screening currents of the reference model. The ${{\mathbf{J}}_{\text{n}}}$ plot in this row shows how the individual tapes were replaced by the homogenous bulks. The $\left| {\mathbf{B}} \right|$ distribution presents a smoother profile due to the homogenous current densities. The accurate estimation provided by this model is confirmed by the following values: ${\text{e}}{{\text{r}}_{\text{P}}} = 1.28{ }\%$ and ${R^2} = 0.9221$. In this case the coefficient ${R^2}$ is computed using the rescaled current density (divided by ${f_{{\text{HTS}}}}$) and considering just the values at the positions of the original superconductor subdomains. The computation time required by the H homogenous model is 36 min 44 s ($\overline {{\text{ct}}} = 1.94{ }\%$).

The second row in figure 12 presents the results of the H densified model. Thicker lines are used to represent the densified tapes. In contrast to the homogenous case, here the $\left| {\mathbf{B}} \right|$ distribution has a jagged profile. This profile degrades the accuracy as can be observed in the values ${\text{e}}{{\text{r}}_{\text{P}}} = 6.67{ }\%$ and ${R^2} = 0.8549$. For the purpose of computing ${R^2}$, the ${\mathbf{J}}$ distributions of the densified tapes are rescaled (divided by $d$) and the ${\mathbf{J}}$ distributions in the removed tapes are approximated with linear interpolation. The computation time is reduced compared to the H full model, but this reduction is not as big as in the case of the H homogenous model, the normalized computation time is $\overline {{\text{ct}}} = 29.03{{\% }}$.

The results of the H multi-scale model are presented in the third row of figure 12. The predefined current density distribution ${{\mathbf{J}}_0}$ is uniform, as can be seen in the first entry of the third row. The uniform ${{\mathbf{J}}_0}$ does not contain any screening current, then it is not a good approximation of the reference ${\mathbf{J}}$ distribution. This fact is also reflected in the low coefficient ${R^2} = 0.0304$. Consequently, the magnetic flux density and the losses exhibit noticeable errors, especially at the upper portion of the pancakes. The red circles in the plot of the losses indicate the position of the analyzed tapes, this is also the case of all the multi-scale models. The losses error is ${\text{e}}{{\text{r}}_{\text{P}}} = - 21.7\%$, the negative sign indicates that the losses are underestimated. The computation time is 27 min 30 s ($\overline {{\text{ct}}} = 1.45\%$). This time is the summation of the time required to run the coil submodel one time (5 min) and the single-tape submodel 30 times, one for each analyzed tape (the average computation time of the single-tape submodel is 45 s).

The results of the H iterative multi-scale model are listed in the last row of figure 12. The convergence criterion is defined as $\varepsilon = 0.01$, this criterion is reached at the 7^th iteration regardless of the interpolation method (linear or ICDF) used to approximate the ${\mathbf{J}}$ distributions in the non-analyzed tapes. Figure 12 presents the results when the ICDF interpolation is used. The results when the linear interpolation is applied are visually indistinguishable from those obtained from the ICDF interpolation, then for clarity only the latter are shown. The error ${\text{e}}{{\text{r}}_{\text{P}}}$ when linear interpolation is used is $- 0.87{ }\%$, when ICDF interpolation is used ${\text{e}}{{\text{r}}_{\text{P}}} = - 0.56\%$. The coefficient ${R^2}$ takes values 0.9796 and 0.9803, for the linear and ICDF interpolations, respectively. The accuracy is marginally better with the ICDF interpolation. The computation time is 3 h 21 min ($\overline {{\text{ct}}} = 10.51{{\% }}$) with linear interpolation, and 3 h 17 min ($\overline {{\text{ct}}} = 10.33{{\% }}$) with ICDF interpolation. These times are approximately seven times the computation time of the H multi-scale model.

The rows of figure 13 compiles the results of the H simultaneous multi-scale, H simultaneous multi-scale homogenous, and H simultaneous multi-scale densified models, respectively. The plots in the first column allow the observation of the different approaches in which the non-analyzed tapes are modeled: tapes with their original geometry, homogenous bulks or densified tapes. The $\left| {\mathbf{B}} \right|$ distribution of the H simultaneous multi-scale is visually indistinguishable from that of the reference model. Moreover, it is easy to find similarities between the distortions in the $\left| {\mathbf{B}} \right|$ distributions of the H homogenous and H simultaneous multi-scale homogenous models, and between the distortions of the H densified and H simultaneous multi-scale densified models. The three models have acceptable and similar accuracies, as demonstrated by the ${\text{e}}{{\text{r}}_{\text{P}}}$ values lower than $1.6{{\% }}$, and the ${R^2}$ values greater than 0.98. The advantage of the homogenization and densification of the non-analyzed tapes is clearly observed in the computation times. The computation times of the H simultaneous multi-scale is 16 h 56 min ($\overline {{\text{ct}}} = 53.7{{\% }}$). This computation time is approximately two times the computation time of the other two models.

The manner in which the homogenization is used in conjunction with the T-A formulation is depicted in figure 14. The magnetic vector potential ${\mathbf{A}}$ is defined all over the entire bounded universe. The stacks of 1D HTS layers are transformed into 2D bulks, and the potential ${\mathbf{T}}$ is exclusively defined inside the bulk. For the purpose of computing ${\mathbf{T}}$, the influence of the component of ${\mathbf{B}}$ parallel to the surface of the tapes is not considered, in figure 14 this parallel component is ${B_x}$. Therefore, from equation (4), it follows that T has only one non-zero component (${T_y}$ in the case of figure 14), which is defined by means of equation (5). Nevertheless, the parallel component ${B_x}$ influences the calculations by means of the definition of ${J_{\text{c}}}$, see equation (11).

Zoom In Zoom Out Reset image size

The densely packed stack (homogenous bulk) must carry a transport current that is the summation of the transport current of the tapes making up the original stack. To impose such transport current, it is necessary to use the values ${T_1}$ and ${T_2}$, defined in equation (6), as Dirichlet boundary conditions along the edges of the bulks perpendicular to the tapes. To compensate for the fact that the current density is computed inside the homogenous bulk, a new homogenized current density ${J_{{\text{zh}}}}$ is defined as

$$\begin{equation}{J_{{\text{zh}}}} = {f_{{\text{HTS}}}}{J_z}\end{equation} \tag{ 22 }$$

where ${f_{{\text{HTS}}}}$ is the ratio defined in equation (17). ${J_{{\text{zh}}}}$ is imposed as a source into the A formulation, then equation (7) is transformed into,

$$\begin{equation}{\nabla ^2}{A_z} = - \mu {J_{{\text{zh}}}}.\end{equation} \tag{ 23 }$$

For the purpose of computing T, the resistivity of the bulk subdomains is considered to be the resistivity of the superconducting material, defined by equations (10) and (11). The losses are computed by integrating the local losses along the lines parallel to the HTS layers at the center of each bulk's subset, then the losses along the rest of the tapes are approximated by interpolation. A more detailed description of the T-A homogenous strategy can be consulted in [50].

The T-A homogenous model of the case study, similarly to the H homogenous model, considers five bulks. Here, the integral constraints (one per bulk's subset) are not required to impose the transport currents. These subsets are used to define the distribution of the elements along the bulk's height. The mesh of the bulks in the T-A homogeneous model is structured and considers one element along the subsets thickness and 60 elements along the tapes' width. Figure 7 also represents the geometry of the T-A homogenous model.

The densification method can also be used in conjunction with the T-A formulation. Here too, the idea is to model the original stack by means of a smaller number of densified tapes.

In a T-A densified model the densified tapes are one-dimensional objects along which the ${\mathbf{J}}$ distributions are computed by means of the current density potential ${\mathbf{T}}$, and its only non-zero component ${T_y}$, as given in equation (5). The resistivity of the HTS material is defined by means of equations (10) and (11). Unlike the H densified strategy, here the critical current density is not modified. In the T-A formulation, the surface current density ${\mathbf{K}} = \delta {\mathbf{J}}$ is imposed into the A formulation, see equations (8) and (9). Thus, to take into account the densification of the tapes and to preserve the electromagnetic behavior of the original stack, the surface current density to be imposed is now scaled so that

$$\begin{equation}{{\mathbf{K}}_d} = d{\mathbf{K}}\end{equation} \tag{ 24 }$$

where $d$ is the number of tapes merged into a single densified tape.

Also the T-A densification process depicted in figure 8 involves two steps. The first step is to remove the tapes $r -$1 and $r + 1$. Second, the magnetic effect of tape $r$ is forced to be three times larger than the magnetic effect of the same tape in the T-A full model. The transport current in the densified tape remains the same as the transport current in the original tapes, see equation (6). The magnetic effect of the densified tapes is incremented by means of the parameter $d$ in equation (24). The losses can be calculated in the densified tapes, these losses are considered as the losses produced by the tapes of the original stack in the position of the densified tapes. The losses in the removed tapes are approximated by interpolation.

The T-A densified model of the case study considers the set of densified tapes of the H densified model, 31 densified tapes per pancake. Then, figure 9 also represents the geometry of the T-A densified model.

The T-A simultaneous multi-scale strategy, as well as the H simultaneous multi-scale strategy, does not require two different submodels. The computation of the background magnetic field and the ${\mathbf{J}}$ distribution are carried out in a single model based on the T-A formulation.

In the T-A full models, the current vector potential ${\mathbf{T}}$ is defined over all the tapes. In the present approach, ${\mathbf{T}}$ is defined only along the analyzed tapes. The ${\mathbf{J}}$ distribution along the analyzed tape is obtained by calculating ${\mathbf{T}},$ see equation (5). The ${\mathbf{J}}$ distributions in the non-analyzed tapes are approximated by linear interpolation using the ${\mathbf{J}}$ distributions of the analyzed tapes. The magnetic potential ${\mathbf{A}}$ is defined over the entire bounded universe. The current density in both analyzed and non-analyzed tapes is multiplied by the thickness of the superconducting layer $\delta$ to obtain a surface current density ${\mathbf{K}}$ to be imposed into the A formulation, see equations (8) and (9).

As it was the case with the H simultaneous multi-scale models, the DOF can be reduced by means of the homogenization or densification of the non-analyzed tapes. Therefore, three T-A simultaneous multi-scale models are presented. The difference between these models is the treatment of the non-analyzed tapes. The three models use the same set of six analyzed tapes per pancake used in the H multi-scale models. The first model, called T-A simultaneous multi-scale, considers the non-analyzed tapes with their original number and geometry.

As justified in [50], the T-A formulation uses first order elements to approximate ${\mathbf{T}}$, and second order elements to approximate ${\mathbf{A}}$. If first order elements are used for both quantities, the computation time can be reduced [48], but this choice produces undesired spurious oscillations in the ${\mathbf{J}}$ distributions [50]. To increase the computational efficiency without compromising the accuracy, the unit cells of the analyzed tapes and their adjacent non-analyzed tapes use second order elements to approximate ${\mathbf{A}}$, while first order elements are used to approximate ${\mathbf{A}}$ throughout the rest of the system. Additionally, 30 elements are considered along most of the non-analyzed tapes, while 60 elements are considered in the analyzed tapes and their adjacent non-analyzed tapes. A more detailed description of the T-A simultaneous multi-scale strategy can be consulted in [50].

The other two T-A multi-scale models modify the geometric description of the pancakes by means of the homogenization or densification of the non-analyzed tapes. These models are called T-A simultaneous multi-scale homogenous and T-A simultaneous multi-scale densified models, respectively. As it was done in section 4.5, the non-analyzed tapes adjacent to the analyzed tapes keep their original geometry. These tapes with their original geometry are used to move the distortions in the magnetic field away from the analyzed tapes. The geometries of the T-A multi-scale models correspond to the geometries of the H multi-scale models, and are shown in figure 11.

The T-A models were validated with the same operating conditions used in the previous sections. The results are compiled in figure 15, using the tabular format of figure 5.

Zoom In Zoom Out Reset image size

The first row in figure 15 shows the results of the T-A homogenous model. Due to the homogenization process, a smoother $\left| {\mathbf{B}} \right|$ distribution can be observed. In the loss plot of this row, the losses in the pancake 1 deviate from the reference results. The losses in the first pancake are two orders of magnitude lower than those of the pancake 5, therefore the deviation does not affect the global result. The total losses are slightly overestimated, the error is ${\text{e}}{{\text{r}}_{\text{P}}} = 0.71\%$. The computation time required by the T-A homogenous model is 14 min 51 s ($\overline {{\text{ct}}} = 0.78{{\% }}$).

The results of the T-A densified model are presented in the second row of figure 15. The densification process causes a more jagged $\left| {\mathbf{B}} \right|$ distribution. The coefficient ${R^2} = 0.8854$ demonstrates that, among the T-A models, the worst approximation of the current density is achieved by the T-A densified model. However, it remains an acceptable strategy for estimating the losses, as demonstrated by the value ${\text{e}}{{\text{r}}_{\text{P}}} = - 2.62\%$. In contrast to the H densified model, here the losses are slightly underestimated. The computation time of the T-A densified model is 61 min ($\overline {{\text{ct}}} = 3.22{{{\kern1.5pt}\% }}$), approximately one third of the computation time of the T-A full model.

The last three rows contain the results of the T-A simultaneous multi-scale models. The ${{\mathbf{J}}_{\text{n}}}$ and $\left| {\mathbf{B}} \right|$ distributions calculated with the T-A simultaneous multi-scale model is indistinguishable from the distributions of the reference model presented in figure 5. The $\left| {\mathbf{B}} \right|$ distributions of the T-A simultaneous multi-scale homogenous and densified models show the respective distortions produced by the homogenization and densification of the non-analyzed tapes. Accurate estimations are achieved with the three T-A simultaneous multi-scale models. In the three cases, the magnitude of the error ${\text{e}}{{\text{r}}_{\text{P}}}$ is less than $1{{\% }}$, and the coefficient ${R^2}$ is greater than $0.98$. The computation time of the T-A simultaneous multi-scale model ($\overline {{\text{ct}}} = 5.06{{\% }}$) is half that of the T-A full model ($\overline {{\text{ct}}} = 10.25{{\% }}$). The densification of the non-analyzed tapes further reduces the computation time by 30 min ($\overline {{\text{ct}}} = 3.46{{\% }}$). Conversely, the homogenization of the non-analyzed tapes produces a noticeable increment of the computation time ($\overline {{\text{ct}}} = 18.41{{\% }}$).