Superconducting Magnetic Bearings Simulation using an H-formulation Finite Element Model

The modeling of superconducting magnetic bearing (SMB) is of great significance for predicting and optimizing its levitation performance before construction. Although lots of efforts have been made in this area, it 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 SMB is summarized. Then its effectiveness for simulating SMBs in 2-D, 2-D axisymmetric and 3-D 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 2-D model for a stack-type SMB modeling the tapes real thickness, implementation of a 2-D axisymmetric SMB model, implementation of a 3-D SMB model, 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 proved, it has now 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.

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 [21]. 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 (CSM) [22] or the E-J power law model [23] 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/nonfree FE softwares is proposed in Table I.   TABLE I  SMB FINITE ELEMENT MODELS  2-D  2-D axi  3-D   A˗V  Homemade Hofmann et al. [24] Dias et al. [25][26][27] Ma et al. [28][29][30] Sugiura et al. [31] Takeda et al. [32] Chun et al. [33] Ruiz-Alonso et al. [34] Wang et al. [35] Sotelo et al. [36] Ueda et al. [38] Software -Li et al. [37] Hauser [39] T˗Ω Homemade Zheng et al. [40] Zhang et al. [41] Gou et al. [42] Uesaka et al. [43,44] Tsuchimoto et al. [45] Tsuda et al. [46][47][48] Ma et al. [49,50] Pratap et al. [ Sass et al. [53] Quéval et al. [54] *this work* Patel et al. [55] *this work* Patel et al. [55] Quéval et al. [54] *this work* 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 2-D [52]. This is probably an evolution of the code reported in [56] for the 3-D 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 [58][59][60] provided interesting guidelines but the authors provided no convincing experimental validation of it. Yu et al. implemented a similar 3-D model to analyze a SMB made of a cylindrical PM and a cylindrical HTS bulk [57]. 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 softwares have also been employed to simulate SMBs using the H-formulation. Actually the groups listed in Table I all used Comsol Multiphysics [61], 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 2-D model [53] 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 [62]. The agreement with measurements for field cooling and zero field cooling was good. A similar model was developed by Quéval et al. [54] 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 3-D model was mentioned in [54] but without details about its implementation. Patel et al. introduced a 2-D axisymmetric H-formulation FEM in [55] 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 2-D 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 K to 77 K in field cooling condition only. The agreement for a SMB with a rolled stack was fair at 20 K and reasonable at 77 K [55]. For a SMB with a stack of annuli [63], the agreement was good. Similarly, a 3-D model was built to study the current pattern for the SMB with the rolled stack with limited discussion and validation [55].
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 2-D, 2-D axisymmetric and 3-D configurations. Key advancement 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 2-D model for a stacktype SMB considering the tapes real thickness, implementation of a 2-D axisymmetric SMB model, implementation of a 3-D SMB model, 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 hope to help focus the effort of the numerical modelling community towards the most relevant approaches [64].

II. 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 and the self-field on the outer boundaries Γ of the HTS assembly model ( Figure 1). 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 [65,53]. 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 [54].

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 materials are non-magnetic. To mathematically model the HTS assembly, the H-formulation is used [66,67], where is the magnetic field strength, is the material resistivity, 0 is the vacuum magnetic permeability, Ω is the computational domain and Γ is the outer domain boundary. Neumann boundary conditions are used for inner boundaries. On Γ, Dirichlet boundary conditions are used to impose the self-field (the one created by the supercurrent) and the external field (the one created by the PM assembly). The current density , the electric field and the magnetic flux density can be obtained from using, The resistivity of the HTS is represented by a power law, where is the field dependent local critical current density, is the critical current criterion and is a material parameter. To impose a transport current in a conductor, an integral constraint on the current density can be used where Ω is the conductor cross section and d is the differential cross-sectional area vector. For the finite element discretization, we use linear edge elements [66].
The external field is obtained from the PM assembly model. The static magnetic field generated by the PM assembly needs to be modified to take the relative movement into account. This is done by ( , , , ) = ( , , ) where 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 [68] and that no macroscopic currents are induced in the HTS [69]. This is experimentally validated by the fact that the forces after cooling but before any movement are null [26]. To simulate the FC case, we can therefore disregard and set 0 = 0 (x, y, z). By doing so, (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 supercurrent does not influence the PM's remanent field.
The self-field is obtained from the HTS assembly model at each time step by numerical integration of the Biot-Savart law. The consideration of is required to make the problem self-consistent since the air/coolant layer around the HTS domain is slim. Indeed is applied on a boundary that is close to the HTS domain.
The force F between the PM assembly and the HTS assembly is obtained with where Ω is the HTS assembly cross section and ds the differential cross-sectional area.

Measurements
The force measurements were carried out using a test rig developed at ASCLab (Figure 2). The 3-D relative motion is obtained by three step motors and screw rods. The 3-D position is recorded by three linear displacement sensors. The 3-D force is measured by a 3-D load cell. The time, the 3-D position and the 3-D force are recorded at 1 kHz. The measured data presented here corresponds to a 500 points 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 [70]. 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. IV. 2D case: linear SMB 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 permanent magnets 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 ( ( ) = 0). More details about such implementation can be found in [71] 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 [72] is used to describe the dependence of the critical current density on the magnetic field, where B // and B ﬩ are the field components parallel and perpendicular to the tape, respectively. 0 , 0 , k and α are material parameters. Equation (11) provides a reasonable description of the anisotropic behavior of HTS coated conductors (without artificial pinning) [73]. To mesh the superconducting layer, we use a mapped mesh [74] with 10 elements distributed symmetrically following an arithmetic sequence in the width and 1 element in the thickness. Such mesh proved to be a good compromise between speed and accuracy. The outer boundary of the HTS assembly model Γ 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 (9) where ( , ) is the time-dependent position of the HTS assembly relative to O. is obtained by 2-D integration of the Biot-Savart law, , ( , , ) = where Ω is the HTS assembly domain. This completes and corrects [53].

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 of the PM. The assumed B-H curve is given in appendix. To obtain the remanent flux density of the PM, we measured the magnetic flux density at several distances above the PM. By a trial-and-error process, we obtained that minimize 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 5 parameters: 0 , n, 0 , k and α. To obtain 0 , it is a common practice to use the maximum levitation force obtained for a zero field cooling sequence [43,56,26,50]. The procedure used here is different. 0 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 0 , k and α are obtained by trial-and-error so that the simulated maximum 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 3 parameters are obtained by trial-and-error. The parameters of the 2-D case are summarized in Table II.     Figure 7). This validates the modeling approach adopted. Similar simulations were performed for the stack-type SMB analyzed in [53] giving similar agreements (not reported here). Note the four main differences between this model and the one previously reported in [53]. First, we obtain using a magnetostatic FEM while Sass et al. used analytical formulas. Second, we model the stack with tapes real thickness whereas Sass et al. used an anisotropic homogenized bulk model [62].  V.

2D axisymmetric case: axisymmetric SMB 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 cylindrical single domain melt-textured YBCO bulk manufactured by Beijing General Research Institute for Nonferrous Metals. The HTS assembly can only move along the z-direction ( ( ) = 0). More details about such implementation can be found in [76] for example. The HTS assembly is a bulk, thus an isotropic Kim-like model [72] is used to describe the dependence of the critical current density on the magnetic field, where 0 and 0 are material parameters. To mesh the HTS bulk, we use the mapped mesh shown in Figure 8 with 8×8 elements distributed following arithmetic sequences in the -plane. The outer boundary of the HTS assembly model Γ is located at 2.5 mm from the HTS bulk, corresponding to half of the minimum levitation gap.  (19) where Ω is the HTS assembly domain, and K and E are the complete elliptic integrals of the first and second kind, 4) Model calibration To obtain the remanent flux density of the PM cylinder, we measured the magnetic flux density at several distances above the PM. By a trial-and-error process, we obtained that minimizes the difference between the measured data and the PM assembly model (Figure 9). To obtain 0 , it is common practice to use the maximum levitation force obtained for a zero field cooling sequence [34,56,26,50]. Accordingly, the critical current density 0 is set here at 2.4•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, 0 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 [77]. The other HTS bulk parameters are set to commonly used values. Note that the value of n weakly affects the calculated results if higher than 15. The parameters of the 2-D axisymmetric case are summarized in Table III.

