Three-dimensional electromagnetic modeling of practical superconductors for power applications

# Three-dimensional electromagnetic modeling of practical superconductors for power applications

M. Kapolka
29th of January 2016
Mid-term PhD report for the program of Physical Engineering in the Faculty of Electrical Engineering in the Slovak Technical University
(FEI STU)
Supervisor:
Enric Pardo
Institute of Electrical Engineering,
Bratislava, Slovakia

### Foreword

This document is the mid-term PhD report from Milan Kapolka. The mid-term examination took place on 25 of February 2016 and the report was submitted by 29th of January 2016, practically as is. Only minor changes have been made before uploading this document, concerning English and expression (only a few corrections). We also removed figure 2.2 in order to avoid copyright issues.

We acknowledge Francesco Grilli for the FEM modeling in figures 4.13 and 4.14.

Enric Pardo,
Bratislava 20th of May 2016

## Chapter 1 Introduction

Thanks to the improvement of industrial production and quality of Hight Temperature Superconductors (HTS) in the last decade, HTS become more important in a wide range of power applications, such as large scale magnets, motors and generators. Superconductors for power applications are produced as a round wires, tapes or bulk samples.

HTS are feasible for hight power applications because of the smaller consumption demand of the cooling system in comparison to the whole power device. Real industral devices contain superconductors with 2D and 3D geometry (coated conductors tapes represent an example of 2D geometry, since the superconducting layer is very thin). In order to know the feasibility and the optimum design of a certain device, there is a need of software tools. These Numerical tools have hight requirements like fast computation, a physical model for any relation of the superconductor and any complex geometry such as coils, motors and generators.

In order to create such tool, the eddy current problem must be solved. There are many methods, mainly based on the Finite Element Method (FEM) and variational principles [1]. Bossavit [2] introduced the variational method for 3D objects. The numerical method is based on a functional that allows smooth relations and include the anisotropic case, where the electric field may not be parallel to the current density when the magnetic field is not perpendicular to (force free situations). However, Bossavit did not solve any example in a 3D geometry. The functional is expressed in terms of the magnetic field, and hence the method needs to solve not only in the superconductor but also a large volume of the surrounding air. A few years later, Prigozhin [3] showed a variational principle only for thin films (2D case), including various samples with perpendicular applied field and holes. The advantage of Prigozhin’s functional is that it is written in terms of , and hence the problem is restricted to the superconductor volume. Recently, Badia [4] provided a deeper insight of the formulation of the functional and proposed the anisotropic power law with specified perpendicular and parallel critical current densities. He solved an approximation for a two-layer counter-wound cable by a stack of the infinite slabs, which is actually a mathematically 2D problem.

In this work we present a 3D variational model based on a new functional that restricts the problem in the superconductor volume. We show the magnetization process of a thin film and a 3D bulk sample. We compare the new model of the thin film geometry with Halse thin film formula[5], reaching a good agreement. We also compare a striated tape, where the filaments are connected by linear material, with a FEM model. We present several results for a thin film with constant critical current density, , magnetic-field dependent , and an anisotropic relation. The last model is the 3D cubic sample.

We find the time dependece of the AC loss dependence for each situation, the current distribution in isotropic and anisotropic cases of the thin film geometry, as well as the bulk sample. In the cubic bulk sample, we found a non-negligible component of the current density in the direction of the applied field.

The next step in our work is to model more complex structures like MgB wires or spiral cables for magnets. We also plan to decrease the computation time by parallel programming.

## Chapter 2 Background

The discovery of superconductivity was 100 years ago, when Kamerlingh Onnes measured zero resistivity in mercury. He obtained this result after liquifying helium and reaching the very low temperature of 3 K in 1908. Superconductivity follows two main hallmarks. The first is perfect conductivity. Persistent current in a superconducting ring without time decay has been measured; Measurements of nuclear resonance showed superconducting current without change during years [6]. The second hallmark of superconductivity is perfect diamagnetism. External applied magnetic field is expelled from the interior of the superconductor, a phenomenon known as Meissner effect. After Kamerlingh Onnes experiment follows years of discovery of a lot of new superconducting materials, first in pure metals and then alloys and composites. Later, Bernortz and Muller discovered compounds with higher critical temperature, around 30 K. These are named as hight temperature superconductors (HTS) and many of them (those used for power and magnets applications) are superconducting above the liquid nitrogen temperature. Nowadays, available commercial superconducting wires or tapes are NbTi, NbSn, BSSCO with two types Bi2212 (BiSrCaCuO) and Bi2223 [(Bi,Pb)SrCaCuO], REBCO (REBaCuO), where RE is Rare Earth usually Y, Gd or Sm and MgB. All of them present high critical current density and high critical magnetic field at temperatures sufficiently below their critical temperature, . Low temperature superconductors are NbTi, NbSn and high temperature superconductors are Bi2212, Bi2223 and REBCO, while MgB is usually classified as ”mid temperature” superconductor.

