# Superconducting magnetic bearings simulation using an H -formulation finite element model

Loic Queval, Kun Liu, Wenjiao Yang, Victor M.R. Zermeno, Guangtong Ma

► To cite this version:

Loic Queval, Kun Liu, Wenjiao Yang, Victor M.R. Zermeno, Guangtong Ma. Superconducting mag- netic bearings simulation using an H -formulation finite element model. Superconductor Science and Technology, 2018, 31 (8), pp.084001. 10.1088/1361-6668/aac55d . hal-01832923

HAL Id: hal-01832923

https://hal.sorbonne-universite.fr/hal-01832923v1

Submitted on 9 Jul 2018

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
PAPER

Superconducting magnetic bearings simulation
using an H-formulation finite element model

To cite this article: Loïc Quéval et al 2018 *Supercond. Sci. Technol.* **31** 084001

View the [article online](#) for updates and enhancements.

Related content

- [Experiment and simulation of superconducting magnetic levitation with REBCO coated conductor stacks](#)
Kun Liu, Wenjiao Yang, Guangtong Ma et al.

- [Superconducting magnetic bearings with bulks and 2G HTS stacks: comparison between simulations using H and A-V formulations with measurements](#)
F Sass, D H N Dias, G G Sotelo et al.

- [H-formulation for simulating levitation forces acting on HTS bulks and stacks of 2G coated conductors](#)
F Sass, G G Sotelo, R de Andrade Junior et al.

<figure><img src="image_2.png" /></figure>

IOP ebooks™

Bringing you innovative digital publishing with leading voices
to create your essential collection of books in STEM research.

Start exploring the collection - download the first chapter of
every title for free.
# Superconducting magnetic bearings simulation using an *H*-formulation finite element model

Loïc Quéval<fnref n="1" />, Kun Liu<fnref n="2" />, Wenjiao Yang<fnref n="2" />, Víctor M R Zermeño<fnref n="3" /> and Guangtong Ma<fnref n="2" />,<fnref n="4" />

E-mail: loic.queval@geeps.centralesupelec.fr, gtma@swjtu.edu.cn and victor.zermeno@kit.edu

Received 16 March 2018, revised 2 May 2018

Accepted for publication 16 May 2018

Published 21 June 2018

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

The modeling of superconducting magnetic bearings (SMBs) is of great significance for predicting and optimizing their levitation performance before construction. Although much effort has been made in this area, there still remains some space for improvements. Thus the goal of this work is to report a flexible, fast and trustworthy *H*-formulation finite element model. First the methodology for modeling and calibrating both bulk-type and stack-type SMBs is summarized. Then its effectiveness for simulating SMBs in 2D, 2D axisymmetric and 3D is evaluated by comparison with measurements. In particular, original solutions to overcome several obstacles are given: clarification of the calibration procedure for stack-type and bulk-type SMBs, details on the experimental protocol to obtain reproducible measurements, validation of the 2D model for a stack-type SMB modeling the tapes' real thickness, implementation of a 2D axisymmetric SMB model, implementation of a 3D SMB model, and extensive validation of the models by comparison with experimental results for field cooling and zero field cooling, for both vertical and lateral movements. The accuracy of the model being having proven, it now has a strong potential for speeding up the development of numerous applications including maglev vehicles, magnetic launchers, flywheel energy storage systems, motor bearings and cosmic microwave background polarimeters.

**Keywords:** superconducting magnetic bearing (SMB), maglev, *H*-formulation finite element model, modeling and simulation

(Some figures may appear in colour only in the online journal)

## 1. Introduction

The relative movement between a permanent magnet (PM) and a high temperature superconductor (HTS) can induce supercurrents in the HTS. By interacting with the PM static magnetic field, these supercurrents produce a force that can be attractive

or repulsive depending on the arrangement and on the operating conditions. It can even provide passive stable levitation. This unique feature motivated the development of superconducting magnetic bearings (SMBs) [1–4]. They have been customized for numerous applications, including maglev vehicles [5–7], magnetic launchers [8, 9], flywheel energy storage systems [10–17], motors [18], and cosmic microwave background polarimeters [19–21].

<footnote n="1">¹Group of electrical engineering-Paris (GeePs), CNRS UMR 8507, CentraleSupélec, UPSud, UPMC, Gif-sur-Yvette, France</footnote>

<footnote n="2">²Applied Superconductivity Laboratory (ASCLab), State Key Laboratory of Traction Power, Southwest Jiaotong University, Chengdu, Sichuan 610031, People's Republic of China</footnote>

<footnote n="3">³Karlsruhe Institute of Technology, Hermann-von-Helmholtz Platz 1, D-76344 Eggenstein-Leopoldshafen, Germany</footnote>

<footnote n="4">⁴ Author to whom any correspondence should be addressed.</footnote>
<table><caption>**Table 1.** SMB finite element models.</caption><thead><tr><th></th><th>2D</th><th>2D axi</th><th>3D</th></tr></thead><tbody><tr><td>A-V</td><td>Homemade<br>Hofmann et al [25]<br>Dias et al [28–30]<br>Ma et al [32–34]</td><td>Sugiura et al [26]<br>Takeda et al [31]<br>Chun et al [35]<br>Ruiz-Alonso et al [36]<br>Wang et al [37]<br>Sotelo et al [38]</td><td>Ueda et al [27]</td></tr><tr><td></td><td>Software</td><td>Li et al [39]</td><td>Hauser [40]</td></tr><tr><td>T-Ω</td><td>Homemade<br>Zhang et al [41]</td><td>Zheng et al [42]<br>Gou et al [45]</td><td>Uesaka et al [43, 44]<br>Tsuchimoto et al [46]<br>Tsuda et al [47–49]<br>Ma et al [50, 51]<br>Pratap et al [52]</td></tr><tr><td></td><td>Software</td><td>—</td><td>—</td></tr><tr><td>E</td><td>Homemade</td><td>—</td><td>—</td></tr><tr><td></td><td>Software</td><td>—</td><td>—</td></tr><tr><td>H</td><td>Homemade<br>Lu et al [53]</td><td>—</td><td>Lu et al [54]<br>Yu et al [55]</td></tr><tr><td></td><td>Software<br>Sass et al [56]<br>Quéval et al [58]<br>This work</td><td>Patel et al [57]<br>This work</td><td>Patel et al [57]<br>Quéval et al [58]<br>This work</td></tr></tbody></table>

There are various analytical and numerical models that are able to predict, more or less accurately, the maglev performances of SMBs. A detailed review is provided by Navau *et al* in [22]. Among them, finite element (FE) models using various formulations are being intensively developed. The formulations are named after the state variables to be solved: A-V-formulation for the magnetic potential vector and the electric potential, T-Ω-formulation for the current potential vector and the magnetic potential, E-formulation for the electric field and H-formulation for the magnetic field. The critical state model [23] or the E-J power law model [24] is then commonly used together with one of these formulations to model the nonlinear resistivity of the superconductor. A summary of the formulation used by independent groups to model SMBs with homemade FE codes and free/paid for FE software is proposed in table 1.

All these models have their own features and limitations. Focusing on the H-formulation, important efforts have been made to simulate SMBs using homemade codes. Lu *et al* wrote a FE code in FORTRAN to estimate the levitation force between a PM and an HTS bulk in 2D [53]. This is probably an evolution of the code reported in [54] for the 3D simulation of a cylindrical HTS bulk over a PM guideway. In those articles, the field of the moving PM, obtained analytically, was applied as a time-dependent Dirichlet boundary condition on the outer boundary of a model including only the HTS domain and a thin air domain. But it is not clear if the self-field of the HTS bulk was included. The model and its extensions to other PM guideways geometries, field cooling and lateral movements [59–61] provided interesting guidelines but the authors provided no convincing experimental validation of it. Yu *et al* implemented a similar 3D model to

analyze a SMB made of a cylindrical PM and a cylindrical HTS bulk [55]. A substantial effort was made there to experimentally validate the model for both zero field cooling and field cooling, but only for vertical displacements. Surprisingly, the simulated levitation force did not go back to zero when the gap increased. And the levitation force loop proved difficult to reproduce for the field cooling case.

FE software has also been employed to simulate SMBs using the H-formulation. Actually the groups listed in table 1 all used COMSOL Multiphysics [62], either with the magnetic field formulation (mfh) physic available in the AC/DC module, or by manually implementing the partial differential equations (PDEs) with the PDE module. Sass *et al* developed a 2D model [56] to obtain the levitation force between a PM and an YBCO bulk or stacks of YBCO tapes. The field of the PM was obtained using analytical equations. To model the movement, the field generated by the PM was applied as a time-dependent Dirichlet boundary condition on a boundary close to the HTS domain. To reduce the computing time, a symmetry axis was used, restricting the movement to vertical displacements. To model the stacks, an anisotropic homogenized model was adopted [63]. The agreement with measurements for field cooling and zero field cooling was good. A similar model was developed by Quéval *et al* [58] to include the PM assembly real geometry and the iron nonlinearity. To do so, the field of the PM assembly was obtained using a magnetostatic FEM. Besides, the model was able to deal with any relative movement, making it possible to optimize the SMB on a realistic displacement sequence. A similar 3D model was mentioned in [58] but without details about its implementation. Patel *et al* introduced a 2D
axisymmetric H-formulation FEM in [57] to estimate the levitation force between a PM and stacks of YBCO tapes. The PM was modeled by a thin current domain approximating the ideal equivalent 2D axisymmetric current sheet. To model the movement, this thin domain was moved along the z-direction by defining it with a time and space dependent current density. With this modeling strategy, the boundary conditions are fixed but many elements are required to mesh the ‘moving’ PM assembly thus limiting the applicability of the model to simple geometries. To model the stacks, an isotropic homogenized model was used. The simulated levitation force, limited to the first magnetization, was compared with measurements from 20 to 77 K in a field cooling condition only. The agreement for a SMB with a rolled stack was fair at 20 K and reasonable at 77 K [57]. For a SMB with a stack of annuli [64], the agreement was good. Similarly, a 3D model was built to study the current pattern for the SMB with the rolled stack with limited discussion and validation [57].

The motivation behind this work is to develop flexible, fast and trustworthy H-formulation FE models able to predict the maglev performances of SMBs in 2D, 2D axisymmetric and 3D configurations. Key advancements with respect to previous models include: clarification of the calibration procedure for stack-type and bulk-type SMBs, details on the experimental protocol to obtain reproducible measurements, validation of the 2D model for a stack-type SMB considering the tapes real thickness, implementation of a 2D axisymmetric SMB model, implementation of a 3D SMB model, and extensive validation of the models by comparison with experimental results for field cooling and zero field cooling, for both vertical and lateral movements. The test cases reported here have been selected to serve as benchmarks, with the aim to help focus the effort of the numerical modeling community towards the most relevant approaches [65].

## 2. Superconducting magnetic bearing model

The SMB model is built by unidirectional coupling between the PM assembly model and the HTS assembly model. The coupling is done by applying the sum of the external field $H_{\text{ext}}$ and the self-field $H_{\text{self}}$ on the outer boundaries $\Gamma$ of the HTS assembly model (figure 1).

### 2.1. PM assembly model

The PM assembly is an arrangement of any number of PMs and ferromagnetic pieces surrounded by air (or any coolant). For simple geometries, analytical formulas could be used [56, 66]. But it is modeled here using a magnetostatic A-formulation FE model. This allows us to include the iron nonlinear B-H curve and to consider complex PM assembly geometries [58].

### 2.2. HTS assembly model

The HTS assembly is an arrangement of any number of normal and superconducting pieces (bulks or conductors) surrounded by air (or any coolant). It is assumed that the

<figure><img src="image_4.png" /><figcaption>Figure 1. Modeling approach. (a) PM assembly model. O is the origin of the coordinate system. The outer boundary is not shown. The time-dependent coordinates of M in the PM assembly reference frame describe the relative movement of the HTS and PM assemblies. (b) HTS assembly model. M is the origin of the coordinate system.</figcaption></figure>

materials are non-magnetic. To mathematically model the HTS assembly, the H-formulation is used [67, 68],

$$ \nabla \times \rho \nabla \times \mathbf{H} = -\mu_0 \frac{\partial \mathbf{H}}{\partial t} \quad \text{in } \Omega, \quad (1) $$

$$ \mathbf{H} = \mathbf{H}_{\text{self}} + \mathbf{H}_{\text{ext}} \text{ on } \Gamma, \quad (2) $$

$$ \mathbf{H}(t_0) = \mathbf{H}_0 | \nabla \cdot (\mu_0 \mathbf{H}_0) = 0, \quad (3) $$

where $\mathbf{H}$ is the magnetic field strength, $\rho$ is the material resistivity, $\mu_0$ is the vacuum magnetic permeability, $\Omega$ is the computational domain and $\Gamma$ is the outer domain boundary. Neumann boundary conditions are used for inner boundaries. On $\Gamma$, Dirichlet boundary conditions are used to impose the self-field $H_{\text{self}}$ (the one created by the supercurrent) and the external field $H_{\text{ext}}$ (the one created by the PM assembly). The current density $\mathbf{J}$, the electric field $\mathbf{E}$ and the magnetic flux density $\mathbf{B}$ can be obtained from $\mathbf{H}$ using,

$$ \mathbf{J} = \nabla \times \mathbf{H}, \quad (4) $$

$$ \mathbf{E} = \rho \mathbf{J}, \quad (5) $$

$$ \mathbf{B} = \mu_0 \mathbf{H}. \quad (6) $$

The resistivity $\rho_{\text{sc}}$ of the HTS is represented by a power law,

$$ \rho_{\text{sc}}(|\mathbf{J}|, \mathbf{B}) = \frac{E_c}{J_c(\mathbf{B})} \left| \frac{\mathbf{J}}{J_c(\mathbf{B})} \right|^{n-1}, \quad (7) $$

where $J_c$ is the field dependent local critical current density, $E_c$ is the critical current criterion and $n$ is a material parameter. To impose a transport current in a conductor, an integral constraint on the current density can be used

$$ I_{\text{tr}}(t) = \int_{\Omega_c} \mathbf{J} \cdot \mathrm{d}\mathbf{s}, \quad (8) $$

where $\Omega_c$ is the conductor cross section and $\mathrm{d}\mathbf{s}$ is the differential cross-sectional area vector. For the FE discretization, we use linear edge elements [67].
<figure><img src="image_5.png" /><figcaption>Figure 2. Test rig used to measure the 3D forces of the SMBs. The PM assembly was fixed to the base moving in the xy-direction. The HTS assembly was fixed to the 3D load cell moving in the z-direction.</figcaption></figure>

The external field $\mathbf{H}_{\text{ext}}$ is obtained from the PM assembly model. The static magnetic field generated by the PM assembly $\mathbf{H}_{\text{PM}}$ needs to be modified to take the relative movement into account. This is done by

$$ \mathbf{H}_{\text{ext}}(x, y, z, t) = T_t \mathbf{H}_{\text{PM}}(x, y, z), \quad (9) $$

where $T_t$ is the translation operator that describes the time-dependent position of the HTS assembly in the PM assembly reference frame.

The HTS is said to be ‘field cooled’ (FC) when the cooling is achieved close to the PM assembly, and ‘zero field cooled’ (ZFC) when the PM is far enough so that the applied field is negligible. We assume that during the cooling all the flux is pinned [69] and that no macroscopic currents are induced in the HTS [70]. This is experimentally validated by the fact that the forces after cooling but before any movements are null [29]. To simulate the FC case, we can therefore disregard $\mathbf{H}_{\text{self}}$ and set $\mathbf{H}_0 = T_{t0} \mathbf{H}_{\text{PM}}(x, y, z)$. By doing so, equation (3) is respected because the divergence of the field generated by the PM is zero. Note that we implicitly make here the hypothesis that the field generated by the super-current does not influence the PM’s remanent field.

The self-field $\mathbf{H}_{\text{self}}$ is obtained from the HTS assembly model at each time step by numerical integration of the Biot-Savart law. The consideration of $\mathbf{H}_{\text{self}}$ is required to make the problem self-consistent since the air/coolant layer around the HTS domain is slim. Indeed $\mathbf{H}_{\text{ext}}$ is applied on a boundary that is close to the HTS domain.

The force $\mathbf{F}$ (in N) between the PM assembly and the HTS assembly is obtained with

$$ \mathbf{F} = \int_{\Omega_{sc}} \mathbf{J} \times \mathbf{B} \, \mathrm{d}\mathbf{s}, \quad (10) $$

where $\Omega_{sc}$ is the HTS assembly cross section and $\mathbf{d}\mathbf{s}$ the differential cross-sectional area.

## 3. Measurements

The force measurements were carried out using a test rig developed at ASCLab (figure 2). The 3D relative motion is

<figure><img src="image_6.png" /><figcaption>Figure 3. SMB geometry and mesh for the 2D case. The point O is located at the center of the PM assembly top surface. The point M is located at the center of the stack bottom surface. The dimension of the PM and HTS assemblies in the x-direction are 240 mm and 100 mm, respectively. The arrows indicate the PM magnetization direction. The mesh of the air/coolant is not shown. Inset: zoom on the stack; the blue lines show the superconductor layers.</figcaption></figure>

obtained by three step motors and screw rods. The 3D position is recorded by three linear displacement sensors. The 3D force is measured by a 3D load cell. The time, the 3D position and the 3D force are recorded at 1 kS s⁻¹. The measured data presented here corresponds to a 500 point moving average.

The PM assembly is at room temperature while the HTS assembly is at liquid nitrogen temperature. A 1 mm sheet of aerogel paper CT200-Z is used to thermally insulate the PM and avoid a shift of its remanent flux density with the temperature during the measurement [71]. The z-direction force recorded by the load cell includes the weight of the HTS assembly: therefore the initial force (i.e. the weight) was subtracted from the measurements to remove any force not produced by the supercurrent in the measured data presented here. The liquid nitrogen container is mounted so that its weight is not measured by the load cell.

## 4. 2D case: linear SMB

### 4.1. Geometry

The linear SMB and the coordinate system adopted in this section are shown in figure 3. The PM assembly is made of
cuboidal Nd-Fe-B PMs and iron slabs arranged in flux concentration. The HTS assembly is a stack of 120 YBCO tapes (SuperPower SCS12050-AP). The HTS assembly can only move along the yz-plane ($x_M(t) = 0$).

## 4.2. Sequences

In this section, we consider three displacement sequences. They are described by the successive positions of $M(y_M, z_M)$ relative to O (in millimeters). The first position of each sequence is the cooling position. The moving speed is $1 \text{ mm s}^{-1}$ representing a quasistatic process.

* ZFC100: $(y_M, z_M) = \{(0, 100), (0, 6), (0, 100)\}$

* FC25: $(y_M, z_M) = \{(0, 25), (0, 6), (0, 25)\}$

* FC25_LD: $(y_M, z_M) = \{(0, 25), (0, 6), (6, 6), (-6, 6), (6, 6)\}$.

## 4.3. Modeling

Equations (1)–(10) are implemented in COMSOL Multi-physics 4.3a PDE mode application in a 2D space. More details about such implementation can be found in [72] for example. The HTS assembly is a stack of YBCO tapes: we model only the superconducting layers taking their real thickness into account. Each tape has a net current enforced to zero by means of an integral constraint. An anisotropic Kim-like model [73] is used to describe the dependence of the critical current density on the magnetic field,

$$J_c(\mathbf{B}) = \frac{J_{c0}}{\left(1 + \frac{\sqrt{k^2 B_{\parallel}^2 + B_{\perp}^2}}{B_0}\right)^{\alpha}}, \quad (11)$$

where $B_{\parallel}$ and $B_{\perp}$ are the field components parallel and perpendicular to the tape, respectively. $J_{c0}$, $B_0$, $k$ and $\alpha$ are material parameters. Equation (11) provides a reasonable description of the anisotropic behavior of HTS coated conductors (without artificial pinning) [74]. To mesh the superconducting layer, we use a mapped mesh [75] with ten elements distributed symmetrically following an arithmetic sequence in the width and one element in the thickness. Such mesh proved to be a good compromise between speed and accuracy. The outer boundary of the HTS assembly model $\Gamma$ is located at a distance of 1.5 mm from the HTS stack. This is less than the minimum levitation gap so that the coupling boundary is always inside the air gap.

From equation (9) in 2D, with the conventions of figure 2, the expression for $H_{\text{ext}}$ becomes

$$H_{\text{ext}}(y, z, t) = H_{\text{PM}}(y + y_M(t), z + z_M(t)), \quad (12)$$

where $(y_M, z_M)$ is the time-dependent position of the HTS assembly relative to O. $H_{\text{self}}$ is obtained by 2D integration of

<figure><img src="image_7.png" /><figcaption>Figure 4. 2D model calibration: magnetic flux density at 2 and 5 mm above the PM.</figcaption></figure>

the Biot-Savart law,

$$H_{\text{self},y}(y, z, t) = \frac{1}{2\pi} \iint_{\Omega_{\text{sc}}} \frac{-J_x(y', z', t) \cdot (z - z')}{(y - y')^2 + (z - z')^2} dy'dz', \quad (13)$$

$$H_{\text{self},z}(y, z, t) = \frac{1}{2\pi} \iint_{\Omega_{\text{sc}}} \frac{J_x(y', z', t) \cdot (y - y')}{(y - y')^2 + (z - z')^2} dy'dz', \quad (14)$$

where $\Omega_{\text{sc}}$ is the HTS assembly domain. This completes and corrects [56]. From equation (10) in 2D, the lateral force $F_y$ and the levitation force $F_z$ (in N) between the PM assembly and the HTS assembly are given by

$$F_y(t) = - \iint_{\Omega_{\text{sc}}} J_x(y', z', t) \cdot B_z(y', z', t) dy'dz' \cdot d_{\text{sc}}, \quad (15)$$

$$F_z(t) = \iint_{\Omega_{\text{sc}}} J_x(y', z', t) \cdot B_y(y', z', t) dy'dz' \cdot d_{\text{sc}}, \quad (16)$$

where $\Omega_{\text{sc}}$ is the HTS assembly domain and $d_{\text{sc}}$ is the dimension of the HTS assembly in the x-direction.

## 4.4. Model calibration

To calibrate the PM assembly model, we need to know the B-H curve of the iron and the remanent flux density $B_r$ of the PM. The assumed B-H curve is given in the appendix. To obtain the remanent flux density $B_r$ of the PM, we measured the magnetic flux density at several distances above the PM. By a trial-and-error process, we obtained $B_r$ that minimizes the difference between the measured data and the PM assembly model (figure 4). To calibrate the HTS assembly model, we need to get the values of five parameters: $J_{c0}$, $n$, $B_0$, $k$ and $\alpha$. To obtain $J_{c0}$, it is a common practice to use the maximum levitation force obtained for a zero field cooling sequence [29, 43, 51, 54]. The procedure used here is different. $J_{c0}$ and $n$ are obtained by fitting the power law to the measured current-voltage curve of a short sample of the same conductor. The measurement was made at 77 K using the four probe method. The other HTS tape parameters $B_0$, $k$ and $\alpha$ are obtained by trial-and-error so that the simulated maximum
<figure><img src="image_8.png" /><figcaption>**Figure 5.** 2D model calibration: levitation force for the ZFC100 sequence.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** 2D model validation: levitation force for the FC25 sequence.</figcaption></figure>

<table><caption>**Table 2.** Parameters for the 2D case.</caption><thead><tr><th>Symbol</th><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td>B<sub>r</sub></td><td>PM remanent flux density</td><td>1.12 T (side)<br>0.975 T (middle)</td></tr><tr><td>E<sub>c</sub></td><td>Critical current criterion</td><td>1 &times; 10<sup>-4</sup> V m<sup>-1</sup></td></tr><tr><td>n</td><td>HTS parameter</td><td>31</td></tr><tr><td>J<sub>c0</sub></td><td>HTS parameter</td><td>3.225 &times; 10<sup>10</sup> A m<sup>-2</sup></td></tr><tr><td>B<sub>0</sub></td><td>HTS parameter</td><td>0.0525 T</td></tr><tr><td>k</td><td>HTS parameter</td><td>0.256</td></tr><tr><td>&alpha;</td><td>HTS parameter</td><td>0.58</td></tr><tr><td>&rho;<sub>air</sub></td><td>Air resistivity</td><td>1 &Omega; m [76]</td></tr><tr><td>&mu;<sub>0</sub></td><td>Air/HTS permeability</td><td>4&pi; &times; 10<sup>-7</sup> H m<sup>-1</sup></td></tr></tbody></table>

levitation force during the ZFC100 sequence is equal to the measured value (figure 5). The procedure used here is applicable for any stack-type SMB with the advantage that only three parameters are obtained by trial-and-error. The parameters of the 2D case are summarized in table 2.

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** 2D model validation: (a) levitation force and (b) lateral force for the FC25_LD sequence.</figcaption></figure>

### 4.5. Model validation

To validate the 2D model, we consider the FC25 and FC25_LD sequences. The force calculated with the 2D model is in good agreement with the measured force (figures 6 and 7). This validates the modeling approach adopted. Similar simulations (not reported here) were performed for the stack-type SMB of [56] giving similar agreements with the measurements, and comparing well with the anisotropic homogenized bulk SMB model [63] implemented by Sass *et al.* Note that, by obtaining $H_{ext}$ using a magnetostatic FEM instead of analytical formulas, and by considering both vertical and lateral displacements, the current model overcomes the limitations of the previous one. Besides, by modeling here the stack with the tapes real thickness, we open the possibility of considering complex HTS assemblies.

## 5. 2D axisymmetric case: axisymmetric SMB

### 5.1. Geometry

The axisymmetric SMB and the coordinate system adopted in this section are shown in figure 8. The PM assembly is a cylindrical Nd-Fe-B magnet. The HTS assembly is a
<figure><img src="image_12.png" /><figcaption>Figure 8. SMB geometry and mesh for the 2D axisymmetric case. The point O is located at the center of the PM top surface. The point M is located at the center of the bulk bottom surface. The arrow indicates the PM magnetization direction. The mesh of the air/ coolant is not shown.</figcaption></figure>

cylindrical single domain melt-textured YBCO bulk. The HTS assembly can only move along the z-direction ($r_M(t) = 0$).

## 5.2. Sequences

In this section, we consider three displacement sequences. They are described by the successive positions of $M(r_M, z_M)$ relative to O (in millimeters). The first position of each sequence is the cooling position. The moving speed is $1 \text{ mm s}^{-1}$ representing a quasistatic process.

* ZFC100: ($r_M, z_M$) = { (0, 100), (0, 5), (0, 100) }

* FC25: ($r_M, z_M$) = { (0, 25), (0, 5), (0, 100), (0, 5) }

* FC5: ($r_M, z_M$) = { (0, 5), (0, 100), (0, 5), (0, 100) }.

## 5.3. Modeling

Equations (1)–(10) are implemented in COMSOL Multiphysics 4.3a PDE mode application in a 2D axisymmetric space. More details about such implementation can be found in [77] for example. The HTS assembly is a bulk, thus an isotropic Kim-like model [73] is used to describe the dependence of the critical current density on the magnetic field,

$$J_c(\mathbf{B}) = \frac{J_{c0}}{1 + \frac{|\mathbf{B}|}{B_0}}, \qquad (17)$$

where $J_{c0}$ and $B_0$ are material parameters. To mesh the HTS bulk, we use the mapped mesh shown in figure 8 with $8 \times 8$ elements distributed following arithmetic sequences in the $rz$-plane. The outer boundary of the HTS assembly model $\Gamma$ is located at 2.5 mm from the HTS bulk, corresponding to half of the minimum levitation gap.

<figure><img src="image_13.png" /><figcaption>Figure 9. 2D axi/3D model calibration: magnetic flux density at 2 and 5 mm above the PM.</figcaption></figure>

From (9) in 2D axisymmetric, with the conventions of figure 8, the expression for $H_{\text{ext}}$ becomes,

$$H_{\text{ext}}(r, z, t) = H_{\text{PM}}(r, z + z_M(t)), \qquad (18)$$

where $z_M$ is the time-dependent position of the HTS assembly relative to O. $H_{\text{self}}$ is obtained by 2D axisymmetric integration of the Biot-Savart law,

$$H_{\text{self},r}(r, z, t) = \frac{1}{4\pi} \iint_{\Omega_{sc}} \frac{-J_{\theta}(r', z', t) \sqrt{m}}{\sqrt{r'r^3}} (z - z') \times \left[ K(m) - \frac{2-m}{2(1-m)} E(m) \right] dr'dz', \qquad (19)$$

$$H_{\text{self},z}(r, z, t) = \frac{1}{4\pi} \iint_{\Omega_{sc}} \frac{J_{\theta}(r', z', t) \sqrt{m}}{\sqrt{r'r^3}} r \times \left[ K(m) + \frac{m(r+r') - 2r}{2r(1-m)} E(m) \right] dr'dz', \qquad (20)$$

$$m = \frac{4rr'}{(r+r')^2 + (z-z')^2}, \qquad (21)$$

where $\Omega_{sc}$ is the HTS assembly domain, and $K$ and $E$ are the complete elliptic integrals of the first and second kind,

$$K(m) = \int_{0}^{\pi/2} \frac{d\alpha}{\sqrt{1 - m \sin^2 \alpha}}, \qquad (22)$$

$$E(m) = \int_{0}^{\pi/2} \sqrt{1 - m \sin^2 \alpha} d\alpha. \qquad (23)$$

From equation (10) in 2D axisymmetric, the levitation force $F_z$ (in N) between the PM assembly and the HTS assembly is given by

$$F_z(t) = - \iint_{\Omega_{sc}} J_{\theta}(r', z', t) \cdot B_r(r', z', t) 2\pi r' dr' dz', \qquad (24)$$

where $\Omega_{sc}$ is the HTS assembly domain.

## 5.4. Model calibration

To obtain the remanent flux density $B_r$ of the PM cylinder, we measured the magnetic flux density at several distances above the PM. By a trial-and-error process, we obtained $B_r$ that
<figure><img src="image_14.png" /><figcaption>**Figure 10.** 2D axi/3D model calibration: levitation force for the ZFC100 sequence.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** 2D axi/3D model validation: levitation force for the FC25 sequence.</figcaption></figure>

<table><caption>**Table 3.** Parameters for the 2D axisymmetric and 3D cases.</caption><thead><tr><th>Symbol</th><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td><i>B</i><sub>r</sub></td><td>PM remanent flux density</td><td>1.27 T</td></tr><tr><td><i>E</i><sub>c</sub></td><td>Critical current criterion</td><td>1 &times; 10<sup>-4</sup> V m<sup>-1</sup></td></tr><tr><td><i>J</i><sub>c0</sub></td><td>HTS parameter</td><td>2.4 &times; 10<sup>8</sup> A m<sup>-2</sup></td></tr><tr><td><i>n</i></td><td>HTS parameter</td><td>21 [79]</td></tr><tr><td><i>B</i><sub>0</sub></td><td>HTS parameter</td><td>0.37 T</td></tr><tr><td>&rho;<sub>air</sub></td><td>Air resistivity</td><td>1 &Omega; m [76]</td></tr><tr><td>&mu;<sub>0</sub></td><td>Air/HTS permeability</td><td>4&pi; &times; 10<sup>-7</sup> H m<sup>-1</sup></td></tr></tbody></table>

minimizes the difference between the measured data and the PM assembly model (figure 9). To obtain $J_{c0}$, it is common practice to use the maximum levitation force obtained for a zero field cooling sequence [29, 36, 51, 54]. Accordingly, the critical current density $J_{c0}$ is set here at $2.4 \times 10^8$ A m$^{-2}$, so that the simulated maximum levitation force during the ZFC100 sequence is equal to the measured value (figure 10). Alternatively, $J_{c0}$ could have been determined beforehand as done for the 2D case, for example by measuring it by VSM (vibrating sample magnetometer) for a small piece from the bulk as reported in [78]. The other HTS bulk parameters are

<figure><img src="image_16.png" /><figcaption>**Figure 12.** 2D axi/3D model validation: levitation force for the FC5 sequence.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** SMB geometry and mesh for the 3D case. The point O is located at the center of the PM top surface. The point M is located at the center of the bulk bottom surface. The arrow indicates the PM magnetization direction. The mesh of the air/coolant is not shown.</figcaption></figure>

set to commonly used values. Note that the value of *n* weakly affects the calculated results if higher than 15. The parameters of the 2D axisymmetric case are summarized in table 3.

### 5.5. Model validation

To validate the 2D axisymmetric model, we consider the FC25 and FC5 sequences. The force calculated with the 2D axisymmetric model is in good agreement with the measured force (figures 11 and 12). This serves as a validation.

## 6. 3D case

### 6.1. Geometry

The 3D SMB and the coordinate system adopted in this section are shown in figure 13. The SMB is the same as that for the 2D axisymmetric case but the HTS assembly can now move along any direction.
## 6.2. Sequences

In this section, we consider six displacement sequences. They are described by the successive positions of $M(x_M, y_M, z_M)$ relative to O (in millimeters). The first position of each sequence is the cooling position. The moving speed is $1 \text{ mm s}^{-1}$ representing a quasistatic process.

* ZFC100: $(x_M, y_M, z_M) = \{(0, 0, 100), (0, 0, 5), (0, 0, 100)\}$

* FC25: $(x_M, y_M, z_M) = \{(0, 0, 25), (0, 0, 5), (0, 0, 100), (0, 0, 5)\}$

* FC5: $(x_M, y_M, z_M) = \{(0, 0, 5), (0, 0, 100), (0, 0, 5), (0, 0, 100)\}$

* ZFC100_Y7.5: $(x_M, y_M, z_M) = \{(0, 7.5, 100), (0, 7.5, 5), (0, 7.5, 100)\}$

* ZFC100_Y15: $(x_M, y_M, z_M) = \{(0, 15, 100), (0, 15, 5), (0, 15, 100)\}$

* FC25_LD: $(x_M, y_M, z_M) = \{(0, 0, 25), (0, 0, 5), (0, 7.5, 5), (0, -7.5, 5), (0, 7.5, 5), (0, -7.5, 5), (0, 7.5, 5), (0, 0, 5)\}$.

The ZFC100, FC25 and FC5 sequences are similar to the 2D axisymmetric case. The ZFC100_Y7.5 and ZFC100_Y15 sequences are similar to the ZFC100 sequences but the HTS bulk is off-axis.

## 6.3. Modeling

Equations (1)–(10) are implemented in COMSOL Multi-physics 4.3a PDE mode application in a 3D space. More details about such implementation can be found in [80] for example. To mesh the HTS bulk, we swept the mesh shown in figure 8 following a 360° circular path to obtain the hexahedral mesh shown in figure 13. The outer boundary of the HTS assembly model $\Gamma$ is here again located at 2.5 mm from the HTS bulk, corresponding to half the minimum levitation gap.

From equation (9) in 3D, with the conventions of figure 13, the expression for $\mathbf{H}_{\text{ext}}$ becomes

$$ \mathbf{H}_{\text{ext}}(y, z, t) = \mathbf{H}_{\text{PM}}(x + x_M(t), y + y_M(t), z + z_M(t)), \quad (25) $$

where $(x_M, y_M, z_M)$ is the time-dependent position of the HTS assembly relative to O. $\mathbf{H}_{\text{self}}$ is obtained by 3D integration of the Biot-Savart law,

$$ H_{\text{self},x}(x, y, z, t) = \frac{1}{4\pi} \times \iiint_{\Omega_{\text{sc}}} \frac{J_y(z - z') - J_z(y - y')}{\sqrt{(x - x')^2 + (y - y')^2 + (z - z')^2}} dx'dy'dz', \quad (26) $$

$$ H_{\text{self},y}(x, y, z, t) = \frac{1}{4\pi} \times \iiint_{\Omega_{\text{sc}}} \frac{J_z(x - x') - J_x(z - z')}{\sqrt{(x - x')^2 + (y - y')^2 + (z - z')^2}} dx'dy'dz', \quad (27) $$

<figure><img src="image_18.png" /><figcaption>Figure 14. 3D model validation: (a) levitation force and (b) lateral force for the ZFC100, ZFC100_Y7.5 and ZFC100_Y15 sequences.</figcaption></figure>

$$ H_{\text{self},z}(x, y, z, t) = \frac{1}{4\pi} \times \iiint_{\Omega_{\text{sc}}} \frac{J_x(y - y') - J_y(x - x')}{\sqrt{(x - x')^2 + (y - y')^2 + (z - z')^2}} dx'dy'dz', \quad (28) $$

where $\Omega_{\text{sc}}$ is the HTS assembly domain. From equation (10) in 3D, the forces $F_x$, $F_y$ and $F_z$ (in N) between the PM assembly and the HTS assembly are given by

$$ F_x(t) = \iiint_{\Omega_{\text{sc}}} J_y(x', y', z', t) \cdot B_z(x', y', z', t) \\ - J_z(x', y', z', t) \cdot B_y(x', y', z', t) dx'dy'dz', \quad (29) $$

$$ F_y(t) = \iiint_{\Omega_{\text{sc}}} J_z(x', y', z', t) \cdot B_x(x', y', z', t) \\ - J_x(x', y', z', t) \cdot B_z(x', y', z', t) dx'dy'dz', \quad (30) $$

$$ F_z(t) = \iiint_{\Omega_{\text{sc}}} J_x(x', y', z', t) \cdot B_y(x', y', z', t) \\ - J_y(x', y', z', t) \cdot B_x(x', y', z', t) dx'dy'dz', \quad (31) $$

where $\Omega_{\text{sc}}$ is the HTS assembly domain.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 15.** 3D model validation: (a) levitation force and (b) lateral force for the FC25_LD sequence.</figcaption></figure>

## 6.4. Model calibration

We use the same parameters as that for the 2D axisymmetric case (table 3).

## 6.5. Model validation

The 3D model should be able to reproduce the results obtained with the 2D axisymmetric model for the ZFC100, FC25 and FC5 sequences. The levitation force calculated with the 3D model has been added to figures 10–12, showing similar results. To further validate the 3D model, we consider the ZFC100_Y7.5 and ZFC100_Y15 sequences. The levitation and lateral forces calculated with the 3D model are in fair agreement with the measured force (figure 14). The calculated forces are somewhat smaller than the measured ones, but globally the force reduction as a function of the off-axis position is predicted correctly. Similar results have been obtained for a field cooling height of 5 mm (not reported here). Finally, we consider the FC25_LD sequence. The levitation and lateral forces calculated with the 3D model are plotted together with the measured data in figure 15. Note the instable behavior of the bearing: when the lateral position increases, the lateral force increases too. Here again, the agreement is fair considering the length of the sequence and

<table><caption>**Table 4.** Degree of freedom, time steps and computing time.</caption><thead><tr><th></th><th>DOF</th><th>Time steps</th><th>Computing time (s)</th></tr></thead><tbody><tr><td>2D ZFC100</td><td>20270</td><td>272139</td><td>392529</td></tr><tr><td>2D axi ZFC100</td><td>539</td><td>252</td><td>13038</td></tr><tr><td>3D ZFC100</td><td>21068</td><td>610</td><td>27777</td></tr></tbody></table>

the small amplitude of the lateral force. This validates the 3D model.

## 7. Discussion

The test cases considered above have been selected carefully to serve as benchmarks. For the 2D case, we selected a stack-type SMB for its true 2D nature. Indeed bulk-type SMB suffer from several factors that make them difficult to be simulated accurately in 2D. In particular, large bulks with homogeneous properties are difficult to obtain. The end effects and the impact of intragrain currents should then be taken into account [58, 81]. For the 2D axisymmetric and 3D cases, we selected a simplistic bulk-type SMB that allows comparison of the results for axial displacement sequences. Finally, we considered on purpose repetitive displacements. This is because simplified models, such as Meissner-limit and frozen-field models, can often estimate the first section of the force loop but generally fail to predict the rest [22, 57].

For the FE discretization, we use linear edge elements [67]. The degrees of freedom of the edge elements being associated with the tangential components along the edges of the elements, it is only possible to impose the tangential component of the field. Nevertheless, in practice a thin layer of air/coolant is sufficient to obtain accurate estimation of the maglev performances as demonstrated in this work.

Melt-textured YBCO bulks have an anisotropic critical current density: it is larger in the *ab*-plane than along the *c*-axis [82]. This is the reason why most of previous 3D SMB models used an anisotropic bulk model. This was either achieved by stacking multiple 2D layers [27, 42–44, 46–48], by superimposing two virtual HTS bulks [54] or by considering a tensor of resistivity [50]. As it is still not clear how to model HTS in 3D to include experimental phenomena such as flux cutting, flux flow and magnetically anisotropic critical current densities [65, 80, 83], we adopted here a simplistic isotropic bulk model. This probably explains the difference between simulation and measurements for the 3D sequences ZFC100_Y7.5, ZFC100_Y15 and FC25_LD. Indeed, for these sequences the bulk is off-axis and a current is induced along the *c*-axis. Nevertheless, the present results show that maglev performance of a bulk-type SMB can be reasonably well predicted using a 3D isotropic bulk model.

To simulate (zero) field cooling, we applied an initial field $H_0$ according to equation (3). But because of inherent numerical approximations (mesh, linear elements, etc), the curl of this field is not perfectly null. In 2D, this would be equivalent to apply a set of Dirichlet’s boundary conditions that does not satisfy the Ampère’s circuital law on the outer boundary $\Gamma$ [76].
<table><caption>**Table 5.** Relative error at the maximum force.</caption><thead><tr><th></th><th>Measured (N)</th><th>Model (N)</th><th>Relative error (%)</th></tr></thead><tbody><tr><td>2D FC25</td><td>110.7</td><td>114.6</td><td>3.5</td></tr><tr><td>2D axi FC25</td><td>41.1</td><td>42.8</td><td>4.1</td></tr><tr><td>3D FC25</td><td>41.1</td><td>42.4</td><td>3.2</td></tr></tbody></table>

As a result, some unphysical induced current might flow in the superconducting domain following equation (4) during (zero) field cooling. A fine mesh was selected here to limit this effect.

The computing time for each model depends on many factors such as: mesh quality, number of time steps, HTS parameters and displacement sequence. All the calculations were performed using COMSOL Multiphysics 4.3a [62] and a standard desktop computer (Intel i7-4770 s, 3.10 GHz, 8 GB RAM). The state variables were scaled to $10^7$, and the relative and absolute tolerances were set to $10^{-2}$ and $10^{-3}$, respectively. Table 4 gives a summary of the computational effort for the ZFC sequences. It can seem prohibitive for some applications, in particular when considering complex 3D SMB geometries. But we used here a rather fine mesh with the goal to obtain good agreements with measurements. As a result, the relative error at the maximum levitation force stayed below 5% for the FC25 sequences (table 5). Actually coarser meshes can often help to decrease the computing time to few seconds for 2D cases, without losing too much information [56, 58].

## 8. Conclusion

We reported here our experience on simulating SMBs with a commercial finite element software using the H-formulation in 2D, 2D axisymmetric and 3D. The main difficulty is linked to the task of modeling a moving magnet. To address this problem, we chose the approach consisting in modeling the movement via time-dependent Dirichlet boundary conditions. It requires (a) only one static solution of the PM assembly finite element model, and (b) a reduced air/coolant domain around the superconducting material in the HTS assembly model. With a proper calibration procedure, we showed that the proposed model can predict accurately the observed behavior of both stack-tape and bulk-type bearings, for various cooling conditions and various displacement sequences. This comprehensive validation is a necessary step before using such models for designing and optimizing realistic bearings. Besides, the test cases have been selected so that they could be used as a benchmark for other models.

Future efforts could be dedicated to reducing the computing time of such models. For stack-type bearings, the anisotropic homogenization proposed in [63] and extended in [84] is a good alternative. But it should be used with caution, and the first validations proposed in [56, 85] should be extended to other geometries and other test conditions. Another necessary step is the coupling of such models with motion equations, in order to predict the dynamic behavior of the loaded bearing. Indeed, here

the relative movement is the input of the simulation but in reality it is a consequence of the efforts exerted on the bearing [86]. Finally, further vetting and refining of the models could help developing and improving lumped parameter SMB models [87, 88], as a mean of drastically speeding up simulations.

## Acknowledgments

This work was supported in part by the National Natural Science Foundation of China under Grants 51475389, and 51722706, in part by the Fundamental Research Funds for the Central Universities under Grant 2682016ZY05, in part by the Sichuan Youth Science & Technology Foundation under Grant 2016JQ0003, and in part by the State Key Laboratory of Traction Power under Grant TPL1712.

## Appendix

### 2D case: iron B-H curve

$$ (B, H) = \{(0.0, 0.0), (0.5, 90.0), (1.0, 270.0), (1.1, 318.25), (1.2, 384.50), (1.3, 479.50), (1.3875, 608.562), (1.45, 755.437), (1.5, 939.185), (1.545, 1188.93), (1.575, 1407.93), (1.6275, 2077.31), (1.67375, 3117.93), (1.7025, 3969.37), (1.7275, 4843.66), (1.75825, 6081.34), (1.80875, 8581.09), (1.85, 11066.4), (1.9025, 14985.7), (2.05, 33003.3), (2.15, 59203.3), (2.22625, 93214.9), (2.27, 118884.0), (2.33375, 163558.0), (2.4075, 220788.0), (2.6, 373973.0), (3.0, 692281.0)\}. B in T, H in A m⁻¹.

## ORCID iDs

Loïc Quéval https://orcid.org/0000-0003-3934-4372

Guangtong Ma https://orcid.org/0000-0002-0249-4310

## References

[1] Moon F C and Chang P Z 1990 High-speed rotation of magnets on high $T_c$ superconducting bearings *Appl. Phys. Lett.* **56** 397–9

[2] Weinberger B R, Lynds L, Hull J R and Balachandran U 1991 Low friction in high temperature superconductor bearings *Appl. Phys. Lett.* **59** 1132–4

[3] Hull J R and Cansiz A 1999 Vertical and lateral forces between a permanent magnet and a high-temperature superconductor *J. Appl. Phys.* **86** 6396

[4] Werfel F N, Floegel-Delor U, Rothfeld R, Riedel T, Goebel B, Wippich D and Schirrmeister P 2012 Superconductor bearings, flywheels and transportation *Supercond. Sci. Technol.* **25** 014007

[5] Wang J S et al 2002 The first man-loading high temperature superconducting maglev test vehicle in the world *Physica C* **378–381** 809–14

[6] Schultz L, de Haas O, Verges P, Beyer C, Rohlig S, Olsen H, Kuhn L, Berger D, Noteboom U and Funk U 2005 Superconductively levitated transport system—the
SupraTrans project *IEEE Trans. Appl. Supercond.* **15** 2301–5

[7] Mattos L S, Rodriguez E, Costa F, Sotelo G G, de Andrade R and Stephan R M 2016 MagLev-Cobra operational tests *IEEE Trans. Appl. Supercond.* **26** 3600704

[8] Yang W J, Wen Z, Duan Y, Chen X D, Qiu M, Liu Y and Lin L Z 2006 Construction and performance of HTS maglev launch assist test vehicle *IEEE Trans. Appl. Supercond.* **16** 1108–11

[9] Wang J S, Wang S, Deng C, Zheng J, Song H, He Q and Zeng Y 2007 Laboratory-scale high temperature superconducting maglev launch system *IEEE Trans. Appl. Supercond.* **17** 2091–4

[10] Bornemann H, Ritter T, Urban C, Paitsev O, Peber K and Rietschel H 1994 Low friction in a flywheel system with passive superconducting magnetic bearings *IEEE Trans. Appl. Supercond.* **2** 439–47

[11] Chen Q Y, Xia Z, Ma K B, McMichael C K, Lamb M, Coolep R S, Fopler P C and Chu W K 1994 Hybrid high $T_c$ superconducting magnetic bearings for flywheel energy storage system *IEEE Trans. Appl. Supercond.* **2** 457–64

[12] Miyagawa Y, Kameno H, Takahata R and Ueyama H 1999 A 0.5 kWh flywheel energy storage system using a high-$T_c$ superconducting magnetic bearing *IEEE Trans. Appl. Supercond.* **9** 996–9

[13] Coombs T, Campbell A M, Storey R and Weller R 1999 Superconducting magnetic bearings for energy storage flywheels *IEEE Trans. Appl. Supercond.* **9** 968–71

[14] Ichihara T et al 2005 Application of superconducting magnetic bearings to a 10 kWh-class flywheel energy storage system *IEEE Trans. Appl. Supercond.* **15** 2245–8

[15] Werfel F N, Floegel-Delor U, Riedel T, Rothfeld R, Wippich D, Goebel B, Reiner G and Wehlau N 2007 A compact HTS 5 kWh/250 kW flywheel energy storage system *IEEE Trans. Appl. Supercond.* **17** 2138–41

[16] Strasik M et al 2007 Design, fabrication, and test of a 5 kWh/ 100 kW flywheel energy storage utilizing a high-temperature superconducting bearing *IEEE Trans. Appl. Supercond.* **17** 2133–7

[17] Mukoyama S et al 2017 Development of superconducting magnetic bearing for 300 kW flywheel energy storage system *IEEE Trans. Appl. Supercond.* **27** 3600804

[18] Nagaya K, Kosugi Y, Suzuki T and Murakami I 1999 Pulse motor with high-temperature superconducting levitation *IEEE Trans. Appl. Supercond.* **9** 4688–94

[19] Hull J R, Hanany S, Matsumura T, Johnson B and Jones T 2005 Characterization of a high-temperature superconducting bearing for use in a cosmic microwave background polarimeter *Supercond. Sci. Technol.* **18** S1–S5

[20] Matsumura T, Kataza H, Utsunomiya S, Yamamoto R, Hazumi M and Katayama N 2016 Design and performance of a prototype polarization modulator rotational system for use in space using a superconducting magnetic bearing *IEEE Trans. Appl. Supercond.* **26** 3602304

[21] EBEX Collaboration 2017 The EBEX balloon borne experiment—optics, receiver, and polarimetry *Astrophys. J. Suppl.* submitted (arXiv:1703.03847)

[22] Navau C, Del-Valle N and Sanchez A 2013 Macroscopic modeling of magnetization and levitation of hard type-II superconductors: the critical-state model *IEEE Trans. Appl. Supercond.* **23** 8201023

[23] Bean C P 1962 Magnetization of hard superconductors *Phys. Rev. Lett.* **8** 250–3

[24] Rhyner J 1993 Magnetic properties and AC-losses of superconductors with power law current–voltage characteristics *Physica C* **202** 292–300

[25] Hofmann C and Ries G 2001 Modelling the interactions between magnets and granular high-$T_c$ superconductor

material with a finite-element method *Supercond. Sci. Technol.* **14** 34–40

[26] Sugihara T, Hashizume H and Miya K 1991 Numerical electromagnetic field analysis of type-II superconductors *Int. J. Appl. Electromagn. Mater.* **2** 183–96

[27] Ueda H, Azumaya S, Tsuchiya S and Ishiyama A 2006 3D electro-magnetic analysis of levitating transporter using bulk superconductor *IEEE Trans. Appl. Supercond.* **16** 1092–5

[28] Dias D H N, Motta E S, Sotelo G G, de Andrade R Jr, Stephan R M, Kuehn L, de Haas O and Schultz L 2009 Simulations and tests of superconducting linear bearings for a maglev prototype *IEEE Trans. Appl. Supercond.* **19** 2120–3

[29] Dias D H N, Motta E S, Sotelo G G and de Andrade R Jr 2010 Experimental validation of field cooling simulations for linear superconducting magnetic bearings *Supercond. Sci. Technol.* **23** 075013

[30] Dias D H N, Sotelo G G and de Andrade R Jr 2011 Study of the lateral force behavior in a field cooled superconducting linear bearing *IEEE Trans. Appl. Supercond.* **21** 1533–37

[31] Takeda N, Uesaka M and Miya K 1994 Computation and experiments on the static and dynamic characteristics of high $T_c$ superconducting levitation *Cryogenics* **34** 745–52

[32] Ma G T 2013 Considerations on the finite-element simulation of high-temperature superconductors for magnetic levitation purposes *IEEE Trans. Appl. Supercond.* **23** 3601609

[33] Ma G T, Liu H, Li X T, Zhang H and Xu Y Y 2014 Numerical simulations of the mutual effect among the superconducting constituents in a levitation system with translational symmetry *J. Appl. Phys.* **115** 083908

[34] Ye C Q, Ma G T and Wang J S 2016 Calculation and optimization of high-temperature superconducting levitation by a vector potential method *IEEE Trans. Appl. Supercond.* **26** 3603309

[35] Chun Y D, Kim Y H, Lee J, Hong J P and Lee J W 2001 Finite element analysis of magnetic field in high temperature bulk superconductor *IEEE Trans. Appl. Supercond.* **11** 2000–3

[36] Ruiz-Alonso D, Coombs T A and Campbell A M 2004 Numerical analysis of high-temperature superconductors with the critical-state model *IEEE Trans. Appl. Supercond.* **14** 2053–63

[37] Wang L, Wang H H and Wang Q L 2006 Finite element analysis of magnetic levitation force in superconducting magnetic levitation system *Cryog. Supercond.* **34** 190–3 (in Chinese)

[38] Sotelo G G, de Andrade R Jr and Ferreira A C 2009 Test and simulation of superconducting magnetic bearings *IEEE Trans. Appl. Supercond.* **19** 1681–6

[39] Li Y L, Fang J, Guo M Z, Xiao L, Zheng M H and Jiao Y L 2008 ANSYS-based analysis of levitation force in the HTS hybrid magnetic bearings *Cryog. Supercond.* **36** 40–4 (in Chinese)

[40] Hauser A O 1997 Calculation of superconducting magnetic bearings using a commercial FE-program (ANSYS) *IEEE Trans. Magn.* **33** 1572–5

[41] Zhang J, Zeng Y, Cheng J and Tang X 2008 Optimization of permanent magnet guideway for HTS maglev vehicle with numerical methods *IEEE Trans. Appl. Supercond.* **18** 1681–6

[42] Zheng X J and Yang Y 2007 Transition cooling height of high-temperature superconductor levitation system *IEEE Trans. Appl. Supercond.* **17** 3862–6

[43] Uesaka M, Yoshida Y, Takeda N and Miya K 1993 Experimental and numerical analysis of three-dimensional high-$T_c$ superconducting levitation systems *Int. J. Appl. Electromagn. Mater.* **4** 13–25
[44] Yoshida Y, Uesaka M and Miya K 1994 Magnetic field and force analysis of high $T_c$ superconductor with flux flow and creep *IEEE Trans. Magn.* **30** 3503–6

[45] Gou X F, Zheng X J and Zhou Y H 2007 Drift of levitated/ suspended body in high-$T_c$ superconducting levitation systems under vibration: I. A criterion based on magnetic force-gap relation for gap varying with time *IEEE Trans. Appl. Supercond.* **17** 3795–802

[46] Tsuchimoto M and Honma T 1994 Numerical evaluation of levitation force of HTSC flywheel *IEEE Trans. Magn.* **4** 211–5

[47] Tsuda M, Lee H and Iwasa Y 1998 Electromaglev (active-maglev)-magnetic levitation of a superconducting disk with a DC field generated by electromagnets: III. Theoretical results on levitation height and stability *Cryogenics* **38** 743–56

[48] Tsuda M, Lee H, Noguchi S and Iwasa Y 1998 Electromaglev (active-maglev)-magnetic levitation of a superconducting disk with a DC field generated by electromagnets: IV. Theoretical and experimental results on supercurrent distributions in field-cooled YBCO disks *Cryogenics* **39** 893–903

[49] Ueda H and Ishiyama A 2004 Dynamic characteristics and finite element analysis of a magnetic levitation system using a YBCO bulk superconductor *Supercond. Sci. Technol.* **17** S170–5

[50] Ma G T, Wang J S and Wang S Y 2010 3D modeling of high-$T_c$ superconductor for magnetic levitation/suspension application: I. Introduction to the method *IEEE Trans. Appl. Supercond.* **20** 2219–27

[51] Ma G T, Wang J S and Wang S Y 2010 3D modeling of high-$T_c$ superconductor for magnetic levitation/suspension application: II. Validation with experiment *IEEE Trans. Appl. Supercond.* **20** 2228–34

[52] Pratap S and Hearn C S 2015 3D transient modeling of bulk high-temperature superconducting material in passive magnetic bearing applications *IEEE Trans. Appl. Supercond.* **25** 5203910

[53] Lu Y and Qin Y 2015 Influence of critical current density on magnetic force of HTSC bulk above PMR with 3D- modeling numerical solutions *Int. J. Mod. Phys. B* **29** 1542038

[54] Lu Y Y, Wang J S, Wang S Y and Zheng J 2008 3D-modeling numerical solutions of electromagnetic behavior of HTSC bulk above permanent magnetic guideway *J. Supercond. Nov. Magn.* **21** 467–72

[55] Yu Z Q, Zhang G M, Qiu Q Q and Hu L 2015 Numerical simulation of levitation characteristics of a cylindrical permanent magnet and a high-temperature superconductor based on the 3D finite-element method *Trans. China Electrotech. Soc.* **30** 32–8 (in Chinese)

[56] Sass F, Sotelo G G, de Andrade R Jr and Sirois F 2015 H-formulation for simulating levitation forces acting on HTS bulks and stacks of 2G coated conductors *Supercond. Sci. Technol.* **28** 125012

[57] Patel A, Hopkins S C, Baskys A, Kalitka V, Molodyk A and Glowacki B A 2015 Magnetic levitation using high temperature superconducting pancake coils as composite bulk cylinders *Supercond. Sci. Technol.* **28** 115007

[58] L'Quéval L, Sotelo G G, Kharmiz Y, Dias D H N, Sass F, Zermeño V M R and Gottkehaskamp R 2016 Optimization of the superconducting linear magnetic bearing of a maglev vehicle *IEEE Trans. Appl. Supercond.* **26** 3601905

[59] Lu Y Y, Lu B J and Wang S Y 2011 The relationship of magnetic stiffness between single and multiple YBCO superconductors over permanent magnet guideway *J. Low Temp. Phys.* **164** 279–86

[60] Lu Y Y and Zhuang S J 2012 Magnetic forces simulation of bulk HTS over permanent magnetic railway with numerical method *J. Low Temp. Phys.* **169** 111–21

[61] Lu Y Y and Dang Q H 2012 Magnetic forces investigation of bulk HTS over permanent magnetic guideway under different lateral offset with 3D-model numerical method *Adv. Mater. Sci. Eng.* **2012** 640497

[62] COMSOL Multiphysics version 4.3a www.comsol.com

[63] Rodriguez-Zermeño V M, Abrahamsen A B, Mijatovic N, Jensen B B and Sørensen M P 2013 Calculation of alternating current losses in stacks and coils made of second generation high temperature superconducting tapes for large scale applications *J. Appl. Phys.* **114** 173901

[64] Patel A, Hahn S, Voccio J, Baskys A, Hopkins S C and Glowacki B A 2017 Magnetic levitation using a stack of high temperature superconducting tape annuli *Supercond. Sci. Technol.* **30** 024007

[65] Sirois F and Grilli F 2015 Potential and limits of numerical modelling for supporting the development of HTS devices *Supercond. Sci. Technol.* **28** 043002

[66] Ma G T, Liu H F, Wang J S, Wang S Y and Li X C 2009 3D modeling permanent magnet guideway for high temperature superconducting maglev vehicle application *J. Supercond. Novel Magn.* **22** 841–7

[67] Brambilla R, Grilli F and Martini L 2007 Development of an edge-element model for AC loss computation of high-temperature superconductors *Supercond. Sci. Technol.* **20** 16–24

[68] Rodriguez-Zermeño V M, Grilli F and Sirois F 2013 A full 3D time-dependent electromagnetic model for Roebel cables *Supercond. Sci. Technol.* **26** 052001

[69] Huang H, Zheng J, Zheng B T, Qian N, Li H T, Li J and Deng Z G 2017 Correlations between magnetic flux and levitation force of HTS bulk above a permanent magnet guideway *J. Low Temp. Phys.* **189** 42–52

[70] Navau C, Sanchez A, Pardo E and Chen D-X 2004 Equilibrium positions due to different cooling processes in superconducting levitation systems *Supercond. Sci. Technol.* **17** 828–32

[71] Lu Y Y, Qin Y J, Dang Q H and Wang J S 2010 Influence of experimental methods on crossing in magnetic force-gap hysteresis curve of HTS maglev system *Physica C* **470** 1994–7

[72] Hong Z, Campbell A M and Coombs T A 2006 Numerical solution of critical state in superconductivity by finite element software *Supercond. Sci. Technol.* **19** 1246–52

[73] Kim Y B, Hempstead C F and Strnad A R 1962 Critical persistent currents in hard superconductors *Phys. Rev. Lett.* **9** 306

[74] Grilli F, Sirois F, Rodriguez-Zermeño V M and Vojenčiak M 2014 Self-consistent modeling of the $I_c$ of HTS devices: how accurate do models really need to be? *IEEE Trans. Appl. Supercond.* **24** 8000508

[75] Rodriguez-Zermeño V M, Mijatovic N, Traeholt C, Zirngibl T, Seiler E, Abrahamsen A B, Pedersen N F and Sørensen M P 2011 Towards faster FEM simulation of thin film superconductors: a multiscale approach *IEEE Trans. Appl. Supercond.* **21** 3273–6

[76] Lahtinen V, Lyly M, Stenvall A and Tarhasaari T 2012 Comparison of three eddy current formulations for superconductor hysteresis loss modelling *Supercond. Sci. Technol.* **25** 115001

[77] Zhang M, Kvitkovic J, Pamidi S V and Coombs T A 2012 Experimental and numerical study of a YBCO pancake coil with a magnetic substrate *Supercond. Sci. Technol.* **25** 125020

[78] Sawamura M and Tsuchimoto M 2000 Numerical analysis for superconductor in sheet and bulk form *Japan J. Ind. Appl. Math.* **17** 199–208

[79] Ainslie M D and Fujishiro H 2015 Modelling of bulk superconductor magnetization *Supercond. Sci. Technol.* **28** 053002
[80] Zhang M and Coombs T A 2012 3D modeling of high-$T_c$ superconductors by finite element software Supercond. Sci. Technol. **25** 015009

[81] Deng Z et al 2012 Trapped flux and levitation properties of multiseeded YBCO bulks for HTS magnetic device applications: I. Grain and current features IEEE Trans. Appl. Supercond. **22** 6800110

[82] Murakami M, Oyama T, Fujimoto H, Gotoh S, Yamaguchi K, Shiohara Y, Koshizuaka N and Tanaka S 1991 Melt processing of bulk high $T_c$ superconductors and their application IEEE Trans. Magn. **27** 1479–86

[83] Badía-Majós A and López C 2012 Electromagnetics close beyond the critical state: thermodynamic prospect Supercond. Sci. Technol. **25** 104004

[84] Quéval L, Rodriguez-Zermeño V M and Grilli F 2016 Numerical models for AC loss calculation in large-scale

applications of HTS coated conductors Supercond. Sci. Technol. **29** 024007

[85] Liu K, Yang W, Ma G T, Quéval L, Gong T, Ye C, Li X and Luo Z 2017 Experiment and simulation of superconducting magnetic levitation with REBCO coated conductor stacks Supercond. Sci. Technol. **31** 015013

[86] Dias D H N, Sotelo G G, Rodriguez E F, de Andrade R Jr and Stephan R M 2013 Emulation of a full scale maglev vehicle behavior under operational conditions IEEE Trans. Appl. Supercond. **23** 3601105

[87] Hearn C S, Pratap S B, Chen D and Longoria R G 2014 Lumped-parameter model to describe dynamic translational interaction for high-temperature superconducting bearings IEEE Trans. Appl. Supercond. **24** 46–53

[88] Hearn C S, Pratap S B, Chen D and Longoria R G 2016 Dynamic performance of lumped parameter model for superconducting levitation IEEE Trans. Appl. Supercond. **26** 3602608