**1. Introduction**

With the increasing power load and increased short-circuit capacity, the short-circuit current of the grid will continue to rise and will gradually approach the limits of the breaking capacity, which will greatly affect the stability of the power system and equipment. Among the current-limiting equipment, the superconducting fault current limiter (SFCL) is considered as one of the most effective methods to solve the problem of exceeding short-circuit current. Thus, simulating and conducting experimental research on it is significant [1].

SFCLs can be divided into resistive, inductive, three-phase reactor, saturated iron core, and bridge types. The principle of resistive SFCL is the transition of the superconductor from superconducting to normal state. They are relatively simple in structure, but continuous losses during rated operation is unavoidable. The inductive SFCL omit the current leads, but in general use an iron core, which makes this SFCL relatively heavy and costly. Recently, the inductive SFCL has received more attention [2–5].

In 1996, the Swiss ABB technology company developed a 1.2 MVA/10.5 kV three-phase inductive high-temperature superconducting fault current limiter (HTS FCL), installed it in a hydropower station in Lausanne, and carried out a durability test for one year. It was the first superconducting electrical equipment tested under the actual operating conditions of a power station [6]. The secondary winding is a superconducting magnetic shielding ring made of BSCCO-2212, cooled by liquid nitrogen. In the test, the peak value of short-circuit current was limited to a range from 60 kA to 700 A. The Warsaw Institute of Superconductivity in Poland and the Karlsruhe Institute of Technology in Germany developed an air-core inductive FCL. In 2016 and 2018, 15 kV/140 A [7] and 10 kV/600 A [8] single-phase inductive FCLs were developed.

The inductive current limiter simulation is divided into a steady-state simulation for rated operation and a transient-state simulation for the current-limiting process. When the SFCL is in rated operation, the temperature variation of superconducting winding can be neglected, because the HTS tapes do not generated large heat quickly. The steady-state simulation can be based on the numerical model using the H formulation [9,10] or the T-A formulation [11,12]. These numerical models usually use the finite element method (FEM) or the finite difference method to solve Maxwell equations (including Gauss's law for magnetism, Faraday's law of induction, and Ampere's law with Maxwell's addition); the detailed derivation process is shown in the above references. The input of the simulation model is current, and the output is the current density distribution, magnetic field distribution, and AC loss in the SFCL. When the SFCL is in the current-limiting process, the temperature of the superconducting winding rises significantly, affected by the short-circuit current. Since impedance of the superconducting winding is affected by current, magnetic field, and temperature simultaneously, multiple factors are combined, so a transient simulation of the inductive FCL is more difficult.

In [13], resistance connected in parallel with a switch was used to simulate the resistive component of the superconducting winding, and the secondary induced current was calculated by the finite element simulation based on an electromagnetic model. The advantage of the model is that it provided fast computation speed and could simulate the macroscopic working characteristics of the current limiter, including voltage and current curves. However, since the model, lumped parameter model, equated the resistive component of tape to the parallel connection of the resistor and the switch, only the equivalent parameter value can be obtained. It was impossible to simulate microscopic working characteristics, such as the distribution of current and temperature. In [14–16], the magnetic field shielding characteristics of a single inductive current limiter unit in air or iron core were measured at different currents, and then equivalent inductance of the entire current limiter was estimated from the angle of the alternating magnetic flux and the winding hinge. The calculation result was used in the circuit model calculation. Since the model combined the measured values, the calculation results were more accurate. However, the measured value of a single current-limiting unit cannot fully represent the electromagnetic environment of the current limiter. At the same time, since the current-limiter impedance was calculated by the overall equivalent method, it was impossible to simulate the microscopic properties of the tape.

In this paper, a transient numerical model for the process of limiting current in the inductive FCL is proposed. This numerical model including three main models: Circuit model, electromagnetic model, and heat transfer model. Then, an optimization method to accelerate convergence velocity of iteration is proposed. Based on the model, several typical examples are discussed, including the following parameters: The current and voltage, magnetic field and electrodynamic force distribution, and temperature distribution. A SFCL prototype is fabricated, and a short-circuit test of the SFCL is compared with the simulation results to verify the reliability of the simulation.