There are two different kinds of superconductors, type I and type II, which can be explained by the Ginzburg-Landau theory. Type I superconductors have a critical magnetic field . Above this level, the material is no more superconducting. In type I superconductors, and are nonzero only in a layer of thickness , which depends on the superconductor. Type II superconductors present current vortices of quantum origin (although the Ginzburg-Landau theory already predicts their main properties), which explains the penetration of magnetic field inside the sample and the creation of two critical magnetic fields and . Below , the superconductor is diamagnetic; it is perfectly shielded by Meissner current. Between and , vortices penetrate the superconductor and a small dynamic resistance occurs. Above , superconductivity disappears.

The main advantage of technological type II superconductors is their high critical current density , with zero resistivity or small dynamic resistivity in high magnetic fields. Standard metals already reached their limits in power applications. Although large scale magnets with normal metal can reach up to 40 Tesla in s small bore when cooled by water under pressure, (such as in copper-based Bitter magnets) such magnets consume prohibitive amounts of energy by Joule heat. On the other hand, superconductors with high current density induce high magnetic fields with negligible Joule heat. Therefore, superconducting large scale magnets are light and small in comparison to normal metals. Large bores for fields above 2T are only possible with superconductors. The world record of generated DC magnetic field of 40 T has been possible by combining resistive and superconducting magnets. Superconducting magnets are commercially used mainly in medicine (Magnetic Resonance Imaging MRI, material characterization, Nuclear Magnetic Resonance Spectroscopy NMRS), magnetic separation and induction heaters. There is the possibility to use superconductors in the public transportation in high speed levitation trains. Another power application are power transmission cables. Their advantage is mainly in big metropolis where there is a big effort to upscale old underground power lines. However this is not possible with normal conductors, which must increase either the cable cross-section or voltage. Superconductor power cables are relatively small, even with cooling shell, and they enable increasing power density without necessarily increasing voltage. Power cables can be used as well as current fault limiters. One property of superconductors is that current density can be smaller or equal to practically without resistance. When the current density is higher than , the superconductor presents resistance. The resistance heats up the superconductor and it transits to the normal state, where the resistance is even higher, acting as open circuit. An interesting application is superconducting magnetic energy storage (SMES). SMES is a superconducting short circuit winding with induced persistent current. SMES has high power capacity and fast response on electric grid demands. Standard applications as motors, generators and transformers have higher efficiency. Superconducting generators are lighter, smaller, with higher magnetic field and with higher torque. Those features make them very useful in wind turbines. Offshore applications need very light devices with the high power density. Superconductors are needed as well for scientific measurements like particle accelerators, such as that in CERN. A promising energy application are fusion reactors, such as ITER (International Thermonuclear Experimental Reactor). Only superconductors are able to create magnetic fields hight enough to confine plasma in the chamber.

The necessity for cooling systems is the reason why superconductors are not used in more commercial applications. Joule heat, which occurs only at the alternating mode, creates heat loss within the cryostat. A problem is that it takes 10-100 times more power to extract the Joule heat to keep the required low temperatures. Therefore, there are a lot of research projects focused on creating the superconducting wires or tapes with reduced AC loss. Nowadays cryocoolers are commercially available with a closed circuit of helium. Recently, the efficiency of these cryocoolers has been increased and they can cool down to almost any temperature. High power devices are more feasible, because the energy demands of the cooling system is smaller compared to the rated power of whole device.

### 2.1 Low and hight temperature superconductors

The superconducting state is reached under three conditions. Current density must be below the depairing current density , the temperature must be below critical temperature and the magnetic field must be below critical magnetic field . Those conditions are exactly valid for type I superconductors. A different situation is for the type II superconductors. Above the critical current , the material presents macroscopic resistance due to flux flow, although superconducting vortices are still present. In addition, there are two critical magnetic fields and , the superconducting state is often extended up to higher magnetic fields than for type I. For , type II superconductors behave mostly the same as type I, while for , type II superconductors present vortices. Type II superconductors can be used in the wide range of the applied magnetic fields for power applications such as large scale magnets and rotating machines.

Low Temperature superconductors (LTS) present up to 20 K and up to 25 T. These values are upper limits since they depend on the other factors. High temperature superconductors (HTS) such as REBCO, present above A/m at liquid He, which is higher than that for LTS, the upper limit of is 134 K (for HgBaCaCuO), and up to 120 T (which is an extrapolated value). HTS opened new possibilities for wider commercial usage because of cheaper cooling systems, although the tape or wire cost are still an issue.

It took long time to develop and optimize LTS materials suitable for magnet applications. Hulm and Blaugher studied a wide range of the alloy superconductors VNb, NbTa, VCr, VMo,NbCr, NbMo, NbW, NbHf, NbZr and NbTi in 1961 [7]. Time showed that the best superconducting alloy was NbTi with =9 K embedded at Cu or Ni stabilization matrix. Multifilamentary wires were prepared by powder in tube (PIT) method, tested and optimized for various applications. Wire specification exceeded = A/m at 5 T and 4.2 K when ferromagnetic pinning centres were inserted. Superconducting magnets of NbTi can provide magnetic fields only up to 12 T. Another low temperature superconductor is NbSn with critical magnetic field up to 23 T. This alloy was discovered in 1961 with 18 K. Wires are prepared by the internal tin process (IT) which reach = A/m at 12 T and 4.2 K. Alternative powder in tube process used Sn rich compounds (NbSn or NbSn). Addition of Ti increased up to around 26 T.