VI.
3D case 1) Geometry The 3-D SMB and the coordinate system adopted in this section are shown in Figure 13. The SMB is the same as that for the 2-D axisymmetric case but the HTS assembly can now move along any direction. The ZFC100, FC25 and FC5 sequences are similar to the 2-D axisymmetric case. The ZFC100_Y7.5 and ZFC100_Y15 sequences are similar to the ZFC100 sequences but the HTS bulk is off-axis.

3) Modeling
Equations (1)-(10) are implemented in COMSOL Multiphysics 4.3a PDE mode application in a 3-D space. More details about such implementation can be found in [79] 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 Γ is here again located at 2.5 mm from the HTS bulk, corresponding to half the minimum levitation gap.
From (9) where Ω is the HTS assembly domain.

4) Model calibration
We use the same parameters as that for the 2-D axisymmetric case (Table III).

5) Model validation
The 3-D model should be able to reproduce the results obtained with the 2-D axisymmetric model for the ZFC100, FC25 and FC5 sequences. The levitation force calculated with the 3-D model has been added to Figure 10, Figure 11 and Figure 12, showing similar results. To further validate the 3-D model, we consider the ZFC100_Y7.5 and ZFC100_Y15 sequences. The levitation and lateral forces calculated with the 3-D 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 FC5_LD sequence. The levitation and lateral forces calculated with the 3-D 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 the small amplitude of the lateral force. This validates the 3-D model.