#### **2. Methodology of Numerical Model**

## *2.1. Model Overview*

When a short-circuit fault occurs in the power grid, the short-circuit current will far exceed the critical current in the current limiter. Meanwhile, part of the current will flow through the metal layer, and a large amount of joule heat will cause a significant temperature rise in the winding. Therefore, to calculate the temperature change inside the winding, the temperature field in a simulation model is needed, which is combined with the electromagnetic model to calculate the current-limiter impedance. In addition, in actual situations, the voltage value of the transformer is known instead of the short-circuit current, so the circuit model is used to calculate the instantaneous short-circuit current based on the power voltage and current-limiter impedance. Since the calculation results of multiple

physical fields involved in the whole model are interdependent (such as current and current-limiter impedance), multilayer iterative calculations are needed to converge the model and simulate the operating characteristics of the actual current limiter [17,18].

The sophisticated commercial FEM software can solve the complex problems quickly and accurately. The ability of interactive post-processing and visualization make the analyzing of the results easier. The FEM software being used in this paper is COMSOL multiphysics 3.0, and both electromagnetic model and heat transfer model are based on the software. The circuit model is calculated by MATLAB.

The overall calculation structure of the model is shown in Figure 1.

**Figure 1.** Schematic diagram of inductive superconducting fault current limiter (SFCL) transient-state simulation model.

## *2.2. Circuit Model*

The SFCL model can be equivalent to a transformer with a secondary short circuit. The circuit model of the SFCL is shown in Figure 2.

**Figure 2.** Schematic diagram of circuit model.

In Figure 2, *U* is the power voltage, *r* is the line resistance, *i*1 and *i*2 are the primary and secondary currents, respectively, *<sup>R</sup>*1(*YBCO*) and *<sup>R</sup>*1(*Shunt*) are the equivalent resistance of the HTS layer and metal layer in primary winding, respectively, *<sup>R</sup>*2(*YBCO*) and *<sup>R</sup>*2(*Shunt*) are the same as above in the secondary winding, *L*1 and *L*2 are self-inductance of the primary and secondary winding, respectively, and *M* is mutual inductance.

In this model, the calculation of resistance is equivalent to the calculation of the resistance voltage drop. The resistance of the primary and secondary windings is in parallel with the resistance of the HTS layer and the metal layer.

For the HTS layer, the resistive voltage drop is related to the temperature, amplitude, and direction of the magnetic field. Because of the nonuniform magnetic field and different positions of the primary and secondary windings, it is necessary to separately calculate the resistive voltage drop of each turn, and then accumulate the entire winding to obtain the total resistive voltage drop, as shown in Equation (1):

$$\Delta U = \sum\_{N} E\_c \Delta l \left( \frac{I\_{YBCO}}{I\_c(B\_rT)} \right)^n \tag{1}$$

where *N* is the total number of turns, Δ*l* is the length of each turn, *Ec* is the critical electrical field, usually taking a constant value 1 × 10−<sup>6</sup> V/m, n value is 29 for the tape used in this study. *IYBCO* is the current of the HTS layer, and *Ic*(*B*, *T*) is the critical current under magnetic field *B* and temperature *T*, as shown in Equation (2) [19]:

$$I\_{\mathfrak{c}}(B, T) = \frac{I\_{\mathfrak{c}}(T)}{\left[1 + \sqrt{\left(kB\_{\parallel}\right)^{2} + B\_{\perp}^{2}}/B\_{\mathfrak{c}}\right]^{\alpha}} \tag{2}$$

In the formula, *B* is the components of the magnetic field parallel to the superconducting surface, *B*⊥ is the components of the magnetic field vertical to the superconducting surface, *k*, α, *Bc* are physical parameter of superconducting tape, which represent the degradation of critical under different background magnetic field, and the value of *k* is between 0 and 1. In this paper, according to the test results of the critical current magnetic field angle-dependency of the YBa2Cu3O7-x (YBCO) tape, *k* = 0.3672, α = 0.6267, and *Bc* = 0.05249 T.