The first hight temperature superconductor wire with above the liquid nitrogen temperature (77 K) is YBaCuO (YBCO), which was discovered in 1987 [8]. It was prepared by the oxide powder in tube (OPIT) method. Another HTS prepared by a similar process is BiSrCaCuO (Bi2212) embedded in a Ag tube found in 1989 with = 85 K and up to 45 T [9]. Few years later, (Bi,Pb)SrCaCuO (or Bi2223) with = 110 K [10] was discovered. Bi2223 multifilamentary composite wires, which are named as the first generation (or 1G) wires, present 200 A at 77 K in self field. Bi2223 wires present higher critical current densities at higher temperatures compared to Bi2212. On the other hand, Bi2212 is a useful superconducting material up to 20 K, specially at 4.2 K. Although the development of this type was slower than Bi2223, state of the art wires show high performance at 4.2 K and high magnetic fields ( A/m at 26 T) [11]. Bi2212 wires are promising for high field magnets, enabling magnetic fields beyond those created by NbSn and showing higher engineering current density than REBCO at 4.2 K. It is expected to use Bi2212 wires in the range of 30-50 T. Only Bi2212 is prepared as round wire, which can be very useful for magnets with complex 3D winding geometries. The highest current density reached is = A/m, at 4.2 K and 10 T [12].

Another evolution step in HTS are coated conductors, named as Second generation (or 2G) tapes. These are using REBCO material, where means RE is a rare-earth (usually Y,Gd or Sm). The general structure of the coated conductor is following from bottom to top: metallic substrate, multifunctional oxide barrier, buffer layer, HTS layers and silver stabilization layer. The superconducting layer is grown on the textured buffer layer. All layers are deposited epitaxially on the textured bi-axiall substrate. Recently, the REBCO structure has been simplified, which enables an easier production and an increase in the production lenght. The production rate is several tens of meters per hour. Fujikura [13] prepared a wire of 1 cm width with 572 A at self-field and 816.4 m[7]. Deposition can by done by a several ways, such as Pulse Laser Deposition (PLD), Metal-Organic Chemical Vapour Deposition (MOCVD), Electron Beam Deposition (EBD), and Metal Organic Deposition (MOD). Their price is high because of the complex production process. The commercial price of 2G wires could be below \$ 50 /kAm at 77 K and 0 T[7].

Overpresure during the production process of Bi2212 tape increased the engineering current density and it reached the highest among all superconducors above 17 T. However, the highest upper critical field is still for REBCO tape.

The last commercial superconducting material is MgB. Superconductivity was discovered in 2001[14] with 39 K and 10 A/m at 4.2 K [15]. MgB presents moderate anisotropy in the perpendicular field compared to the parallel field. Additives based on carbon decreased the anisotropy and enhanced . In wires and tapes, suitable barrier materials between MgB and the stabilizing layer (usually a Cu-based alloy) are Nb, Ti and Ta, because of small chemical reaction between them and MgB. There are three ways of fabrication of MgB. For the powder in tube method, the wire is prepared in two ways, either ex situ or in situ. The ex situ route uses already reacted powder, which is later inserted in a metal tube. To this follows mechanical deformation, which creates mono-core wires and 10 A/m at 13 T and 4.2 K [16]. On the other hand, the in situ process is different only by using non-reacted powder mixture of Mg and B. Comparing to the ex situ method, is the highest, A/m at 13.2 T and 4.2 K[17]. The last way of preparing MgB is the internal Mg diffusion (IMD) process. This method uses an Mg rod and B powder. The powder is between the rod and the outer metal tube. The highest current densities achieved are A/m at 3 T/20 K and A/m at 10 T/4.2 K. The higher is because of higher connectivity of MgB. MgB is a very promising material for power applications with low applied magnetic fields. This superconducting material is cheaper than any other HTS and it is easy to fabricate. In addition, MgB can be prepared with twisted filaments, which highly reduces the AC loss.

### 2.2 Vortices and vortex pinning

Abrikosov predicted the existence of superconducting vortices by means of the Ginzburg-Landau theory. The vortex has roughly tube structure with a normal core zone at the center and superconducting current flowing around the tube. The vortex line is parallel to the induced magnetic field lines and carries the value of one fluxon of magnetic flux, where is Planck’s constant, is the speed of the light and is the electron charge. For applied magnetic field , the magnetic field partially penetrates into the material through vortices, while type II superconductors are still superconducting. A net induced superconducting current creates a driving force on the vortex . This driving force causes the vortices to start to move, which induces electric field and causes dissipation in the superconductor with power density [18]. Movement of vortices can be avoided by pinning centers, which create a force in opposite direction to the driving force. Pinning centers are regions of non-superconducting material, such as voids, nano-particles, precipitates or in-homogeneities; or defects in the crystal lattice, such as dislocations or twin planes, where superconductivity is partially suppressed. When the driving force is smaller than the pinning force, vortices do not move; therefore there is no dissipation.

### 2.3 E(J) relations

The relation of type II superconductors can be deduced by microscopic theories such as collective thermal flux creep or vortex glass. The depends on the activation energy as [19], [20]

 E(J)=Ece−u(J)kT (2.1)

and

 (2.2)

When , as often seen in experiments, and becomes

 E(J)=Ec(|J|Jc)nJ|J|. (2.3)

