Field–circuit model of the radial active magnetic bearing system

Paper presents a mathematical model of the radial active magnetic bearing, which is indented for the simulation of the magnetic bearing dynamic response. The circuit model of the bearing is based on differential equations. The circuit model has incorporated results of magnetic field analysis, which led to the creation of the field–circuit model. Presented model of the magnetic bearing takes into account the necessary control system. The experimental results are presented to validate the proposed model.


Introduction
Magnetic bearings (MBs) represent an alternative support of the rotor in comparison with traditional bearings, i.e., ball or journal ones. MBs have found applications in many industrial devices, for example, in high-speed turbines, energy storage flywheels, turbomolecular pumps, turbogenerators, machine tool spindles and compressors [1][2][3]. The benefits of using magnetic bearings are well known [1]. Owing to the contactless operation of the rotor, the bearing provides a lack of friction, absence of lubricating substance, good vibration damping, online monitoring of the operation and reduced maintenance and operation costs.
Magnetic suspension dedicated to electric machine usually consists of two radial and one axial electromagnetic actuators and a control system. The actuator of the radial active magnetic bearing (RAMB) comprises two elements-a stator and rotor. The interaction between the stator and rotor is based on the principle of the electromagnetic interaction. The current flowing in the windings causes the pull of the movable ferromagnetic material. Unfortunately, the stable levitation of the RAMB rotor is only achievable by using position controllers. In this paper, a field-circuit model of the RAMB system dedicated to the simulation of the transient states is described. The model is based on a set of the differential equations implemented in MATLAB/Simulink software. The main parameters of the RAMB were obtained from the magnetic field analysis. The model also includes the necessary control system with PID controllers for the rotor position and PI controllers for currents excited in windings. The presented simulation model was compared with the real object.
The aim of this paper is to present an effective and fast model of the RAMB system, which can be used to test various controllers as well as determine its parameters.

Structure of the active magnetic bearing
The construction of the RAMB actuator consists of a stator and rotor. In order to significantly reduce the eddy currents effects, the stator and rotor were fabricated from 0.5-mmthick silicon steel M600-50A. The cross section of the RAMB actuator is presented in Fig. 1. The winding of the stator has wound twelve coils, which allows creating several configurations of the magnetic field excitation. It is possible to obtain three sections, four sections and six sections of the winding. In this paper, we consider four sections of the winding, which permits to implement the differential driving mode [1]. For such mode, currents excited in windings are equal to:

c)
where I b denotes the bias current and i cy and i cx indicate the control currents for y-and x-axis, respectively. The bias current is used for linearization of the magnetic force characteristic within the operating range of the magnetic bearing [1]. The main parameters of the RAMB actuator are presented in Table 1.

