A numerical model to introduce students to AC loss calculation in superconductors

A numerical model implemented in the open-source finite-element method FreeFEM program is presented, with the aim of introducing students to the calculation of AC losses in superconductors. With this simple approach, students can learn about the critical state model used to describe the macroscopic electromagnetic behavior of superconductors and the importance of different factors influencing the AC losses of superconductors.


Introduction
Despite being materials with no electrical resistance, superconductors exhibit power loss dissipation when subjected to time-varying magnetic fields. This is a situation met often in several applications. For example, power cables carrying AC transport current or DC superconducting magnets being charged or discharged.
Such power dissipation can represent a significant burden, because of the low temperatures at which superconductors operate. For example, 1 W dissipated at cryogenic temperature has a much higher cost in terms of the energy necessary to remove the generated heat. On top of 2 Enrico Rizzo works now with a different employer. 1 Author to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. that, the efficiency of machines used to obtain cryogenic temperatures is very far from the ideal one [1]. In short, such power dissipation-which goes under the general name of 'AC loss'-can limit the efficiency and increase the cost of superconducting applications, thus limiting the possibility of their penetration in the market. Therefore, for the design of efficient superconducting applications, numerical models able to accurately predict the AC losses are extremely important [2].
In this contribution, we propose a numerical model aiming at introducing students to AC loss calculation in superconductors. The model is based on the finite-element method (FEM) and is implemented in the open-source program FreeFEM [3,4], which can be easily installed as a stand-alone application on different platforms. This enables students to run the scripts on their computers also outside of lectures and to review the problems at their own pace. The code is available on the HTS Modelling Workgroup website [5] and as supplemental material of this article.
The article is organized as follows. Section 2 concisely recalls the origin of AC losses in type-II superconductors and introduces the critical state model (CSM) for modelling the electrodynamics of 'hard' superconductors. Section 3 describes the implementation of the CSM into the program FreeFEM. Section 4 contains a series of exercises related to AC loss calculation in superconductors. Section 5 discusses the relation of the results found in the exercises to cases of practical interest. Finally, section 6 draws the main conclusions.

AC losses and critical state model
Superconductors used in large-scale applications belong to the group of type-II superconductors, which can maintain the superconducting state while being subjected to large magnetic fields. They operate in the so-called mixed state: the magnetic field penetrates the superconductor in a quantized form: each quantum carries a very small quantity of magnetic flux, approximately equal to 2.07 × 10 −15 Wb. These flux quanta are characterized by a normal zone and by superconducting screening currents, and they are also referred to as 'vortices'. The spatial gradient of magnetic flux inside the material gives rise to macroscopic currents (figure 1). A detailed description and a visual representation of the mechanism are available in [6].
The vortices can anchor (or 'pin') to imperfections of the materials, avoiding the movement that would cause energy dissipation. This pinning mechanism is what enables type-II superconductors to carry large currents in the presence of magnetic fields without dissipation. The materials used in applications have a strong pinning force and are called 'hard' superconductors. When the magnetic field varies with time, the magnetic flux quanta have to rearrange themselves in order to follow the time-varying excitation: as they move, they dissipate energy. This is the mechanism behind AC loss dissipation in type-II superconductors.
For estimating the AC losses, one needs to determine how the macroscopic currents evolve as a result of time-varying magnetic fields. This can be done by means of the critical state model (CSM), developed by Bean in the early 1960s [7,8]. The model assumes that the magnitude of the macroscopic current density is either null or equal to the so-called critical current density J c .
The evolution of the current density in a superconductor is determined by the magnetic history, according to the following principles: • no current flows in the regions that have not been previously reached by the flux vortices; • in the rest of superconductor, the flux vortices arrange with a gradient of density of magnitude J c . In mathematical terms, this is expressed as follows: The direction of the current density is determined by the magnetization history of the superconductor. An example is shown in figure 1 for an infinite superconducting slab subjected to a time-varying magnetic field parallel to its face, for which the current density and magnetic field profiles can be easily calculated. For this particular case, the slab is infinitely large in the y and z directions and the problem to solve is one-dimensional: one has to find the current and field profiles along the thickness of the slab (x direction). In addition, for this infinite geometry, the field at the vacuum-superconductor interface is known and equal to the applied one.
The general relation Based on (1), the current density is different from zero only where the magnetic field is different from zero and the gradient of the magnetic field is equal to ±J c , according to (3). Figure 1 shows the magnetic field profiles across the slab's thickness as the field is increased from zero to a maximum, reduced to zero, and then inverted to a negative value. One remarkable feature is that when the field is brought back to zero, some magnetic field is trapped in the superconductor. The profile of the trapped field depends on the amplitude of the field at which the slab has been exposed. With the CSM, the rate with which the field changes is irrelevant. The profiles of figure 1 depend only on the magnetization history of the sample, not on the rate with which the magnetization occurs.
From the current density and magnetic field distributions shown in figure 1, it is possible to calculate the cyclic losses of a superconducting slab subjected to an AC magnetic An infinite superconducting slab subjected to a magnetic field parallel to its face. Left: As the external field is increased from zero, it penetrates further inside the superconductor; the slope of the front is constant and equal to −J c , according to (3). Right: As the external field is decreased, oppositely sloped fronts (slope +J c ) appear near the surface of the slab; at t = t 3 , the external field is null, but some magnetic field is trapped inside the superconductor; at t = t 5 the external field is completely reversed with respect to t = t 2 and so are the field fronts inside the superconductor. Reproduced from [9]. CC BY 3.0.
field. Other geometries, such as superconducting wires with circular cross section, are more challenging and finding the current density and magnetic field distributions with pen-and-paper calculations can be very difficult. Numerical simulations provide an easier and much more flexible approach for finding such distributions. In the following section, we describe a possible implementation of the CSM into a finite-element method program. Such implementation will be used to calculate the current and magnetic field profiles and to evaluate the AC losses in cases of increasing complexity.