is the critical electric field for and is the critical current density, the Factor is and is Boltzmann’s constant. Equation (2.3) is a power law for a non-linear conductor. This smooth relation also allows to model other dependences such as , and . The limit cases fot the factor are , which is the Critical State Model (CSM) and , which is the linear Ohm’s law. The smooth power law is more realistic than CSM. Equation (2.3) is an isotropic relation, where electric field and current density are parallel.

Recently published anisotropic relation [21] with the force-free regions is

 E(J)=Ec⎡⎣J2∥J2c∥+J2⊥J2c⊥⎤⎦n−12⋅(Jc⊥Jc∥J∥Jc∥e∥+J⊥Jc⊥e⊥), (2.4)

where can be non-parallel to , is the critical current density when is parallel to and is the critical current density when is perpendicular to , is the parallel current component and is the perpendicular component. Where is

 J||=B⋅J|B| (2.5)

and is

 |J⊥|=|J×B||B|. (2.6)

### 2.4 Critical state model (CSM)

The Critical state model proposed by Bean [22] is a macroscopic theory that describes the response of superconductors under an applied magnetic field , a transport current or a combination of both, although in this work we concentrate on applied magnetic fields only. The most important cases for this work are infinite slabs and thin films. The original CSM statement by Bean is:“Any electromotive force induces a current with constant critical current density ”. That is the same as the multi-valued relation in figure (2.2).

In this work, we consider only an alternating applied magnetic field. If the applied field is smaller than the penetration field , then the superconductor is completely shielded by Meissner currents. There is no superconducting current flowing inside the sample except in . The magnetic field created by the current is on the surface of the superconductor with opposite sign to the applied field. Further increase of the applied field above creates superconducting vortices, which creep into the superconductor. This causes that current starts flowing deeper than inside the material. For infinitely long samples, current flux lines are parallel to the edges of the sample. Thanks to the infinitely long geometry, the magnetic field that creates the region with creates a uniform magnetic field at the rest of the sample. As a consequence of the critical state model, the current density is always equal to the critical current density for slabs and infinite cylinders. For thin films, there are regions with . In principle, this also applies for general 3D shapes.

Assuming negligible , the induced current in a infinite slab or a cylinder starts flowing from the surface and with increasing field is penetrating more towards the sample centre (figure 2.3a). Current has positive sign at the left side and negative on the right side [22] [23]. The current penetration depth increases with increasing the applied magnetic field. For a large enough magnetic field, the superconductor is fully saturated by the current density. The main difference between the slab and the thin film is the non-zero current density in the space with zero magnetic field inside the sample (figure 2.3 b). After ramping up the magnetic field and ramping down to zero, the current density changes the sign at the surface and starts to penetrating again into the sample (figure 2.3 e). At zero applied field, the magnetic flux is trapped inside the superconductor with the same sign as the applied field at the ramping up. Further ramp down causes penetration of negative flux (figure 2.3 g). The positive trapped flux is completely rewritten with negative flux at the minus peak of applied field [24][25].

The basic equations of the critical state model for the two geometries of thin film and slab (figure 2.4) are the following. The current density for a slab [24] is

 Jy(x) = Jc,−W

where is the slab width and . Its AC loss per cycle and unit length (in the direction) is [23]

 Q = Dw4μ0H3m3Hp,Hm>Hp.

where D is the height of the sample, is the applied magnetic field amplitude and is the penetration field, . The current density for the thin film is[26, 5, 24]

 Jy(x) = 2Jcπarctancy√(b2−y2),|y|

where

 b=wcoshHaHc, (2.10)
 c=tanhHaHc, (2.11)

and

 Hc=Jcdπ. (2.12)

The AC loss for the thin film is [5]

 Q=8μ0J2cw2π[lncosh(πHmJc)−πHm2Jctanh(πHmJc)]. (2.13)

This is a basic example of the critical state model. For a rotating magnetic field, crossed magnetic fields and combination of applied field and transport current, we can get more complex situations, where it is difficult to find an analytical solution. General 3D shapes and 2D plates also require numerical solving. The Bean critical state model above assumes that does not change with applied magnetic field. The Kim model assumes the same relation as the original CSM except a magnetic field dependence of [27]

 Jc(B)=Jco(1+BB0)m, (2.14)

where is the magnitude of the magnetic field.

### 2.5 Eddy current problem

The electromagnetic quantities in superconductors are usually obtained by solving the eddy current problem. The eddy current problem is the task to find the unknown variables of the electric field , current density , magnetic flux density , and magnetic field , which solution must agree with Maxwell’s equations and constitutive equations, as follows

 D=ϵE (2.15)
 B=μH (2.16)
 E=ρJ, (2.17)

where is electric flux density, is the permittivity, is the permeability and is the resistivity. The properties of the material are , , , which are generally non-linear tensors of , and (or , and ), respectively. In standard conductors, the resistivity or conductivity is linear and the permeability of usual magnetic materials is highly non-linear. Superconductors have exactly the opposite situation. Superconductivity is non-linear and, the constant permeability of vacuum is assumed () for the superconductor. The eddy current problem can be solved by any non-linear relation, such as the power law. The critical state model corresponds to the limit of an n-factor approaching to infinite (n-factor above 100 are usually enough to reproduce the CSM).