For the metal layer, resistivity is only related to temperature, but since it is composed of various metals, the resistivity of each metal varies with temperature. The resistivity of superconducting tape as a function of temperature is measured when it is above the critical temperature. Below the critical temperature, a linear cubic spline fit is made at 92 K based on the measured data.

In summary, the overall circuit model formula can be obtained, and discrete iterations are used in the actual calculation, as shown in Equation (3):

$$\begin{cases} \begin{aligned} \boldsymbol{L}\boldsymbol{l}\_{(m)} &= \boldsymbol{i}\_{1(m)}\boldsymbol{r} + \sum\_{\boldsymbol{N}\_{1}} E\_{\boldsymbol{c}} \Delta \Big( \frac{\boldsymbol{i}\_{1(\boldsymbol{Y}\boldsymbol{R}\boldsymbol{C}\boldsymbol{O})(\boldsymbol{m})}}{\boldsymbol{I}\_{\boldsymbol{c}}(\boldsymbol{B},\boldsymbol{T})} \Big)^{\boldsymbol{n}} + L\_{1} \frac{\boldsymbol{i}\_{1(m)} - \boldsymbol{i}\_{1(m-1)}}{\Delta t} - M \frac{\boldsymbol{i}\_{2(m)} - \boldsymbol{i}\_{2(m-1)}}{\Delta t} \\ \boldsymbol{0} &= \sum\_{\boldsymbol{N}\_{2}} E\_{\boldsymbol{c}} \Delta \Big( \frac{\boldsymbol{i}\_{2(\boldsymbol{Y}\boldsymbol{R}\boldsymbol{C}\boldsymbol{O})(\boldsymbol{m})}}{\boldsymbol{I}\_{\boldsymbol{c}}(\boldsymbol{B},\boldsymbol{T})} \Big)^{\boldsymbol{n}} + L\_{2} \frac{\boldsymbol{i}\_{2(m)} - \boldsymbol{i}\_{2(m-1)}}{\Delta t} - M \frac{\boldsymbol{i}\_{1(m)} - \boldsymbol{i}\_{1(m-1)}}{\Delta t} \end{aligned} \end{cases} \tag{3}$$

In the formula, *<sup>i</sup>*1(*YBCO*)(*m*) is taken as an example, the number 1 represents the primary side, (*m*) represents the *m*th step iteration, and (*YBCO*) represents the HTS layer, so *<sup>i</sup>*1(*YBCO*)(*m*) represents the current of the HTS layer in the primary winding at the *m*th iteration; *<sup>i</sup>*2(*<sup>m</sup>*−1) represents the total current in the secondary winding at the *(m*−1*)*th iteration; *N*1 and *N*2 is the total number of primary winding and secondary winding, respectively, and Δ*t* is the size of the iterative step (1 × 10−<sup>5</sup> s).

## *2.3. Electromagnetic Model*

The electromagnetic model is based on the COMSOL AC/DC module, and it is based on Maxwell equations, as shown in Equation (4):

$$\begin{cases} \begin{array}{l} \nabla \cdot B = 0 \\ J = \nabla \times H \\ \nabla \times E = -\frac{\partial B}{\partial t} \end{array} \end{cases} \tag{4}$$

where *B* is the magnetic flux density, *J* is the current density, *H* is the magnetic field strength, *E* is the electric field strength, and *t* represent time.

Figure 3 shows a schematic diagram of the electromagnetic model, which is a two-dimensional axisymmetric graph. Figure 3b is a diagram with mesh generation in simulation, the thickness of the superconducting layer is magnified 10 times to accelerate the simulation. About mesh generation in superconducting layer, 20 grids in z direction and 1 grid in r direction are generated.

**Figure 3.** Schematic diagram of electromagnetic model. (**a**) Structure; (**b**) mesh.

In order to simplify the model and speed up the calculation, *B* (Equation (2)) is the average of parallel magnetic field in one turn, and *B*⊥ (Equation (2)) is the average of vertical magnetic field in one turn. Di fferences in current density at di fferent points are ignored for the same turn. The electromagnetic model is used to calculate the *IYBCO* and *Ishunt* when the total current is determined.