Mathematical model of the active magnetic bearing
Magnetic field distribution was obtained from the 2D finite element method (FEM). The field analysis is based on the calculation of the magnetic vector potential A [4]. Although The geometry of 2D FEM model and a quarter of the mesh the 2D FEM model neglects the end effects of the magnetic field, they have a minor impact on integral parameters of the magnetic field [5]. In numerical model, the eddy currents effects were neglected and the magnetic field was considered to be stationary. However, the impact of the manufacturing process was taken into account in the model, which changes the magnetic properties of narrow layers of the stator and rotor sheets [6,7]. Hence, the air gap in the simulation model was increased in comparison with the real object. The calculation area includes the whole geometry of the RAMB actuator cross section. Dirichlet's boundary conditions were assumed at its edges (Fig. 2).
The finite element mesh was carefully selected, in order to obtain the balance between accurate results and short time of computation. The stator and rotor subregions were discretized with the fine mesh. The discretization of the air gap has a significant impact on the force calculation results. Therefore, the air gap was divided into two subregions.
Based on the FEM algorithm, two cases of the magnetic field distribution in the RAMB actuator were calculated. In the first case, the central position of the rotor and the bias current I b excitation in windings were assumed. For such assumption, the magnetic forces F x and F y are equal to zero (Fig. 3). Magnetic field density in each electromagnet is approximately 0.74 T. That value constitutes half of the maximal field density (B max 1.5 T) assumed for the RAMB.
In the second case, the central position of the rotor, the control current i cy equal to 5 A was assumed. Thus, the current excited in the first winding is equal to 10 A, whereas in the third winding it is equal to zero. In this case, the magnetic force F y in the y direction is equal to 68.9 N. This value rep-  Magnetic field density in the first electromagnet is equal to 1.42 T (Fig. 4).
The calculated magnetic flux distribution was used to estimate the magnetic force generated by all four electromagnets. The magnetic force was calculated from Maxwell stress tensor. The calculations were executed for various cases, over the entire operating range of the magnetic bearing currents I 1 , I 2 , I 3 , I 4 ∈ (0, 10 A) and the various positions of the rotor shaft x, y ∈ (− 400, 400 µm).
The calculated values of the magnetic force generated by the first electromagnet, as a function of the winding current I 1 and the shaft position y, are given in Fig. 5. The flux linkage was calculated from the following equation [8]: where N is the turn number of the stator winding and φ k is the magnetic flux linking with one turn of the winding. The flux linkage of the first electromagnet as a function of the winding current I 1 and the shaft position y is given in Fig. 6.
The magnetic flux linked with the windings was used for calculating the dynamic inductances and the velocityinduced voltage coefficient [9]. The dynamic inductance L d1 was calculated from the following expression: whereas the velocity-induced voltage coefficient h v1 was calculated from the expression: Figure 7 presents values of the dynamic inductance L d1 in function of the rotor position and the current intensity. It is noticeable that the value of the dynamic inductance increases with increasing the coordinate y. This effect is due to decreasing value of the air gap. The exception is the range of high values of the current intensity and high values of the rotor deflection when the values of the dynamic inductance slightly decrease. The reason of that is saturation of the magnetic material. Figure 8 presents the velocity-induced voltage coefficient in function of the rotor position and the current I 1 . The value of the velocity-induced voltage coefficient increases in accordance with the current value and position of the rotor in the y-axis. Also, there is a region, where the value of the velocityinduced voltage coefficient is decreasing. The reason for this is the saturation of the magnetic material, because of the proximity of the stator and rotor.
The main parameters of MBs are current k i and position k s stiffness [10]. These parameters are often included in the equations of motion. The parameters are calculated from the following expressions:  The mathematical model of the RAMB actuator is described by a set of the ordinary differential equations: The voltage and current balance Eqs. (6.a, 6.b, 6.c) govern the electrical characteristics of the RAMB actuator. The parameter R i (where, i 1, 2, 3, 4) denotes the winding resistance. Parameters L di and h vi (where, i 1, 2, 3, 4) indicate the dynamic inductance of the windings and velocity-induced voltage coefficients, respectively. Equations (7.a, 7.b) govern the rotor motion and include the magnetic forces, F i (where, i 1, 2, 3, 4) generated by four electromagnets. Also, these equations incorporate static unbalance of the rotor with the eccentricity e S , which denotes a difference between the rotor center and the center of the rotor mass. Symbols m and ω indicate the mass and the angular velocity of the rotor, respectively. The equations of the motion include the gravity force acting on the rotor, which is equally divided between both axes. Equations (6.a, 6.b, 6.c) and (7.a, 7.b) were implemented in Simulink/MATLAB software. The values of constant parameters of the mathematical model are listed in Table 2. The winding resistances and mass of the rotor were measured and constitute input data. The value of the eccentricity e S was assumed in order to achieve a similar, to the physical object, dynamic response of the rotor during rotation.
The other parameters, like dynamic inductances, velocityinduced voltage coefficients and magnetic forces are nonlinear functions of the currents excited in windings and position of the rotor. Therefore, these parameters were obtained from previously presented finite element model and were incorporated in the model as lookup tables. Figure 11 presents the simulation model of the first electromagnet that corresponds to Eq. (6.a).
The RAMB is an inherently unstable device that means without the control system stable levitation of the rotor is impossible. Therefore, in order to obtain dynamic responses, the control system for considered MB was designed. It consists of four switching power amplifiers, four PI current controllers and two PID controllers of the rotor position. Figure 12 presented the layout of the RAMB mathematical model with the control system.
The pole placement method was used to obtain values of parameters for the current and position controllers [11].
The coefficients values for the position controllers PID are given in Table 3, whereas the coefficients for the current controllers PI are given in Table 4. The value of the control current, for both axes in a steady state, is equal to 1.3 A. For this state, the electromagnets generate only the forces to balance the weight of the rotor. The settling time t s of the MB system equals 36.8 ms.      One can find a good agreement between magnetic flux density values obtained from the field simulation and from measurements. Differences between the calculated and measured results do not exceed a few percents.

Simulation results
For the central position of the rotor, the magnetic force values, in function of the control currents i cx and i cy , are nearly linear, which are presented in Figs. 16 and 17. However, a difference between the calculated and measured values of the magnetic force is noticeable. The largest differences occur for the maximal control current values, and they are equal to 5.0% for x-axis and 8.8% for the y-axis.
In Fig. 18, an outline of the test stand which contains the RAMB system is presented. The RAMB actuator is supplied by the switching power amplifier with supplying voltage U dc 35 V and the switching frequency f PWM of the PWM system equal to 20 kHz. All control tasks are carried out by the magnetic bearing controller (MBC), which was based on the ARM7 microcontroller. An additional motor supplied with the variable frequency drive was employed for the rotation of the RAMB rotor. An analog input card was used for measuring of all signals. Figure 19 presents a photograph of the test stand. Figure 20a, b presents time responses of measured currents excited in electromagnets and the rotor displacement in xand y-axis for the rotor lifting.
The values of control currents i cy , i cx in a steady state are approximately equal to 1.5 A, and they are slightly higher for the real object in comparison with those obtained from the computer simulation. This difference could be the result of omitting the eddy currents in the finite element model. The eddy currents are responsible for decreasing the values of the current and position stiffness. The settling time of the MB system is equal to 28.8 ms and its value is lower than in the computer simulation.
In Fig. 21, measurement results during the rotor rotation are presented. The maximal deviation of the rotor is equal to an electric motor, (3) the radial magnetic bearing, (4) switching power amplifier and magnetic bearing controller and (5) a measurement card 73.9 µm. In contrast to the simulation, the rotor center path is not a circle which may indicate an asymmetry of the force generation. Table 5 consists of gathered values of the RAMB system performance indicators. The parameters t Sx and t Sy denote the rotor settling time for the rotor lifting, respectively, in x and y direction. Parameters J 1x and J 1y denote the integral of the square error e(t) calculated during rotor lifting: where t k denotes the final time of the system response observation.
Parameter J 2 denotes maximal deviation of the rotor center.
One can notice that the settling time and the integral of the square error took higher values for the simulation model than for the real object. It could be due to the higher value of damping in the physical model than in the numerical simulation.