The case of simultaneous alternating applied magnetic field and alternating current (AC-AC case) is a complicated situation, where the standard eddy current problem cannot be easily solved in superconductors. There are several ways to solve the eddy current problem. There exist three formulations for the Finite Element Method (FEM): -- (where and are the vector and scalar potential respectively), (where and are the current vector and scalar potentials, respectively) and the formulation [28].

The following three formulations are made to solve the differential equations by FEM. In section (2.5.4), we present an alternative way to solve the Eddy current problem by a variational principle.

#### 2.5.1 A-ϕ-J formulation

This formulation is solving the eddy current problem by taking the current density , vector potential and scalar potential as unknown variables. The following equation is general for the formulation of non-linear materials

 ρ(J)J=−˙A−∇ϕ. (2.18)

This relation can be found by Maxwell equations. We start with the magnetic field expressed by the vector potential

 B=∇×A, (2.19)

which is a consequence of . The electric field is then specified according to Faraday’s law

 ∇×E=−˙B. (2.20)

After substituting (2.19) into (2.20) we get

 ∇×E=−(∇×˙A) (2.21)

and then

 E=−˙A−∇ϕ, (2.22)

since for any scalar function . Because of the smooth non-linear relation, we can write it as well as equation (2.18). This is the first equation which must be solved. Now, we start from Ampere’s law (2.23) and Coulomb’s gauge (2.24),

 ∇×H=J+˙D (2.23)
 ∇⋅A=0 (2.24)

In quasimagneticstatics, we neglect the second term of Ampere’s law,

 ∇×H=J, (2.25)

this correspond to neglecting the electromagnetic radiation. The magnetic field follows

 H=μ−1(B)B. (2.26)

Now, we express the current density from (2.25) and (2.26).

 J=∇×[μ−1(B)B], (2.27)

where is a tensor for anisotropic magnetic materials. After substituting (2.19), we have the second equation of the formulation. The eddy current problem reduces to solve equations (2.18) and (2.28).

 J=∇×((∇×A)μ−1(∇×A)) (2.28)

The current is specified by the following constrain

 ∫DiJ⋅dS=Ii, (2.29)

where is a sub-domain cross-section, is the net current and is the cross-section surface differential.

In order to solve the differential equations (2.18) and (2.28), we must exactly set the boundary conditions and constrains and find the solution of the vector potential and electric potential. The air is assumed as non-conducting space. This set of the equations can be written by a matrix and solve approximately with the finite element method. Finally, we have the formulation of unknown variables in space and time .

 M˙X=f(t,X) (2.30)

#### 2.5.2 T-ψ formulation

Another type of extended eddy current problem for superconductors is the - formulation. It is based on the current vector potential and current scalar potential . The initial definitions for both unknown variables are

 J=∇×T (2.31)
 H=T−∇ψ (2.32)

Equation (2.31) comes from and equation (2.32) is consequence of (2.25) and (2.31). When we use Faraday’s law (2.20) and the material constitutive equations (2.16) and (2.17), and

 ∇⋅B=0, (2.33)

we get the following equations,

 ∇×(ρ(J)J)=−∂t[μ(H)⋅H] (2.34)
 ∇⋅[μ(H)H]=0. (2.35)

When we substitute the initial definitions (2.32),(2.31) into them, we obtain

 ∇×[ρ(∇×T)(∇×T)]=−∂t[μ(T−∇ϕ)⋅(T−∇ϕ)] (2.36)

and

 ∇⋅[μ(T−∇ϕ)(T−∇ϕ)]=0. (2.37)

These two equations must be set again with proper boundaries of the domain and resistivity of the air. The net current is

 ∫∂ΩT⋅dl=I. (2.38)

Then, we can combine equations (2.36) and (2.37) again into a matrix, as we did previously in (2.30), and solve the eddy current problem by the finite element method.

#### 2.5.3 H formulation

The last formulation for the finite element method is the formulation. This formulation uses directly the magnetic field. Recently, it became a widely used approach to solve the superconducting problem, since Comsol [29] software has the differential algebraic equation solver. It is possible to enter the differential equation as input to the solver. The formulation is not using any specific gauge of the vector potential and scalar potential. The space discretization is based on edge elements. The formulation [4] starts with the non-linear Ohm’s law,

 E=ρ(J)J+Ec, (2.39)

where is a critical constant voltage in the sub-domains. After substitution of Ampere’s law

 E=ρ(∇×H)∇×H+Ec. (2.40)

Another substitution of Faraday’s law and (2.26) makes

 −∂t[μ(H)H]=∇×[ρ(∇×H)⋅∇×H+Ec], (2.41)

which is the final equation for the formulation.

The condition for net current in sub-domains is based on Ampere’s law,

 ∫∂ΩH⋅dl=I. (2.42)

#### 2.5.4 Variational principles

The eddy current problem can be solved by a variational methods, which is a different numerical method from FEM. The variational method is based on a functional, which is solved by minimization. We call the numerical method in this work as Minimum Electro Magnetic Entropy Production (MEMEP) since the functional minimizes the entropy production [30] . MEMEP uses the formulation based on the current density and the scalar potential for fully 3D objects. The formulation presents the advantage compared to the other functionals in the formulation [2, 4] that it reduces the unknown variables to the sample volume, avoiding to solve unknown variables in the air. Different for integral methods, MEMEP avoids time-consuming matrix inversions.