The specifications of the SFCL model are shown in Table 1.



#### *2.4. Heat Transfer Model*

The heat transfer model is based on the COMSOL heat transfer in solids module. As shown in Figure 4, a heat flux is applied to the outer boundary to simulate the process of boiling liquid nitrogen heat transfer, so that the simulated liquid nitrogen region can be directly deleted.

**Figure 4.** Schematic diagram of heat transfer model.

The heat transfer model is divided into two parts. The first part is solid heat transfer, including heat transfer inside the superconducting tape, between tapes, and between tape and insulating material. The basic equation is as follows [20]:

$$
\rho \mathbb{C}\_p(T) \frac{\partial T}{\partial t} = \nabla \cdot (\lambda(T) \nabla T) + Q \tag{5}
$$

where ρ is density, *T* is temperature, *Cp*(*T*) is specific heat capacity, λ(*T*) is thermal conductivity, *t* represent time, and *Q* is the heat source in this model, which is the density loss *E*·*J* of the winding.

In the model, the simulation space is divided into three parts: Superconducting, metal, and insulating layers. The metal layer is composed of multiple layers of different materials, and the specific heat capacity *Cp*(*T*) and thermal conductivity λ(*T*) as functions of the temperature of each layer are different. To simplify the model, uniform thermal parameters are needed. *Cp*(*T*) and λ(*T*) of each layer are equivalently calculated, as shown in Equations (6) and (7):

$$\mathbf{C}\_{eq} = \frac{\sum \mathbf{C}\_n(T)\rho\_n t\_n}{\sum\_N \rho\_n t\_n} \tag{6}$$

$$
\lambda\_{eq} = \frac{\sum\_{N} \lambda\_{\text{ll}}(T) t\_{\text{m}}}{\sum\_{N} t\_{\text{m}}} \tag{7}
$$

where, respectively, *Ceq* and <sup>λ</sup>*eq* are equivalent specific heat capacity and thermal conductivity and *Cn* and λ*n* are specific heat capacity and thermal conductivity of each metal layer, and ρ*n* is the density and *tn* is the thickness of each metal layer. *N* represents the types of mental materials.

The second part of the heat transfer model is the heat exchange between the outer boundary of the winding and the liquid nitrogen, which is dominated by thermal convection. However, since the increased temperature of the winding causes the liquid nitrogen to boil, the heat transfer model is complicated. In order to simplify the calculation, a heat flux is applied to the outer boundary to

simulate the process of boiling liquid nitrogen heat transfer, so that the simulated liquid nitrogen region can be directly deleted. The basic equation is as follows:

$$Q = h(T\_s - T\_0) \cdot A \cdot (T\_s - T\_0) \tag{8}$$

where *Ts* is the temperature of the heat flux boundary, *A* is the boundary cross-sectional area, *T*0 is the liquid nitrogen temperature (77 K), and *h (Ts* − *T*0*)* is the effective nonlinear steady-state convection coefficient (W/(K·<sup>m</sup>2)). The value of *h (Ts* − *T0)* is liquid nitrogen boiling heat flux curve in the simulation, which is cited in [21].

#### *2.5. Model Convergence and Optimization*

#### 2.5.1. Middle Voltage Loop Optimization

As shown in Figure 1, the iterative calculation is divided into three layers: Inner, middle, and outer. The middle voltage loop is used to iterate in the circuit model (calculating the primary and secondary currents) and the electromagnetic model (calculating the resistive voltage drop of the superconducting winding) under a specified supply voltage. The intermediate variable is the superconducting winding resistive voltage drop *Usc\_p*, *Usc\_s*, and the updating method is shown as Equation (9):

$$\begin{cases} \mathcal{U}\_{\text{sc\\_p}(k)} = \mathcal{U}\_{\text{sc\\_p}(k-1)} - \Theta \Delta \mathcal{U}\_{\text{sc\\_p}(k)}\\ \mathcal{U}\_{\text{sc\\_s}(k)} = \mathcal{U}\_{\text{sc\\_s}(k-1)} - \Theta \Delta \mathcal{U}\_{\text{sc\\_s}(k)} \end{cases} \tag{9}$$