The Campbell model
In 2007, Prof. A M Campbell presented a numerical model for the solution of the critical state in superconductors using the magnetic vector potential as state variable [10]. The basic idea of the model is to solve where A is the magnetic vector potential and J is the critical current density. For the 2D problems considered in this paper, the equation becomes According to the CSM, the E-J relationship takes the form of the thick black line in figure 2: any electric field drives a current of magnitude J c . For the implementation in a numerical model, the sharp E-J relationship of the CSM needs to be smoothed. This can be done in different ways, here we use where erf is the error function [11] and E 0 is a parameter controlling the steepness of the curve.
Other functions based on exponential [10] or hyperbolic tangent [12] can be used. Note that, differently from the 'pure' CSM, with these approximations the current density can take values in the interval [−J c , J c ].
The electric field can be expressed as In this paper we consider individual tapes or multiple tapes in parallel, so the voltage gradient term ∇φ can be ignored. The combination of equations (5)-(7) allows solving general time-dependent problems and calculating the power dissipation as product J · E. The cyclic losses of an AC excitation (e.g. a sinusoidal external magnetic field) can be calculated by integrating J · E on the superconductor's cross section and averaging over the AC cycle. However, it is possible to demonstrate that in certain AC cases, the cyclic losses can be calculated as four times the losses dissipated in the first quarter of the AC cycle (starting from a virgin sample, i.e. with no trapped magnetic flux) [13][14][15][16]. In addition, as long as the penetration of the field inside the superconductor is monotonic, the time integration of the power dissipation from t = 0 to the peak can be performed in one step. In that case, the AC losses can be computed as where Ω is the superconductor's cross section, and J p and A p are, respectively, the current density and the magnetic vector potential at the peak of the AC cycle. From the practical point of view, this means that the cyclic AC losses can be computed with a single one-step simulation, solving for the magnetic vector potential at the peak of the AC cycle 3 .