The functional is derived from the vector potential and and Coulumb’s gauge:

 E(J)+˙A+∇ϕ=0 (2.43)
 ∇⋅A=0 (2.44)

Next, we use the Ampere’s law with neglecting the displacement current , and relations and . Then, Ampere’s law is

 ∇×∇×Aμ0=J. (2.45)

According the vector formula for rotor of rotor

 ∇×(∇×A)=∇(∇⋅A)−∇2A, (2.46)

we get, by substituting Coulomb’s gauge,

 −∇2Aμ0=J. (2.47)

The solution of this differential equation is

 A(r)=μ04π∫rd3r′J(r′)|r−r′|. (2.48)

Then, the derived equations (2.43) and (2.44) are

 E(−∇2Aμ0)+˙A+∇ϕ=0 (2.49)
 ∇⋅A=0 (2.50)

Since we use the Coulomb’s gauge, the scalar potential is the electrostatic potential. Taking a time discretization such that and is a certain initial time, and , where . Then, equations (2.49) and (2.50) become

 E(−∇2(A0+ΔA)μ0)+ΔAΔt+∇ϕ=0 (2.51)
 ∇⋅(A0+ΔA)=0. (2.52)

Now, we start with a functional and we will find that minimizing this functional is the same as solving equations (2.43) and (2.44). The functional is

 L=∫Vd3rf=∫Vd3r[12ΔAΔt⋅ΔJ+U(J0+ΔJ)+∇ϕ⋅(J0+ΔJ)], (2.53)

where is at time , and the dissipation factor is

 U(J)=∫J0E(J′)⋅dJ′. (2.54)

The Euler-Lagrange equations are for the whole functional. Then, we can write the Euler-Lagrangian equation of the functional for (page 192 of [31])

 ∂f∂ΔAi−∑j∂j[∂f∂(∂jΔAi)]+∑jk∂j∂k[∂f∂(∂j∂kΔAi)]=0, (2.55)

and the Euler-Lagrangian equation for is

 ∂f∂ϕ−∑j∂j[∂f∂(∂jϕ)]=0. (2.56)

The Euler-Lagrange equation for results in

 ∇2[ΔAΔt+E(J0+ΔJ)+∇ϕ)]=0. (2.57)

The Euler-Lagrange equation for yields

 ∇⋅(J0+ΔJ)=0. (2.58)

Applying (2.47) to the equation above, then it becomes

 ∇2[∇⋅(A0+ΔA)]=0 (2.59)

Since equation (2.57) and (2.59) apply for the whole space, they imply, (2.51) and (2.52), respectively. The Euler-Lagrange equations show that and are extremes of the functional (2.53). Since the vector potential is not restricted to the sample volume, we rather use and variables with condition , which sets the boundaries directly on the sample surface. For this formulation, is satisfied. The unknown variables are then, and instead of and with condition (2.47). The vector potential could come from external sources like a coil or from the internal currents of the sample. Then, we can write and the functional becomes

 L=∫Vd3r[12ΔAJΔt⋅ΔJ+ΔAaΔt⋅ΔJ+U(J0+ΔJ)+∇ϕ⋅(J0+ΔJ)], (2.60)

where we removed all terms independet on . The boundary condition is set at the surface sample that , where is the normal vector to the sample surface.

## Chapter 3 Model and Numerical method

### 3.1 Mathematical model

The model is based on the variational principle of the Minimum Electro Magnetic Entropy Production (MEMEP). The method is solved by minimizing the functional (2.60) using an original numerical method. The model presented in this work assumes the magnetization case without a transport current. We derived a new formulation of the functional, which contains the vector based on the functional (2.60) for and .

We define the last term of the functional (or Lagrangian) as

 Lϕ=∫Vd3r∇ϕ⋅J. (3.1)

According the following vector formula

 ∇⋅(ϕJ)=J⋅∇ϕ+ϕ∇⋅J, (3.2)

the last term of the Lagrangian becomes

 Lϕ=∮SdSϕJ−∫Vd3rϕ∇⋅J, (3.3)

where at the last equation we applied the Divergence theorem.

 ∫Vd3r∇⋅A=∮SdS⋅A. (3.4)

For the magnetization case, the current density is induced inside the sample only and at infinite . Then, the first component of the Lagrangian is equal to zero and

 Lϕ=−∫Vd3rϕ∇⋅J. (3.5)

In this Lagrangian, the condition is not explicitly satisfied. We can use the vector to describe as . Then , and hence is always satisfied. The vector is substituted in the Lagrangian

 Lϕ=−∫Vd3rϕ∇⋅(∇×T)=0 (3.6)

The last component of the total Lagrangian in (2.60) drops out for the same reasons. Then, we have a new Lagrangian with unknown variable .

 LT=∫Vd3r[12ΔAJΔt⋅(∇×ΔT)+ΔAaΔt⋅(∇×ΔT)+U(∇×T)] (3.7)

### 3.2 Discretization