VII. Discussion
The test cases considered above have been selected carefully to serve as benchmarks. For the 2-D case, we selected a stack-type SMB for its true 2-D nature. Indeed bulk-type SMB suffer from several factors that make them difficult to be simulated accurately in 2-D. 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 [80,54]. For the 2-D axisymmetric and 3-D 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 [21,55].
For the FE discretization, we use linear edge elements [66]. 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 caxis [81]. This is the reason why most of previous 3-D SMB models used an anisotropic bulk model. This was either achieved by stacking multiple 2-D layers [43-47, 38, 40], by superimposing two virtual HTS bulks [56] or by considering a tensor of resistivity [49]. As it is still not clear how to model HTS in 3-D to include experimental phenomena such as flux cutting, flux flow and magnetically anisotropic critical current densities [79,82,64], we adopted here a simplistic isotropic bulk model. This probably explains the difference between simulation and measurements for the 3-D 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 3-D isotropic bulk model.
To simulate (zero) field cooling, we applied an initial field 0 according to (3). But because of inherent numerical approximations (mesh, linear elements, etc.), the curl of this field is not perfectly null. In 2-D, this would be equivalent to apply a set of Dirichlet's boundary conditions that doesn't satisfy the Ampère's circuital law on the outer boundary Γ [75]. As a result, some unphysical induced current might flow in the superconducting domain following (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 [61] and a standard desktop computer (Intel i7-4770S, 3.10 GHz, RAM 8 GB). The state variables were scaled to 10 7 , and the relative and absolute tolerances were set to 10 -2 and 10 -3 , respectively. Table IV gives a summary of the computational effort for some sequences. It can seem prohibitive for some applications, in particular when considering complex 3-D SMB geometries. But we used here a rather fine mesh with the goal to obtain good agreements with measurements. Actually coarser meshes can often help to decrease the computing time to few seconds for 2-D cases, without losing too much information [53,54].

VIII. Conclusion
We reported here our experience on simulating superconducting magnetic bearing with a commercial finite element software using the H-formulation in 2-D, 2-D axisymmetric and 3-D. 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 permanent magnet 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-types bearings, the anisotropic homogenization proposed in [62] and extended in [83] is a good alternative. But it should be used with caution, and the first validations proposed in [53,84] 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 [85]. Finally, further vetting and refining of the models could help developing and improving lumped parameter SMB models [86,87], as a mean of drastically speeding up simulations.