where *Usc\_p*(*k*) represents the primary resistive voltage drop at the *k*th iteration and *Usc\_s*(*k*) is the secondary voltage drop, where *k* refers to the number of iterations in the middle voltage loop; <sup>Δ</sup>*Usc\_p*(*k*) is the difference between *Usc\_p*(*k*) and *Usc\_p*(*k*−1); and θ is the learning rate. In order to accelerate convergence in the early stage, a method of exponential decay plus constant is used to calculate θ, such as Equation (10):

$$
\theta(k) = \beta \epsilon^{-\frac{k}{\tau}} + \xi \tag{10}
$$

where β is the exponential decay coefficient, and the value is large (0.2), which ensures that the model can be quickly converged at the beginning; τ controls the decay rate; and ξ is the constant learning rate, which ensures that the learning rate of the model does not decay to a too small value after multiple iterations and cannot update. In the simulation, β = 0.2, τ = 30, and ξ = 0.01.

#### 2.5.2. Inner Voltage Loop Optimization

The inner voltage loop is used to calculate the current distribution in the superconducting and metal layers under a specified total winding current. In the actual calculation, due to the characteristics of the HTS, its current will enter a large nonlinear region when it reaches the critical current, as shown in Equation (1). At this time, if only *Usc\_p* and *Usc\_s* are used as the convergence conditions, the model will oscillate easily and cannot converge.

To solve this problem, when the total winding current is close to or greater than the critical current, *Umetal\_p* and *Umetal\_s* are used as convergence conditions. The superconducting layer has nonnegligible resistance at this time, a part of current flow through the metal layer, which results in a nonnegligible voltage drop in the metal layer, and this voltage drop can be directly calculated. Because the resistance of the metal layer is constant at the same temperature, there is no oscillation.

Figure 5 shows the number of iterations per step in two cycles (40 ms) before and after optimization. The ordinate represents the iterations of the inner current loop when the outer voltage loop is iterated once; the abscissa is the iterations of the outer voltage loop. Figure 5 shows that the number of iterations at some special points can be greatly reduced after optimization.

**Figure 5.** Number of iterations at each step before optimization and after optimization.

#### **3. Results and Discussion**

## *3.1. Typical Examples*

#### 3.1.1. Current and Voltage

The primary and secondary current curves at different supply voltages were obtained and used to analyze the macroscopic operating characteristics of the SFCL in the simulation. The operating conditions under two typical voltages (voltage amplitudes of 0.5 V, 30 V) are listed, including primary current, primary superconducting layer current, secondary current, secondary superconducting layer current, and voltage.

Figure 6 shows the simulation results of the voltage and current curves at a voltage magnitude of 0.5 V. The current flows from the superconducting layer before the primary current reaches 60 A for the first time; after the current exceeds 60 A, the shunting phenomenon begins to appear. It shows that the superconducting layer begins to show obvious resistance because its current exceeds the critical current. However, compared with the metal layer resistance, this resistance is still small, so only a few amperes of current flow through the metal layer.

**Figure 6.** Simulation results of voltage and current curves at voltage magnitude of 0.5 V.

When the waveform is stable, the primary current amplitude is 53.9 A, the secondary current amplitude is 51.2 A, and the phase di fference between the two is 174.6◦. In order to explain the shielding e ffect of the secondary winding on the magnetic flux generated by the primary winding, the sum of the primary and secondary currents is defined as "net current," which generates magnetic flux (the ratio of primary to secondary is 1:1). The purple curve in Figure 6 is net current; the amplitude is 5.1 A and bias is 6.0 A, which is 9.5% of the primary current. The secondary winding shields most of the magnetic flux generated by the primary winding. Only 9.5% of the primary current builds an alternating magnetic field, so the overall leakage inductance is close to zero. The primary and secondary currents flow through the HTS layer, so the overall impedance of the SFCL is close to zero. The phase di fference between the primary current and the voltage waveform is 84.6◦. The leakage inductance of the primary and secondary windings makes the SFCL inductive.