The entire superconducting sample is divided into small cells with shape of a rectangular prism (figure 3.2). The mesh is uniform in the whole volume. The Cartesian coordinate system is appropriate to model the thin film geometry and bulk samples. Each cell contains 3 independent variables of the current density vector . We assume that the current component normal to each side is constant in the whole surface of the cell. With this dicretization, at any point can be obtained by linear interpolation of and at their respective surfaces. The components are set as a constant value along the entire length of the edges at the appropriate cell edge (figure 3.4a). Again, at any point can be obtained by linear interpolation of at their respective edges. The relation between and is in the coordinate system (figure 3.4 b). In integral form this is

 Jx=1S∫SxdS⋅(∇×T). (3.8)

Following Stoke’s theorem and then integrating along the length around the surface, at a surface with constant is

 Jx,ijk=1Sx,ijk[lz(Tz,i,j+1,k−Tz,i,j,k)+ly(Ty,i,j,k−Ty,i,j,k+1)], (3.9)

where is surface area, and are the edge lengths of the surface, and and are components of vector at the edges, respectively. The components and can be similarly obtained from at the edges. The indexes are given according coordinate system in figure (3.1).

Next, we transform the Lagrangian (3.7) from integral form to sum form. The discretized Lagrangian as a function of crossing all surfaces is

 LT=12Δtns∑i3∑sΔAJisVis⋅ΔJsi+ns∑i3∑s1ΔtΔAaisVis⋅ΔJsi+L3 (3.10)

where

 L3 = 13nx∑iVxiU(Jxiex+Jyi(ri)ey+Jzi(ri)ez) + 13ny∑iVyiU(Jxi(ri)ex+Jyiey+Jzi(ri)ez) + 13nz∑iVziU(Jxi(ri)ex+Jyi(ri)ey+Jziez),

where are the unit vectors in the directions, respectively. For a power law realtion, can be calculated for (2.3) and (2.54) as

 U(J)=EcJcn(|J|Jc)n+1. (3.12)

For the anisotropic relation of (2.4), the relation is [4]

 U(J)=EcJc⊥n+1⎡⎣(J||Jc||)2+(J⊥Jc⊥)2⎤⎦n+12. (3.13)

For our discretization, the total vector potential generated by is calculated as the sum of the contribution of the current density at all surfaces

 AJsij=ns∑jJjVsjasij, (3.14)

where labels the surface where the vector potential is evaluated and is a surface that generates vector potential, is the element of the interaction matrix between all surfaces of kind

 asij=μ04π1VsiVsj∫d3r∫d3r′hsi(r)hsj(r′)|r−r′|, (3.15)

where and are the volumes of influence defined as , and and are interpolation functions, which reach value at surface and at the other surfaces and decay linearly, between surface and neighboring ones of the same kind.

### 3.3 Minimization

The minimization is a numerical method, which solve the functional by finding its minimum. After discretizing the functional (3.10) and setting variables at the sample geometry, we can solve it. The functional is minimized by the following method that we developed. The minimization requires initial conditions, according which the algorithm is searching for a minimum of the functional. These initial conditions must be close to the minimum, in order to reduce the calculation time. The present case is looking for non-linear the eddy current in the superconductor. We set all variables to zero at and the vector potential from external source is calculated by for each surface according the applied magnetic field, where the applied magnetic field follows the direction. The vector contains 3 components , but the thin film case with 1 layer of the cells in the z direction uses only the component, because we assume outside the sample and the boundary condition , where is the unit vector perpendicular to the sample surface. Therefore, the minimization takes approximately 3 times more variables in a bulk sample than in a thin film with the same number of cells.

Now, we describe the case of a thin film, where only the component is present. For change of by at one edge produces a change in and ( and , respectively) in 4 neighbouring surfaces (figure 3.4 b). For this change of , we calculate the change in the Lagrangian. The algorithm is checking all edges and as well . For the 3D case, where all 3 components of are present, we also check and . Then, we find the edge where changing by or reduces the functional the most. The selected edge is later updated by or , respectively. The sequence is repeated until changing in any edge increases the Lagrangian instead of decreasing it. For uniform mesh, magnitude of associated to is

 |δJs|=|δTr|1lt, (3.16)

where with and is the cell dimension in the t direction. Since the magnitude of the current density is around in the solution, we can start with directly and then . Afterwards, in order to increase the accuracy of the solution, we decrease and repeat the whole algorithm until the solution satisfies the input tolerance of 0.1-0.001 of current.

## Chapter 4 Results and discussion

We modeled a superconducting thin film with square shape. We assumed several relations; isotropic relations with constant and , an anisotropic relation, and coupling effects in striated tapes. Then, a bulk sample with cube shape was analyzed. Cube dimensions are mm and an uniform mesh of cells for an isotropic relation with constant .

### 4.1 Comparison to the long thin film formula with constant Jc

In this work, we use the power law relation for the superconductor (2.3) with finite N-factor, which correspond to the slope of curve in log-log scale for superconductors. The next calculations are testing the N-factor dependence. We modeled the cases of long samples with width to length ratio 1:2 (dimensions mm) and ratio 1:4 (dimensions mm) with critical current density A/m, applied field 20 mT, and several n-factors (30,200,1000). Figure (4.1 a) and (4.1 b) show the component of the current density at the mm and mm planes, respectively at the peak of applied field. The result (figure 4.1 a,b) shows a clear dependence between increasing n-factor from 30 to 1000 and increasing sharpness of curve. The current density close to the edges starts to decrease to the magnitude as it is according the Bean model. The higher than in the sample with finite n is because of the smooth curve. At the edges , and hence .