Exercises
The model aims at calculating the current density and magnetic field distributions in the transverse cross sections of an infinitely long superconductor. In this respect, the computational domain is a two-dimensional space consisting of one air region in which one or more superconductor cross sections are immersed. The magnetic flux density has only the two components B x and B y , whereas the current density and magnetic vector potential have only the z component and are therefore scalar (see figure 3 for the definition of the coordinates). The current density distribution in the superconductor is calculated by solving the static equation where A r is a parameter resulting from the combination of E 0 and the time it takes to reach the peak of the AC excitation (for a sine at 50 Hz, this is 5 ms). In this work, we used A r = 1 × 10 −7 Wb m −1 . An external uniform magnetic field is applied by setting the boundary conditions for A on the outer air boundary. In 2D Cartesian coordinates, B = ∇ × A becomes B x = ∂A/∂y, B y = −∂A/∂y. For example, in order to obtain a magnetic field B 0 along y, it is sufficient to impose A = −B 0 x as boundary condition. In a similar way, a transport current is imposed by setting a constant value for the magnetic vector potential on the boundary. The actual value of the current can be calculated in the post-processing by integrating the current density J over the superconductor's cross section.
The solution is sought on a multi-region domain (air and superconductor) discretized by means of a mesh with triangular pattern. As test functions, second order polynomial functions are used.
The weak formulation of the problem is [ where Ω is the computational domain (air and superconductor) and v are second-order polynomials.
Due to the non-linearity of the problem, the achievement of a stable solution occurs by iteratively solving the weak formulation using the values computed in the previous iteration as input for the next one until the difference among two consecutive iterations drops below a certain value, in this case 1 × 10 −7 .
The cyclic losses can be computed by means of equation (8). In the exercises, we will refer to magnetization and transport losses. The former are caused by an externally applied AC magnetic field, the latter by an AC transport current.
The commented code for solving the exercises of sections 4.1 and 4.2 and a video tutorial are provided as supplemental material (https://stacks.iop.org/EJP/41/045203/mmedia). The code can also be used to solve the exercise with transport current (section 4.3) by simply changing the boundary conditions.

Round superconductor in external field
This exercise consists in calculating the magnetization losses of a round superconducting wire subjected to an external magnetic field. The superconductor is defined as an ellipse with semiaxes a = R √ η and a = R/ √ η, so that the area of the ellipse is always πab, independently of the value of η. This definition will be useful for the exercises of section 4. For now, we consider a circle and then set η = 1. We choose R = 1 mm and a critical current density J c = 1 × 10 8 A m −2 . The critical current of the wire is therefore about 314.15 A. Figure 4 shows the distributions of magnetic vector potential, magnetic flux density and current density in a round superconductor subjected to a magnetic field of 20 mT. 4 Note the bending of the magnetic vector potential lines around the superconductor (figure 4 (a)), which results in a magnetic flux density higher than the applied one ( figure 4(b)). In the superconductor, a non-zero current density is present only in the regions where there is magnetic flux density, and its value is ±J c , as expected according to the CSM ( figure 4(c)). Figure 5 shows the calculated cyclic losses as a function of the magnitude of the external field B 0 . The losses agree well with the approximate analytical expression by Carr [19]: where H 0 = B 0 /μ 0 and H p = 2RJ c /π. On a log-log scale like that of figure 5, one can observe a change of slope, from Q ≈ B 3 to Q ≈ B 1 , occurring when the field fully penetrates the superconductor. This exercise confirms the accuracy of AC loss calculation by means of equation (11) and a general important feature of the magnetization losses of superconductors: when plotted as a function of the applied field, initially they increase very rapidly, then, after the field fully penetrates the superconductor, more gradually.

Elliptical superconductor in external field
The first part of this exercise consists in repeating the previous exercise keeping the external field at a given value (and applied along the y direction) and varying the aspect ratio of the ellipse by changing the parameter η ( figure 6). It should be noted that even if the aspect ratio of the superconductor changes (because of η), both the amount (i.e. the area of the ellipse πab) and the 'quality' (i.e. the critical current density J c ) of the superconductor material remain the same. In other words, this exercise investigates solely the effects of a geometrical deformation of the superconductor's cross section on the AC losses. Figure 7 shows the magnetization AC losses of an elliptical superconductor for an applied magnetic field of 20 mT as a function of the parameter η. This parameter modifies the aspect ratio of the ellipse, and values of η > 1 correspond to progressively flatter ellipses 5 . The magnetic field is applied parallel to the ellipse's minor axis. As a consequence, as the ellipse become flatter, the magnetic field lines bend and concentrate more around the edges of the ellipse. The field there is progressively higher (an example for η = 10 is shown in figure 7) and this explains why the losses increase with η.
One could expect the AC losses to decrease for η < 1 (when the field becomes parallel to the ellipse's major axis), but this happens only to a certain extent. This is probably due to the fact that as η becomes very small, the ellipse becomes very elongated in the direction  parallel to the field. It is therefore very easy for the field to penetrate the superconductor and hence the losses do not decrease.
The second part of this exercise deals with another important factor affecting the values of the AC losses: the direction of the magnetic field. Figure 8 shows the AC losses for an elliptical superconductor with η = 10 subjected to a magnetic field of 20 mT of various orientations. When the field is parallel to the minor axis (θ = 0 • ), the field near the edges is greatly enhanced (see figure 7) and the losses are high. When the field is parallel to the major axis (θ = 90 • ), the distortion of the field is minimal, the field is almost equal to the applied one and the losses are low.
As it will be discussed later, the shape of the superconductor and the orientation of the magnetic field have important consequences for the design of superconducting applications.

Elliptical superconductor with transport current
This exercise consists in calculating the AC losses in a superconductor of elliptical cross section caused by an AC transport current of magnitude I 0 . A transport current can be imposed by applying a constant value of the magnetic vector potential on the boundary of the air domain. The value of the generated current is calculated in the post processing as integral of J over the superconductor's cross section. Figure 9 shows the distributions of magnetic vector potential, magnetic flux density and current density in a round superconductor carrying transport current. Similarly to the case of magnetization, a non-zero current density is present only in the regions where there is magnetic flux density. Figure 10 shows the transport AC losses as a function of the current magnitude for η values equal to 1 and 10. The results are rather similar. This is in accordance with the result of the analytical formula by Norris [13]: which predicts that the losses do not depend on the aspect ratio of the ellipse. This is a substantial difference with respect to the magnetization losses.

Stack of tapes in external field
Superconducting tapes are often assembled together in cable structures, in order to obtain conductors with higher current-carrying capacity. The simplest tape assembly is a stack of tapes. In this and the next subsections, we investigate the magnetization and transport AC losses of a stack of five tapes. We use elliptical tapes with η = 10. Figure 11 shows the current density distribution in a five-tape stack subjected to a magnetic field of 50 mT parallel to the ellipses' minor axes. Due to the small separation between the tapes (1 mm center-to-center), the tapes have a strong electromagnetic coupling. As a result, the stack behaves similarly to a unique larger conductor, with a width-to-thickness ratio close to unity. This is clearly visible from the shape of the current-free region, which resembles that of the circular superconductor shown in figure 4(c).
From figure 11 one can see that the top and bottom tapes have a larger current penetration. This is because they are exposed to a larger magnetic field. As a consequence, they are expected to have larger AC losses than the tapes in the center of the stack, which are partially 'protected' from the penetration of the magnetic field. This is confirmed from the data shown in table 1, which show the magnetization AC losses of the five-tape stack for different center-to-center Figure 10. Transport AC losses as a function of the magnitude of transport current I 0 for elliptical superconductor wires of different aspect ratio. Figure 11. Current density distribution in a five-tape stack subjected to a magnetic field of 50 mT, parallel to the minor axes of the ellipses. The center-to-center distance between the tapes is 1 mm. distances. As the distance increases, the tapes tend to behave more like individual ones, and the AC losses increase.
One can calculate the AC losses of an individual elliptical tape (with η = 10) subjected to the same field of 50 mT and find them to be 4.03 × 10 −2 J m −1 . Then, for reference, one can estimate the losses of five 'independent' tapes as five times that amount, i.e. 2.02 × 10 −1 J m −1 , and observe that the losses of the stack are much lower, especially when the tapes are closely packed. Figure 12. Current density distribution in a five-tape stack carrying a total transport current of 930 A. The center-to-center distance between the tapes is 1 mm. The current is not evenly distributed between the tapes.  One can conclude that assembling the tapes into a stack has a positive effect for the magnetization losses. This is due to the fact that, if the tapes are closely packed, the stack behaves as a monolithic conductor with a smaller width-to-thickness ratio, so that the concentration of the magnetic flux at the edges is reduced. As the distance between the tapes is increased, the magnetic flux penetrates between the tapes, which tend to behave as individual ones. The monolithic-conductor behavior gradually disappears and the losses increase.

Stack of tapes with transport current
In the case of a stack of tapes, the current is injected at the ends where the tapes are usually soldered together. Due to the different impedance of the current paths offered by the different tapes, the current is not equally divided between them, but tends to flow in the tapes situated at the top and the bottom of the stack. This case can be easily considered with the current model by simulating several tapes on top of each other and applying the transport current as done in section 4.3 for an individual tape. It is important to note that the boundary condition for the magnetic vector potential A = A 0 only sets the total current flowing in the stack, not the current flowing in each tape, which is determined by the impedance offered by the different current paths. Figure 12 shows the distribution of current density in a five-tape stack. It is clearly visible that most of the current flows in the top and bottom tapes. In the simulated case, a magnetic vector potential of 5 × 10 −4 Wb m −1 generates a total current of about 930 A, but the currents in the various tapes and the corresponding cyclic losses are very far from equal, as reported in table 2. One can also easily verify that the losses of the individual tapes are higher than those calculated for a single tape by means of equation (12), with I 0 equal to the values in the second column of the table. This is because the tapes carry transport current while they are also subjected to the field generated by the other tapes.
In fact, the tapes are electromagnetically coupled, and they approximately behave in the same way as a larger, unique superconductor. If one uses equation (12) with the 'total' values of I 0 and I c , i.e. I 0 = 930 A and I c = 5× 314.15 A, one obtains an AC loss value of 5.04 × 10 −3 J m −1 , which is not too far from the 5.95 × 10 −3 J m −1 numerically calculated for the stack (see table 2). As it will be discussed later in section 5, this loss increase due to electromagnetic interaction between the tapes in a stack limits the possibility of using stacks for long-distance transmission of AC power.

A look at real applications
The simplicity of the used model notwithstanding, some results of the proposed exercise highlight important features that are at the core of present day's research aiming at the commercialization of superconducting applications. This applies particularly to high-temperature superconductor (HTS) tapes, which are manufactured in two forms: as multi-filamentary tapes embedded in a silver matrix or as coated conductors where the superconductor material is a thin film a few millimeters wide and only one or few micrometers thick. In the first case, the area of the superconducting filaments has an approximate elliptical shape, so that the approach of the proposed exercise can be directly applied. In the second case, the geometry is more 'extreme' and as such more challenging to simulate. But because of the extremely large widthto-thickness ratio, some of our findings (such as the influence of the orientation of the field on the magnetization losses) are even amplified.
For example, superconducting transformers have the tape wound in the form of a solenoid. In most of the solenoid, the field is oriented parallel to the flat face of the tape (i.e. θ = 90 • with the definitions of figure 8), which results in low AC losses. At the two ends of the solenoid, however, there is a substantial perpendicular component of the field (see figure 13), which causes very large AC losses. In [20] the authors report some calculation for a solenoid used as a shielded-core superconducting fault current limiter (SFCL) where the last turn has  Magnetic field lines of the cross section of a cable carrying current. Due to the interaction between the fields generated by the various tapes, the field is mostly oriented in the 'good' direction, i.e. parallel to the flat face of the tapes. The 'bad' perpendicular component is almost completely canceled. The cancelation is more effective as the gap between them is reduced. a power dissipation ten times higher than that of the turns situated in the central part of the solenoid. In general, measures can be taken to avoid excessive power dissipation at the ends, e.g. by using ferromagnetic shields or superconductors made of narrower transposed strands.
Resistive-type SFCL are usually manufactured in the form of bifilar windings [21]. The current of adjacent turns flows in opposite directions, thus strongly reducing the 'bad' field component perpendicular to the flat face of the tape (i.e. θ = 0 • with the definitions of figure 8) and then reducing the losses. Bifilar windings have the additional advantage of reducing the coil's inductance virtually to zero.
For transporting AC power over long distance, one could be tempted to use compact tape stacks, because of the extremely reduced cross section. However, this would not be an efficient solution. From the electromagnetic point of view, the stack behaves like a monolithic conductor, and this has much larger losses than the same number of independent tapes, as shown in table 2. This is why AC power cables are manufactured by putting the tapes side-by-side around a cylindrical former: the 'bad' component of the magnetic field generated by the tapes near the gap between them is strongly reduced (see figure 14), and so are the AC losses. In this way, the AC losses of each tape are even lower than those of as many individual tapes carrying the same current. This was confirmed by measurements in [22].

Conclusion
With this work, we presented a finite-element model, implemented in the open-source program FreeFEM, that can be used to introduce students to calculation of AC losses in superconductors. Estimation and reduction of AC losses are important topics for designing and manufacturing efficient and cost-competitive superconducting applications. With the presented model, students can learn about the influence of the superconductor's geometry and of field orientation on the AC losses as well as on the effects of assembling superconducting tapes in bigger cable structures.
The model is implemented in an open-source program and runs on several platforms as a stand-alone application. The solution of the proposed exercises is very fast (tens of seconds), so that the students can quickly simulate different cases, both in the classroom and on their own.