Figure 7 shows the simulation results of the voltage and current curves at a voltage magnitude of 30 V. The first peak of the primary current is 652.7 A, and the current wave becomes steady state in the fifth cycle with the peak value of 427.4 A. After the current of the secondary winding becomes stable, the peak value is 127.1 A, and the amplitude of net current generating the alternating magnetic field is 399.5 A, which is 93.5% of the primary current. It indicates that most of the primary current is used to construct the alternating magnetic field, and the secondary current can only shield a fraction of the magnetic flux generated by the primary current. The SFCL has large inductance; meanwhile, the resistive component of the winding is also increased due to the large primary current. The phase di fference between the primary current and the voltage is 57.6◦ in steady state.

**Figure 7.** Simulation results of voltage and current curves at voltage magnitude of 30 V: (**a**) Overall curves, (**b**) net current that generates magnetic flux, (**c**) details of curves at first two cycles.

At this moment, due to the large primary and secondary current, a considerable amount of joule heat is generated, and the influence of the increased winding temperature on the HTS strip cannot be ignored. As shown in Figure 7a, the heat accumulation results in decreased saturation currents

in the primary and secondary superconducting layers, because the critical current decreases as the winding temperature rises, and the equivalent impedance increases. The temperature rise curve and temperature distribution of the winding will be analyzed in detail later.

#### 3.1.2. Magnetic Field and Electrodynamic Force

By analyzing the distribution of the magnetic field, the distribution of critical current can be derived, then the distribution of power density can be inferred. In addition, the maximum electrodynamic force acting on the limiter depends on the first peak of primary current under short-circuit conditions. Figure 8 shows the magnetic field distribution at the first current peak of the primary winding under a voltage magnitude of 30 V.

**Figure 8.** Magnetic field distribution at first current peak of primary winding under voltage magnitude of 30 V: (**a**) Flux density norm value distribution, (**b**) magnetic field direction diagram.

In the area between primary and secondary windings, the entire secondary winding, and the inside and ends of the primary winding, the magnetic field is largest. Among them, the secondary winding and inside and ends of the primary winding are mainly influenced by the magnetic field parallel to the tape. According to Equation (2), this magnetic field has little effect on the critical current, which results in a low power loss density; the ends of the primary winding are mainly affected by the magnetic field vertical to the tape, and this results in low critical current density and high power loss density. To improve the recovery speed after cutting off the short fault, the above problem should be taken into consideration in the design of a cooling structure. If necessary, tapes with higher critical current and better cooling structure can be used in the ends of the primary winding.

Figure 9 shows the electrodynamic force distribution at the first current peak of the primary winding under a voltage magnitude of 30 V. The electrodynamic force affecting the secondary winding is small and the direction is negative; in the primary winding, the inner side is influenced by positive force and the outer side, and the winding is squeezed in the r-axis direction. For electrodynamic force in the z-axis direction, the secondary winding has a small influence and the direction is positive; in the primary winding, the end is greatly affected by a negative force and the middle is influenced by a positive force, and the winding is squeezed in the z-axis direction.

As shown in Figure 9c, different parts of the winding are influenced by electrodynamic force in different directions. Electrodynamic force is a macro performance of the Lorentz force; the magnetic field of the conductor is generated by its own current, so the direction of electrodynamic force does not change with the direction of the current.

**Figure 9.** Electrodynamic force distribution at the first current peak of primary winding under voltage magnitude of 30 V: (**a**) R-axial component, (**b**) z-axial component, (**c**) electrodynamic force direction.

To ensure that the SFCL will not be destroyed in the structure under short-circuit current, the electrodynamic force affecting different parts of the winding should be considered in the design. The superconducting coils in the primary and secondary windings of the upper half of the SFCL are numbered and calculated by the electrodynamic force they affect; the numbers are shown in Figure 9c and the calculation results are shown in Table 2.


**Table 2.** Stress of pancake winding under first peak current of primary winding.

The force acting on the secondary winding is not large, and the external shedding force in the z-axis should be considered in the design; the force acting on the primary winding is large, more attention should be paid to the outward force, and it is necessary to enhance the structure because of the large force.