Next, the MEMEP model is verified directly by the formula for current density at a infinite thin film derived from the critical state model [5, 26], equation (2.4). The result agrees very well with the current profile for infinite thin film for n=1000 (figure 4.1 c). However, the current profile at the center of the sample is not completely identical. This can be explained by a not sufficiently long sample, since the longer the sample, the closer the result to the infinite thin film formula.

### 4.2 Finite superconducting thin film with constant Jc

In this section, we present results for the dimensions mm. The thin film contains cells with uniform mesh (figure 4.2). The applied field is parallel with the axes with amplitude 50 mT and sinusoidal shape of frequency 50 Hz. The critical current density is 1 A/m, which corresponds to 360 A of the 12 mm wide tape with thickness of 3 m. For the power law, we use V/m and 30. The tolerance for the critical current density is . We calculated the penetration of current density, current flux lines and AC losses. To calculate the AC loss, several time steps are needed. We modelled 1.25 AC cycles with 40 time steps per cycle. One cycle period is 20 ms, corresponding for 50 Hz frequency (figure 4.3).

Since we start the computations at the zero applied field, the current density with magnitude around progressively penetrates into the sample, until the applied field reaches the first peak at (figure 4.4). The sample is almost completely saturated with critical current density except of the centre of the sample, where . The rest of the sample contains smaller current density, as predicted by the critical state model results for a thin film. The decrease of applied field that follows causes a swap of the current density at the edges and current with opposite sign penetrates into the centre. At zero applied field, around half of the previous penetrated current density is re-written by negative current. When the applied field is at the negative peak, the current is penetrated into the same depth as for the positive peak, but the current flows in the opposite direction. After this, the scenario is repeated for each half cycle practically without any change. The important regions are only the first peak and first valley of the applied field.

The AC power loss per cycle is calculated between the first positive and negative peak of the applied field and in (figure 4.5). The total calculated loss is 1.03 mJ. The computation time of the model was 51 hours for 50 time steps using a model computer with CPU Intel core i7-4771@3.5Ghz, 8GB RAM. The Software is written in C++ and it is using only one core of the CPU. The program does not assume any symmetry.

### 4.3 Finite superconducting thin film with Jc(B) dependence

In this section, we model a superconductor with dependence according to the Kim model (equation 2.14). We used the same size of the sample mm with uniform mesh cells (figure 4.2). The applied field is in the Z direction and presents a sinusoidal shape with amplitude 50 mT and frequency 50 Hz. The critical current density at zero magnetic field is higher than the case with constant , A/m, factor for Power law is 30, the factor for Kim model is 0.5. The tolerance for the current density is 0.1 of . We chose the magnitude of such that the self-field critical current of an infinitely long tape is the same as for the case with constant . The calculations where done by the cross-sectional model in [30].

The current penetration presents a different behaviour than that with constant (figure 4.6). The increase of the applied field creates penetration of the current to the sample in roughly the same shape as for constant . The difference is in the value of the current density, since the magnetic field decreases the critical current density. A low applied field at around 20 mT (figure 4.6 b) makes the same penetration depth of the current, but the current density has a value of approximately 0.8. At the applied field of 35.4 mT, the sample is almost completely penetrated by current of 0.7 (figure 4.6 c). At the moving front line, the current density is higher with value 0.9. After further increasing the applied field to the first peak, current with penetrates the entire sample and decreased to around 0.5 (figure 4.6 d). When decreasing the applied field to zero, the current swaps the sign almost in the entire sample (figure 4.6 e). The current is higher closer to the edge with magnitude 0.9-1 and smaller at the centre, with modulus 0.7. The higher average compared to the peak of the applied field is because the applied field is zero and is not decreased.

At the negative peak of the applied field, the current swapped the sign in the whole sample compared to the positive peak and decreased in the whole sample to the same level of around 0.7 (figure 4.6 f). The magnitude of the current decreased because the increase of the applied field decreases . For the results in figure (4.6), we took 50 time steps in total, which required 120 hours of computing time.

The power AC loss decreases with increasing when the sample is not penetrated. This causes a local minimum at in the power loss curve (figure 4.7). This is because at this time mT, and hence is maximum. The total loss per cycle is 0.9 mJ.

### 4.4 Finite superconducting thin film with anisotropic E(J) relation

The next step is to model the anisotropic case with the relation of equation (2.4). In this case, the electric field and the current density are not parallel. In the computations, we used the same uniform grid of cells. The applied field is with magnitude 70.7 mT () and 287.9 mT () in order to keep the same perpendicular component of the applied field of 50 mT as for the model with constant . The applied field is a sinusoidal signal of 50 Hz. We modelled two cases with 45°(figure 4.8) and 80°(figure 4.8). The parallel critical current density is A/m and the perpendicular current density is A/m, which is three times lower than the parallel current density and corresponds to the value of the previous case with constant . We kept the same dimension of the sample