Modeling and identification of the hysteresis nonlinear levitation force in HTS maglev systems

High-temperature superconducting (HTS) maglev has great potential in the field of high-speed transportation due to its capability for passive stabilization. The levitation force between the bulk HTSs and the permanent magnet guideway is a significant parameter relating to operational safety and comfort. This force has an obvious hysteresis nonlinear characteristic, which can be represented by nonlinear stiffness and damping. The stiffness and the damping are functions of vertical displacement and velocity, respectively. The vibration velocity of a HTS maglev vehicle can at times exceed 100 mm s−1, but the existing levitation force test methods are almost quasi-static. These methods are unable to accurately measure the damping characteristic of the maglev system. In this paper, a viscoelasticity model is introduced to describe the dynamic force. The parameters in the model are identified using the least square method based on the vibration response of the HTS maglev system. Meanwhile, the effectiveness of the model and identification method are tested by numerical simulations. The hysteresis loops derived from the motion theory coincide with the practical ones. Finally, the method is applied to identify the parameters of hysteresis nonlinear levitation force in a previous experiment with dampers. Based on the established hysteretic model, the dynamic characteristics of the HTS maglev system can be well presented.

The levitation force between the HTSs and the PMG is an important factor affecting the operation quality. The strong nonlinear E(J) power law of the HTS [21,22] makes the levitation force appear ot have an obvious hysteretic nonlinear characteristic, which suggests that the force is related not only to the displacement but also to the velocity [16,23,24]. In quasi-static measurements, the round-trip curves of levitation force cannot coincide with each other and a hysteresis loop is formed with the two curves [25][26][27][28]. In experiments, the bulk HTSs are usually driven near or away from the PMGs by electric motors. And the levitation force will be collected by a force sensor in real time.
In practical operation, the vibration velocity of a HTS maglev vehicle sometimes exceeds 100 mm s −1 . The experiments however, are always quasi-static and the relative velocity is less than 10 mm s −1 . Additionally, the range of displacement in quasi-static experiments is large, and can even reach 50 mm [7,28] in most cases. However, the displacement amplitude is usually no more than 5 mm in the practical vibration of HTS systems [14,15,18,20]. Hence, the existing quasi-static measurements are unable to show the damping characteristics (related to the vibration velocity) of the maglev system accurately.
Compared with experiments, simulations are more efficient because the parameters can be reset more easily. In recent years, finite element method (FEM) models have been widely used to investigate the electromagnetic behavior between HTSs and permanent magnets [23,24,[29][30][31][32][33][34][35][36]. Ma [23] and Huang [24] showed that with the increase of vibration velocity, the levitation force exhibits a remarkable increasing tendency. In their simulations, the maximum velocities are 25 mm s −1 and 10 mm s −1 , respectively. Liao [36] investigated the dynamic levitation force and its decay under varying external magnetic fields. These studies reflected the damping characteristic of the hysteresis nonlinear levitation force.
In FEM simulation, the process of the HTS maglev vibrating at a certain velocity is usually simulated. However, in practical operation, the displacement, velocity and acceleration of the vibration are time-varying. On the other hand, if the dynamic equations are coupled with the electromagnetic field equations to solve the vibration process, the computation will be greatly increased, resulting in low computational efficiency. Therefore, it is necessary to explore other methods with higher precision and efficiency for dynamics study.
Some mathematical models have been proposed for studying dynamics of the HTS maglev. Hikihara [37] introduced a state variable model to describe the hysteresis levitation force and numerically investigated the typical nonlinear dynamic issues like quasi-periodic, period doubling and chaotic motion of the HTS maglev system. Then, Zhuo [38] modified the state variable model with better velocity dependence. By combining a linear damping term with a nonlinear stiffness term, we proposed a modified levitation force model. Based on this model, we have analyzed the nonlinear vibration through both numerical simulations and the approximate analytic solution [39]. The outcomes matched well with the experimental results, which denoted that the levitation force can be characterized by a nonlinear skeleton function and a linear damping. The stiffness term and damping term can be described simultaneously by the aforementioned models in a large stroke. However, the parameters in these models should be determined by fitting the quasi-static experimental results or the FEM numerical results. So, these models have the same disadvantages as the quasi-static experiments and the FEM simulations in dynamics.
The operation of a maglev vehicle is always accompanied with the vibration, so the displacement, the velocity and the dynamic levitation force are all time-varying. Meanwhile, the range of displacement is always small (less than 10 mm), which means that the vibration process can be regarded as an elastic vibration. Hence, a reasonable mathematic model describing the dynamic levitation force is needed and would make the analyses more accurate. In a recent study, we incorporated an electromagnetic shunt damper (EMSD) into the HTS maglev system [40]. The experiments showed that the vibrations were reduced in varying levels with the employment of the EMSD. However the specific damping values under different conditions were not identified or analyzed theoretically since there was no reasonable mathematic model to characterize this process.
In this article, we introduced a nonlinear viscoelasticity model to describe the levitation force in a small stroke. Based on this model, we identified the parameters through dynamic response according to least square method. Firstly, the quasistatic levitation force was obtained through experiments and simulations. Secondly, the logic of parameter identification was also reported. Several dynamic simulations were executed to verify the effectiveness of the identification method. Thirdly, based on the mentioned experiments [40], the proposed model was used to identify the dynamic levitation force and its damping term.

Experiments and simulations of levitation force
In this paper, the small and large range of movement are defined by displacements less than 10 mm and greater than 10 mm, respectively. In the large range case, the levitation force between HTSs and PMG can be measured by experiments. Due to the limitation of equipment precision, it is difficult to measure the force in a small range case, as shown in figure 2(c). In addition, when the motion velocity gets faster, obvious errors will occur because of the vibration of the motor, as shown in figure 2(a). Hence, the simulations were used in the small range case to complement the measurements.

Experimental setup and procedure
We used the SCML-01 device [41] to measure the levitation force. Figure 1 shows the schematic diagram of the measurement device-SCML-01. This device includes a liquid nitrogen vessel, a PMG, a data collection/processing system and a mechanical drive/control system. In experiments, four 64 mm×32 mm×13 mm multi-seeded rectangular YBa-CuO bulk HTSs, made by ATZ GmbH, Germany, were fixed in the vessel and placed above the PMG. This kind of bulk has superior superconducting performance across the grains [42] and was used in the ring test line [7]. The PMG used in the experiments has the same interface and material as that used in the ring test line.
The YBCO bulk HTSs were cooled at a height of 30 mm with liquid nitrogen [43] for 15 min. After that, these bulks were moved vertically between 60 mm and 10 mm above the PMG driven by step motors at different velocities (1 mm s −1 , 2 mm s −1 , 5 mm s −1 and 8 mm s −1 , respectively).

Simulation process
The simulations are based on the vibration experiments in which the EMSD was employed to suppress the vibration for the HTS-PMG system. In the vibration experiments, the levitation body weighing 4.72 kg vibrates in the vertical direction at a height of around 18 mm. For a detailed description of the experimental setup and process refer to [40]. The HTSs and PMG are the same as those in the aforementioned experiments in section 2.1.
The field-cooling height (FCH) was also set at 30 mm in the simulations. In one simulation case, the velocities of movement were set as 1 mm s −1 , 2 mm s −1 , 5 mm s −1 and 8 mm s −1 respectively, which were the same as those in section 2.1, but the range of displacement was set as 8 mm. In another case, the velocity of the relative motion was selected as 1 mm s −1 and the ranges of displacement were set as 1 mm, 2 mm, 4 mm and 8 mm respectively around the height of 18 mm above the PMG (static levitation height in the EMSD experiment). The calculation process is as follows.
The model of the PMG was established in the FEM software COMSOL Multiphysics 5.3a with AC/DC module selected. The remanence B r of the PMG was set as 0.91 T. And the magnetization direction of the PMG is distributed according to the Halbach array [44]. The electromagnetic properties were simulated using the H-formulation. The governing equations were derived from Maxwell's equations, known as Faraday's and Ampere's laws.
where H is the magnetic field intensity; E is the electric field intensity; and J is the current density. The E-J relationship is considered as [45].
where E 0 is the critical current criterion; J c is the critical current density; and n is the power law exponent. The parameters in the simulations are shown in table 1.
The equations above are solved in the magnetic field formulation (MFH) module. The levitation force can be obtained by equation (4).
where B x is the external magnetic flux density in the xdirection; J y is the current density in the y-direction; S is the cross section of the HTSs. Figure 2 shows the hysteresis loops of the levitation force. , it can be found that the levitation force is slightly larger when the velocity is bigger. This is because the magnetic force is determined by the external field and the shielding current induced in the bulk HTSs during the motion. As the moving velocity increases, the external field applied to the bulk HTSs changes faster, which increases the induced shielding current, and the magnetic force is also enhanced. This phenomenon will act as a damping term in the dynamics.

Results and discussion
In figure 3, we fixed the relative moving velocity to 1 mm s −1 . Then we simulated the hysteresis loops of the levitation force with the ranges of displacement as 1 mm, 2 mm, 4 mm and 8 mm to study their dependence on the relative displacements. In the vibration test [40], the static levitation height is 18 mm, and the gravity of the suspension body is 46.3 N. Therefore, the point (18 mm, 46.3 N) was set as the coordinate origin in figure 3. We got four hysteresis loops, and the loops look like a single curve because of the low velocity and small displacements. By regarding them as a single curve, it is not hard to find that the slope of the curve is negative and it is an obvious concave curve whose graph is always below its tangents. A quadratic polynomial, F Lev =−9.07z+0.65z 2 was used as a skeleton function to fit the loops and they matched well. Where F Lev is the relative levitation force which is defined as the actual levitation force minus the gravity and z is the vertical displacement. The hysteresis loops can be ignored in figures 2(b) and 3, because the maximum velocity is 8 mm s −1 and the maximum displacement is 8 mm. However, it is worth noting that the hysteresis loops could not be ignored in practical operation because the moving velocity may exceed 100 mm s −1 . As discussed in the introduction, a linear damping term should be added into the quadratic polynomial. Hence, the dynamic levitation force model is written as follows: where z is the motion velocity. The levitation force can be simplified as a viscoelastic mechanical model shown in figure 4. Equation (5) is a classical viscoelasticity model, which is called the Kelvin model and is often used to describe the hysteresis. As shown in equation (6), F Lev consists of a nonlinear elastic force F 1 and a linear damping force F 2 .
q q q q

Identification of the hysteresis nonlinear levitation force
The dynamic response of a structure highly depends on the dissipation ability of energy by means of hysteretic behavior. Hence, parameters in the model can be identified under the given dynamic excitation, dynamic response and a certain model [46]. In this part, we will identify the parameters in equation (5) from the dynamic response using the least square method.

Nonlinear least square method
Supposing that we have n observations (w i , u j ), i=1, 2, K, n, from a fixed-regressor nonlinear model with a known functional relationship f. Thus, where z i is random noise vector whose mean is 0; θ is an undetermined parameter vector with p dimensions; the leastsquares estimate of q* (the true value), denoted by q minimizes the error sum of square Assuming that q is the result of identification, when each q f u , i ( ) is differentiable with respect to θ, q will satisfy For most nonlinear models, equation (9) cannot be solved analytically, so iterative methods are necessary. In this paper, the Newton method [47] is used, in which q S ( ) is expanded using a quadratic Taylor expansion, let where q  is the transpose vector of θ. Assuming that q m is in the neighborhood of q . * Then we have the quadratic approximation q q q q q q q q q q The right-hand side is minimized with respect to θ when This suggests that, given a current approximation q , m the next approximation should be After several steps of iteration, q can be obtained.

Numerical verification of parameter identification based on least square method
For a single degree of freedom system, according to Newton's second law, its motion function can be written as where m is the mass of the levitation module and f (t) is the external excitation force. In this part, Function (15) will be solved by the Runge-Kutta method based on MATLAB.
Then the method discussed in section 3.1 can used to identify the parameters in equation (5). In the simulations, m equals 4.72 kg. The calculation time is set to 2 s and the time step of integration is 10 −4 s.

Sine excitation.
Assuming that the excitation displacement is z 0 =1×10 −3 sin(5t). The high-speed HTS maglev trains generally run above elevated bridges, and the excitation caused by the bridge deformation can be regarded as a sine curve. The excitation and response are shown in figure 5. In engineering, the measured data is often polluted by noise. To get closer to the reality, we added 5% white noise to reduce the accuracy of the data. Without loss of generality, the set values of q 1 and q 2 are randomly generated, which are close to the corresponding values in the skeleton function shown in figure 3. The set value of q 3 is a random value based on the experimental experiences. The set values and the identification values are shown in table 2. It can be found that the parameters are identified with small errors in the table. The accuracy of q 2 is slightly lower than the other two because its value is much larger than the others, causing the reduction of the accuracy in iteration.
The levitation force was rebuilt with the identified results and compared with the actual values, shown in figure 6. The actual values were obtained by solving equation (15) with set values and the rebuilt values were obtained by solving equation (15) with identification values.  The shape of the hysteresis loops shown in figure 6(a) are similar to those simulated from FEM shown in figure 3, while the hysteresis loops in figure 6(a) are more obvious. That is because the velocity in the vibration is much larger. Specifically, the maximum relative velocity (the velocity of levitation device minus the velocity of excitation) reaches 96.7 mm s −1 , as shown in figure 6(b). At the same time, in figure 6(a), the value of maximum active displacement (about 2.7 mm) is larger than the negative one (about −2.5 mm) even if under the symmetric excitation shown in figure 5(a). The reason is that the stiffness becomes weak as the HTSs moving away from the PMG, and it can be seen in figure 3, a typical phenomenon in nonlinear vibration.
3.2.2. Random excitation. The HTS maglev system is mainly composed of onboard bulk HTSs and the PMG equipped with NdFeB permanent magnets and yoke. Moreover, the assembly tolerances such as surface peeling and joint gap will lead to the unevenness on the guideway surface, which can be manifested as the irregularity in the maglev system. The irregularity acts as an excitation to the vehicle system. The excitation is random like the wheel on rail transit systems. Thus, a random excitation is chosen in this section.
The excitation and response are shown in figure 7 and a 5% white noise is also added. The identification values are shown in table 3. In figure 7, the vibration is more severe and irregular. The maximum amplitude of the velocity exceeds 140 mm s −1 , as shown in figure 7(c), which means the damping force in this condition reached 5.6 N (0.14 m s −1 × 40.35 Ns m −1 =5.6 N), while this phenomenon does not appear in the quasi-static experiments and FEM simulations.
The set values in table 3 are generated using the same approach as in table 2. In order to verify the wide applicability of the identification method, the different set values were employed in sections 3.2.2 and 3.2.1. Table 3 shows that the parameters are identified with small errors. The accuracy of q 2 is also slightly lower than the other two in this condition. Under random excitation, the levitation force was rebuilt with the identified results and compared with the actual values, shown in figure 8.
In figure 8, the motion becomes irregular instead of periodic or quasi-periodic. The change in levitation force also becomes more dramatic, but the rebuilt results match well with the actual values. The applicability of the identification method was further proved.

Application of the parameter identification in vibration test
To enhance the damping effects, we incorporated an EMSD into the HTS maglev system [40]. The EMSD, which contains coils, capacitors and resistors, can convert some of the vibrational kinetic energy into electrical energy and dissipate it with resistors in the form of Joule heat. In this section, a set of experimental results were chosen to verify the parameter identification method proposed in section 3.1.
In the case of resonance, the frequency of external excitation was 8 Hz. The vibration response is shown in figure 9. The FCH in this test was 30 mm. In figure 9, it can be seen that within the scope of the experiments, the damping effect can be improved with decreasing of the resistance R.
Using the method of parameter identification discussed above, the parameters of dynamic levitation forces with   In the test, the acceleration was collected by a B&K vibration analyzer. To reduce the experimental error, we started to collect data after the vibration became stable. The duration of data collection was 8 s, and a total of 4096 points were collected. After that, the vibration with no damper was rebuilt using the identified parameters. The test acceleration of the vibration actuator was used as the excitation. The results are shown in figure 10 and they are compared with the EMSD vibration experimental data and they match well with each other. In this way, the rationality of the dynamic levitation force model and identified method are further demonstrated.

Conclusion
To better understand the dynamic behavior, a viscoelasticity model was introduced to describe the dynamic levitation force of the HTS maglev system. At the same time, the nonlinear least square method was successfully used to identify the parameters in the viscoelasticity model. It has been proved that the proposed method can be applied to test the dynamic levitation force of the HTS maglev system with good applicability and rationality. Based on the studies mentioned above, conclusions can be drawn as follows: 1. The dynamic levitation force of the HTS maglev system can be modeled by a nonlinear stiffness term and a linear damping term. 2. The viscoelasticity model was applicable to describe the damping term in the dynamic levitation force of the HTS maglev system.