Simulation of a Synchronous Planar Magnetically Levitated Motion System Based on a Real-Time Analytical Force Model

Abstract

The existing simulation method for the control of linear or planar magnetically levitated actuators always ignores the characteristics of the real physical object, which deteriorates the accuracy of the simulation. In this work, the proposed emulator for the magnetically levitated actuator is developed to consider both the force characteristics and the control algorithm. To model the real controlled object, the mathematical model for 1D (one-dimensional) and 2D (two-dimensional) magnetic arrays is derived where the yaw angle is taken into consideration using the coordinate transformation. The solution of the mathematical model is compared with the commercial BEM (boundary element method) software and the measurements from a force and torque testing setup to highlight the accuracy of the proposed mathematical model. Compared with the traditional simulation method of the motion control systems founded on the simplified system transfer function, the proposed simulation method has higher consistency and is closer to reality. The accuracy and efficiency of the proposed magnetic force model are further verified by the emulator based on the numerical force model and the testing data of the real setup.

1. Introduction

6-DoF magnetically levitated actuators (MLAs) have been increasingly used in micro-assembly, micro-fabrication, and semiconductor manufacturing [1,2,3,4] because the moving part is capable of undertaking the large motion range in each direction benefiting from the characteristics of non-connect. The performance of the MLA system is mainly determined by the physical properties [5,6] and the performance of the control system [7,8,9,10]. The physical properties containing the maximum output magnetic force and motion range are decided by the structure, and the control strategy affects the dynamical performance such as the frequency bandwidth, the resolution, and the responding speed. Hence, a large number of researchers have focused their research on the design of the physical structure and the control algorithm for the MLA system, and plenty of research efforts have been published. Before these efforts are employed in the MLA system, a series of the test is needed to verify the design. However, these designs cannot be tested straight in the real plant. Before the MLA is fabricated, the integration of the controller and managerially levitated physical system may cause unpredictable problems in the process of testing. These uncertain issues may cause damage to the MLA whose manufacturing process is complicated and expensive. In this situation, simulation for the MLA system is proposed. Instead of the real plant, in the simulation, the control strategy is evaluated on the virtual object, which could decrease the verification cost and unknown risks. Through analyzing the simulation results, the control algorithm can be improved and the control parameters can be adjusted to get a better control effect. Therefore, simulation plays an important role in the design and testing of the MLA system. However, in the traditional simulation methods of the MLA, the virtual controlled object is just viewed as a second-order system where the desired force and torque produced by the control strategy are regarded as the driving force and torque acting on the MLA. In the actual control system of the MLA, the required magnetic force and torque are the output of the control algorithm, and the driving magnetic force and torque are generated by the coils of the input current. This difference causing the error between the real control system and the simulation model will affect the reliability of the final simulation results. It is ignored in the traditional simulation. To avoid this case, the co-simulation model reflecting not only the performance of the control strategy, but also the influence of the physical structure is proposed.

To establish the co-simulation model of the MLA, the physical process should be well understood and generalized by the mathematical equations that describe all major features of the real MLA system. These mathematical equations are utilized to create the simulation model of the real system. As a consequence, the co-simulation model with a real MLA control strategy and the emulated plant model is built to test and verify the MLA system in real time. This kind of simulation model has been proposed, and some efforts have been published. K. B. Wipke utilized the look-up table to model the vehicle integrals in the co-simulation model for an automotive system [11]. It is time consuming to get the accurate static table for vehicle components. M. Cipek used the bond graph methodology to build a co-simulation model for a two mode power-split HEV [12]. E. Vinot built the co-simulation model for hybrid drivetrains and integrals in an automotive hybrid system [13]. S. C. Carpiuc came up with the co-simulation model for a permanent magnet synchronous machine (PMSM) [14]. Analyzing the existing co-simulation models, they concentrate on the components of the automotive vehicle, hybrid system, and electrical machine, and none of them focuses on the MLA. However, this idea can be referred to in the simulation model of the magnetically levitated positioning system.

As mentioned before, to successfully build the real-time co-simulation model, the mathematical equations for the MLA are essential. For the MLA, the physical process is described by the magnetic force model, and plenty of studies have been reported [15,16,17]. Finite element analysis (FEA), the surface charge model, the harmonic model, and the analytical model are widely used magnetic force models [18,19,20]. These models can accurately calculate the magnetic force and torque in a numerical or analytical approach. However, they cannot meet the requirement of real-time simulation because of the large computational burden, causing a big time delay, which decreases the simulation efficiency. The analytical model derived from the harmonic model is suitable for physical model simulation. The shortcoming of this kind of model is that it is hard to get accurate results unless the yaw angle is $0^{\circ }$, ${45}^{\circ }$, and $-{45}^{\circ }$. However, in the process of controlling the actuator, it is inevitable that the magnetic array will un-align the coil, which will cause the yaw angle between the magnetic array and the coil. We can conclude that the yaw angle of the mover complicates the magnetic force and torque expressions, and it is hard to reduce the computational complexity considering the yaw angle, which will decrease the simulation efficiency. J. Peng proposed a real-time analytical model, which considers the corner of the coils, for real-time control [21]. This wrench model has good computational accuracy and lower computational expenditure, but it does not take the rotation angle into account. V. H. Nguyen put forward a real-time wrench model derived from the harmonic model [22]. This physical model can realize rotation, but the rotational range is limited. X. Lu designed a magnetically levitated rotatory table [23]. The proposed wrench model can calculate the magnetic force (and torque) in real-time and realize unlimited rotation around the axis, but the stroke in the radial direction is limited. Similar to [23], X. Lu utilized a semi-analytical model to model the planar magnetically levitated actuator in which the relative translation of the mover and stator is not taken into consideration [24]. F. Liu proposed a semi-analytical model based on the circle coils for real-time control [25,26]. This model uses the symmetry of the circle coils to view the yaw angle as a variable in the calculation of the magnetic force (and torque). This wrench model can achieve large translation and rotation, but the circle coil does not have universal applicability compared with the racetrack coil. In summary, the universal real-time magnetic force and torque, which has a large stroke in the translation and rotational axes for 6-DoF MLAs, have not been studied in the existing magnetic force and torque models.

This paper presents a simulation model focusing on the MLA system. In this aspect, a new real-time co-simulation model of the physical model and controller is proposed for the 6-DoF MLA system, which is developed in MATLAB/Simulink and can be used for desktop simulation. This simulation model consists of the magnetic force model for the MLA and the tested controller. The desired force and torque calculated by the tested controller are transformed into the input current value of the magnetic force model of the MLA via the transformation matrix. The driving force is the output of the magnetic force model. The other contribution of this paper is the universal real-time mathematical magnetic force model for the MLA, which can be used in either one-dimensional or two-dimensional magnetic fields.

Section 2 shows the real-time co-simulation model for the MLA. The mathematical model is briefly introduced in Section 3. The experiment and validation of the simulation model and mathematical model are shown in Section 4, and the conclusion is presented in Section 5.

2. Design and Implementation of the Simulation Model

As one of the two achievements of this paper, the real-time co-simulation model is designed and implemented with MATLAB/Simulink. The structure of this simulation model is shown in Figure 1. The MLA model block is the greatest difference between the traditional simulation method and the proposed co-simulation model for the MLA. In the traditional simulation model, the MLA model is just viewed as a second-order system, and the influence of the physical properties on the control algorithm is ignored, which influences the simulation result. To avoid this case, the traditional second-order model should be replaced by the MLA physical model. In other words, the tested control algorithm and the physical model of the MLA are combined together in the co-simulation model. In the co-simulation model, the input of the physical model is the current value, and the output is the magnetic force.

This simulation model consists of four main blocks, the Tested Controller Block, the Force-current Block, the MLA Model Block, and the Integrator Block. The specific control algorithm is implemented in the Tested Controller Block, and the output of this block is the requiring driven force and torque for the MLA calculated by analyzing the desired and feedback position information. As mentioned before, the input of the MLA model is the current value. It is necessary to transform the desired force and torque into the current value. The Force-current Converter can achieve this goal. In this block, the input is the desired force, torque, and feedback position information of the MLA. The current value in each coil is calculated by the relationship between the force, the current value, and the relative position of the coil and magnetic array, obtained by curve fitting of the measurement data in this simulation model. The MLA Model Block implements the magnetic force model for the MLA where the real driven force and torque on the magnetic levitation are calculated in this block. Through the Integrator Block, the real-time position of the MLA is calculated based on the magnetic force and torque and fed back to the tested controller, force-current converter, and mathematical model for the MLA.

For real-time feedback of the position information, the driven magnetic force and torque should be accurately calculated by the magnetic force model of the MLA in real time. To realize the real-time calculation of the magnetic force, the mathematical model for the MLA is proposed in this paper. A brief introduction to the proposed mathematical model is given in Section 3.

3. Mathematical Model

3.1. Structure of the Maglev Actuator

There are two mainly kinds of structure in existing MLAs: one is based on a 1D Halbach array, and the other is based on a 2D magnetic array. Their typical design is respectively shown in Figure 2a,b. The MLA based on 1D Halbach arrays is shown in Figure 2a. This MLA is made up of two parts: the stator and translator. The stator consists of eight racetrack coils, and the translator consists of four 1D Halbach array. This design has the advantages of a simple structure and strong representativeness where the positional relationship of each part is relatively simple, and the structure of each part is universal such as 1D Halbach array and racetrack coil. The force and torque generated by each coil are controlled by the value of the current in the coil, and these eight coils are controlled separately. The total force and torque on the translator is the vector sum of the force and torque generated by each coil on the four 1D Halbach arrays. Therefore, the motion of the actuator is controlled by these eight racetrack coils, and each coil is controlled separately. Therefore, the actuator can be controlled to achieve a specific motion including translation along the axes and rotation around the axes by applying specific currents to the eight coils. Figure 2b shows the MLA based on a 2D magnetic array. Different from the 1D Halbach array, in this structure, the translator consists of only one 2D magnetic array, and the resultant force and torque are the vector sum of the force and torque generated by the eight coils on this 2D magnetic array.

In the prototype of the MLA based on the 2D magnetic array, the 2D magnetic array utilizes the 2D superimposed magnetic field proposed by V. H. Nguyen [22,27]. This superimposed magnetic field can be viewed as the superimposing of two single-axis Halbach array shown in Figure 3b.

To facilitate the analysis of the force and torque acting on the translator, every two coils placed side-by-side and one 1D Halbach array upon them form a linear magnetic actuator unit, which is shown in Figure 3a, and this MLA consists of four units. Since the total force and torque on the translator are the sum of each unit and every actuator unit is controlled separately, this mathematical model can be derived from only one unit. In each actuator unit, the force and torque on the magnetic array are also the vector sum of the force and torque generated by the two side-by-side coils, and the current value in a coil has no influence on the other. As for the 2D magnetic array, it also has a periodic arrangement, and the total magnetic force on the 2D magnetic array is also the vector sum of the eight coils. Consequently, the mathematical model is finally derived from one racetrack coil only.

With the translation and rotation of the translator, the relative position of the coil and magnetic array changes in each actuator unit. The force and torque acting on the magnetic array are different when the relative position relationship of the coil and magnetic array changes. To calculate the force and torque produced by the coil, two coordinate systems are created, as shown in Figure 3a. The coil coordinate system’s origin is the geometry center of the coil. It is denoted by the superscript $\left\{c\right\}$. Because the coil is set on the stator, this coil coordinate system can be viewed as the global coordinate system. The magnetic coordinate system utilizing the geometry center of the magnetic array as the coordinate origin point is denoted by the superscript $\left\{m\right\}$. This coordinate system can be viewed as a local coordinate system. $L_c$, $W_c$, $H_c$, $R_c$, $L_m$, $W_m$, and $H_m$, which denote the design sizes of the coils and Halbach arrays, are depicted in Figure 3a. $L_c$, $R_c$, and $H_c$ represent the length, witch, and height of the long straight side of a coil, respectively. $W_c$ is the length of the short straight side of a coil. $L_m$, $W_m$, and $H_m$ are the length, width, and height of a magnet. In this prototype, the long straight side of the coil is longer than the length of the magnet, which can decrease the calculation complexity of force and torque.

3.2. Magnetic Flux Density

The magnetic flux density distribution can be described by many approaches [16,17], each of which has its shortcomings and advantages in different situations. In this paper, the harmonic model is derived to decrease the computational complexity. The magnetic flux density of the 1D magnetic array ${}^m{\mathbf{B}}_1$ at the point ${}^m\overrightarrow{x}={{[}^m{x}_0{,}^m{y}_0{,}^m{z}_0]}^T$ can be expressed as:

$${}^m{\mathbf{B}}_1=\sum _{n=1}^{1,5,9,\dots }·\frac{2B_0}{n\pi }·(e^{\tau m_t}-e^{\tau m_b})·e^{{\tau }^m{z}_0}·\left[\begin{array}{c}-\sin \left(\tau {}^m{x}_0\right)+\cos \left(\tau {}^m{x}_0\right)\\ 0\\ +\cos \left(\tau {}^m{x}_0\right)+\sin \left(\tau {}^m{x}_0\right)\end{array}\right],$$

(1)

where:

$$\tau =\frac{n\pi }{2W_m}$$

(2)

$m_t$ and $m_b$ are the top and the bottom of the 1D or 2D magnetic array. $B_0$ is the remaining magnetization of the permanent magnet. n represents the harmonic numbers for the ${}^{m}x$-direction. The magnetic flux density is solved by the method of the separation of variables, as the magnetic field is the gradient of the scalar magnetic potential. Considering the boundary conditions provided by the space periodic magnetic array, the 3rd, 7th, 11th...harmonic components are zero, and the even harmonic numbers are canceled by the trigonometric transformation. In this method, the magnetic flux density distribution of the 1D Halbach array is just along the x-direction and z-direction, and the y-direction end effect of the 1D Halbach array is neglected [28].

The 2D superimposed field is the vector sum of the single x-axis Halbach array and the single y-axis Halbach array, so the distribution of the 2D superimposed field ${}^m{\mathbf{B}}_2$ can be obtained:

$${}^m{\mathbf{B}}_2=\sum _{n=1}^{1,5,9,\dots }·\frac{2B_0}{n\pi }·(e^{\tau m_t}-e^{\tau m_b})·e^{{\tau }^m{z}_0}·\left[\begin{array}{c}-\sin \left(\tau {}^m{x}_0\right)+\cos \left(\tau {}^m{x}_0\right)\\ -\sin \left(\tau {}^m{y}_0\right)+\cos \left(\tau {}^m{y}_0\right)\\ +\cos \left(\tau {}^m{x}_0\right)+\sin \left(\tau {}^m{x}_0\right)+\cos \left(\tau {}^m{y}_0\right)+\sin \left(\tau {}^m{y}_0\right)\end{array}\right].$$

(3)

3.3. Force and Torque

3.3.1. Coordinate Transformation

Within the obtained magnetic flux density composition, the force and the torque between the magnets and the coils can be calculated by the Lorenz integral. However, the obtained magnetic flux density composition is built in the magnetic (local) coordinate system, and the force and the torque should be solved in the coil (global) coordinate system. Therefore, the coordinate transformation is employed to build the connection.

The schematic diagram of the actuator after rotating in the vertical direction is shown in Figure 4. The position of a random point in the coil is described as ${}^c\overrightarrow{x}={{[}^c{x}_0{,}^c{y}_0{,}^c{z}_0]}^T$, and the center points of the magnetic arrays shown in Figure 4a,b can be expressed in the global coordinate system as:

$${}^c\overrightarrow{p}={{[}^c{x}_p{,}^c{y}_p{,}^c{z}_p]}^T.$$

(4)

According to the coordinate transformation matrix, substituting (1) and (3), the magnetic flux density distribution along the x-axis ${}^c{\mathbf{B}}_x$ can be obtained:

$$\begin{array}{cc}\hfill {}^c{\mathbf{B}}_{xx}=& \cos \gamma ·\frac{\sqrt{2}{B}_0}{n\pi }·(e^{\frac{H_m\tau }{2}}-e^{\frac{-H_m\tau }{2}})·e^{\tau {(}^c{z}_0{-}^c{z}_p)}·\sin [\tau ·\cos \gamma ·{(}^c{x}_0{-}^c{x}_p)\hfill \\ \hfill & +\tau ·\sin \gamma ·{(}^c{y}_0{-}^c{y}_p)+\frac{3}{4}\pi ]\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {}^c{\mathbf{B}}_{xy}=& \sin \gamma ·\frac{\sqrt{2}{B}_0}{n\pi }·(e^{\frac{H_m\tau }{2}}-e^{\frac{-H_m\tau }{2}})·e^{\tau {(}^c{z}_0{-}^c{z}_p)}·\sin [\tau ·\cos \gamma ·{(}^c{x}_0{-}^c{x}_p)\hfill \\ \hfill & +\tau ·\sin \gamma ·{(}^c{y}_0{-}^c{y}_p)+\frac{3}{4}\pi ]\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill {}^c{\mathbf{B}}_{xz}=& \frac{-\sqrt{2}{B}_0}{n\pi }·(e^{\frac{H_m\tau }{2}}-e^{\frac{-H_m\tau }{2}})·e^{\tau {(}^c{z}_0{-}^c{z}_p)}·\cos [\tau ·\cos \gamma ·{(}^c{x}_0{-}^c{x}_p)\hfill \\ \hfill & +\tau ·\sin \gamma ·{(}^c{y}_0{-}^c{y}_p)+\frac{3}{4}\pi ]\hfill \end{array},$$

(7)

where $\gamma$ is the rotation around the z-axis shown in Figure 4.

In the same way, the magnetic flux density distribution along the y-axis ${}^c{\mathbf{B}}_y$ can be described as:

$$\begin{array}{cc}\hfill {}^c{\mathbf{B}}_{yx}=& -\sin \gamma ·\frac{\sqrt{2}{B}_0}{n\pi }·(e^{\frac{H_m\tau }{2}}-e^{\frac{-H_m\tau }{2}})·e^{\tau {(}^c{z}_0{-}^c{z}_p)}·\sin [\tau ·\cos \gamma ·{(}^c{y}_0{-}^c{y}_p)\hfill \\ \hfill & -\tau ·\sin \gamma ·{(}^c{x}_0{-}^c{x}_p)+\frac{3}{4}\pi ]\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {}^c{\mathbf{B}}_{yy}=& \cos \gamma ·\frac{\sqrt{2}{B}_0}{n\pi }·(e^{\frac{H_m\tau }{2}}-e^{\frac{-H_m\tau }{2}})·e^{\tau {(}^c{z}_0{-}^c{z}_p)}·\sin [\tau ·\cos \gamma ·{(}^c{y}_0{-}^c{y}_p)\hfill \\ \hfill & -\tau ·\sin \gamma ·{(}^c{x}_0{-}^c{x}_p)+\frac{3}{4}\pi ]\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill {}^c{\mathbf{B}}_{yz}=& \frac{-\sqrt{2}{B}_0}{n\pi }·(e^{\frac{H_m\tau }{2}}-e^{\frac{-H_m\tau }{2}})·e^{\tau {(}^c{z}_0{-}^c{z}_p)}·\cos [\tau ·\cos \gamma ·{(}^c{y}_0{-}^c{y}_p)\hfill \\ \hfill & -\tau ·\sin \gamma ·{(}^c{x}_0{-}^c{x}_p)+\frac{3}{4}\pi ]\hfill \end{array}.$$

(10)

3.3.2. Force and Torque Computation for the 1D Magnetic Field

In existing MLAs, the shapes of current-carrying coils are various such as circular, square, and rectangle coils. To obtain a universal wrench model of MLAs, the racetrack coil is selected as the analysis object for which the circular-, square-, or rectangle-shaped coils all can be indicated as a racetrack topology structure. The computational burden of the wrench model for the real-time simulation method of the MLA also should be considered, which affects the simulation efficiency. For this reason, the computational burden of this mathematical model should be as small as possible. In the 1D magnetic field, it is assumed that the magnetic array always moves in the area above the long straight sides of the racetrack coils in this model. In this assumption, the arc region and short sides of the coils are ignored as they are far away from the permanent magnet. The area where the integrated part of the coil is just the area covered by the magnetic array on the long straight sides. Compared with the existing magnetic force model [3], this way can significantly simplify the computation complexity and decrease the time delay for real-time simulation efficiently. According to the required range of the translation stroke and yaw angle of the MLA, the size of the coil and permanent magnetic array can be adjusted to meet this assumption.

According to the Lorentz law, the force and torque generated by the coil can be calculated. As the force on the 1D Halbach magnetic array is opposite the force on the coil, the force and torque can be expressed as:

$${}^c\overrightarrow{F}=-\int \int \int {-}^c\mathbf{J}{\times }^c\mathbf{B}dV,$$

(11)

$${}^c\overrightarrow{T}=-\int \int \int {-}^c\mathbf{r}{\times }^c\mathbf{J}{\times }^c\mathbf{B}dV.$$

(12)

where ${}^c\mathbf{J}$ is the current density in the coil; ${}^c\mathbf{B}$ is the magnetic flux density created by the 1D Halbach array; ${}^c\mathbf{r}$ is the arm of force, which can be expressed as:

$${}^c\mathbf{r}{=}^c\overrightarrow{x}{-}^c\overrightarrow{p}.$$

(13)

In the process of the 1D Halbach array’s motion, the coils’ area covered by the 1D Halbach array changes with the translation and rotation of the magnet. Hence, the volume integral of the covered area needs to be changed synchronizing with the change of the covered area. The changing area can be obtained by:

$$f\left(x\right)=\tan \gamma x-\tan \gamma (x_p-\frac{L_m}{2}\sin \gamma )+y_p+\frac{L_m}{2}\cos \gamma$$

(14)

and:

$$g\left(x\right)=\tan \gamma x-\tan \gamma (x_p+\frac{L_m}{2}\sin \gamma )+y_p+\frac{L_m}{2}\cos \gamma$$

(15)

which shows the expressions of the upper and lower edges of the magnetic array in the c coordinate system.

In this model, the $L_m$ is less than the $L_c$, which means the motion of the 1D Halbach array is always within the area directly above the long sides of coils shown in Figure 4a. Hence, the influence of the short and the arc region of the coil is neglected in the volume integration. In preview works, the coils were fully covered by the magnetic array, and the calculation of the force and torque needs to consider the arc and short region of the coils, which complicates the calculation and increases the computational burden [29]. Therefore, compared with existing models, this mathematical model for the 1D Halbach array not considering the arc and short region of the coil can not only have good calculation accuracy, but also decrease the computational burden.

Substituting (5)–(7), (14), and (15) into (11), the ${}^c\overrightarrow{F}$ for 1D magnetic array can be expressed as follows:

$$\left[\begin{array}{c}{}^c{F}_x\\ {}^c{F}_y\\ {}^c{F}_z\end{array}\right]=\sum _n^{1,5,9\dots }\Xi ·\frac{e^{{-}^c{z}_p}}{\pi }\left[\begin{array}{c}1\\ 0\\ \cos \tau \end{array}\right]\left[\begin{array}{c}\sin R\\ 0\\ \mathrm{cosR}\end{array}\right]J_y$$

(16)

where:

$$\begin{array}{cc}\hfill & \Xi =\frac{16\sqrt{2}{B}_0}{{\tau }^3\tan \gamma }\left(1-e^{-H_m\tau }\right)\left(1-e^{H_c\tau }\right)·\sin \left(\tau \frac{R_{\mathrm{c}}}{2\cos \gamma }\right)·\sin \left(\tau \frac{L_{\mathrm{m}}}{2}\sin \gamma \right)·\sin \left(\tau \frac{W_{\mathrm{c}}+R_{\mathrm{c}}}{2\cos \gamma }\right).\hfill \end{array}$$

(17)

$$R=\frac{\tau }{\cos \gamma }(\frac{W_c}{2}+R_c-x_p)+\frac{3}{4}\pi ,$$

(18)

and:

$$J_y=\frac{NI}{R_c{H}_c}$$

(19)

N represents the winding numbers of a coil, and I is the current in the coil.

The torque can be calculated as:

$${}^c\overrightarrow{T}={}^c\mathbf{r}\times {}^c\overrightarrow{F}.$$

(20)

3.3.3. Force and Torque Computation for the 2D Magnetic Field

Different from the former proposed magnetic force model for a 1D magnetic field, for the 2D magnetic field, the racetrack coil is completely covered in the magnetic field, as shown in Figure 4b, divided into four parts, as shown in Figure 5. Every part of the coil is subject to a magnetic force. The magnetic force acting on the coil is the resultant force of the magnetic field acting on four districts of the coil, and the magnetic field can be viewed as the superposition of the magnetic fields distributed along the x-axis and y-axis. Hence, the three components of the magnetic force for the 2D magnetic field can be expressed as:

$${}^c{\mathrm{F}}_z{=}^c{\mathrm{F}}_{xz1,2}{+}^c{\mathrm{F}}_{xz3,4}{+}^c{\mathrm{F}}_{yz1,2}{+}^c{\mathrm{F}}_{yz3,4}$$

(21)

$${}^c{\mathrm{F}}_x{=}^c{\mathrm{F}}_{xx1,2}{+}^c{\mathrm{F}}_{yx1,2}$$

(22)

$${}^c{\mathrm{F}}_y{=}^c{\mathrm{F}}_{xy3,4}{+}^c{\mathrm{F}}_{yy3,4}.$$

(23)

${}^c{\mathrm{F}}_{xz1,2}$ and ${}^c{\mathrm{F}}_{xz3,4}$, respectively, are the magnetic force in the z-direction generated by the magnetic array arranged along the x-axis acting on Part 1,2and Part 3,4. ${}^c{\mathrm{F}}_{yz1,2}$ and ${}^c{\mathrm{F}}_{yz3,4}$ are also the magnetic force in the z-direction, but they are generated by the magnetic array distributed along the y-direction. ${}^c{\mathrm{F}}_{xx1,2}$ and ${}^c{\mathrm{F}}_{xy3,4}$ are the force in the x-direction and y-direction generated by the magnetic array arranged along the x-axis. ${}^c{\mathrm{F}}_{yx1,2}$ and ${}^c{\mathrm{F}}_{yy3,4}$ are the force generated by the magnetic array arranged along the y-axis.

Through the Lorentz integral law, the magnetic force can be obtained:

$${}^c{F}_{xx1,2}=\sum _n^{1,5,9\dots }{K}_{xr1,2}·e^{{-}^c{z}_p\tau }·\sin [\tau ·\cos \gamma {·}^c{x}_p+\tau ·\sin \gamma {·}^c{y}_p-\frac{3}{4}\pi ]$$

(24)

$${}^c{F}_{xz1,2}=\sum _n^{1,5,9\dots }{K}_{xr1,2}·\cos \gamma ·e^{{-}^c{z}_p\tau }·\cos [\tau ·\cos \gamma {·}^c{x}_p+\tau ·\sin \gamma {·}^c{y}_p-\frac{3}{4}\pi ]$$

(25)

$${}^c{F}_{yz1,2}=\sum _n^{1,5,9\dots }{K}_{yr1,2}·\sin \gamma ·e^{{-}^c{z}_p\tau }·\cos [\tau ·\cos \gamma {·}^c{y}_p-\tau ·\sin \gamma {·}^c{x}_p-\frac{3}{4}\pi ]$$

(26)

$${}^c{F}_{yx1,2}=\sum _n^{1,5,9\dots }{K}_{yr1,2}·e^{{-}^c{z}_p\tau }·\sin [\tau ·\cos \gamma {·}^c{y}_p-\tau ·\sin \gamma {·}^c{x}_p-\frac{3}{4}\pi ]$$

(27)

$${}^c{F}_{xz3,4}=\sum _n^{1,5,9\dots }{K}_{xr3,4}·\sin \gamma ·e^{{-}^c{z}_p\tau }·\cos [\tau ·\cos \gamma {·}^c{x}_p+\tau ·\sin \gamma {·}^c{y}_p-\frac{3}{4}\pi ]$$

(28)

$${}^c{F}_{xy3,4}=\sum _n^{1,5,9\dots }{K}_{xr3,4}·e^{{-}^c{z}_p\tau }·\sin [\tau ·\cos \gamma {·}^c{x}_p+\tau ·\sin \gamma {·}^c{y}_p-\frac{3}{4}\pi ]$$

(29)

$${}^c{F}_{yz3,4}=\sum _n^{1,5,9\dots }{K}_{yr3,4}·\cos \gamma ·e^{{-}^c{z}_p\tau }·\cos [\tau ·\cos \gamma {·}^c{y}_p-\tau ·\sin \gamma {·}^c{x}_p-\frac{3}{4}\pi ]$$

(30)

$${}^c{F}_{yy3,4}=\sum _n^{1,5,9\dots }{K}_{yr3,4}·e^{{-}^c{z}_p\tau }·\sin [\tau ·\cos \gamma {·}^c{y}_p-\tau ·\sin \gamma {·}^c{x}_p-\frac{3}{4}\pi ]$$

(31)

where:

$$\begin{array}{cc}\hfill & K_{xr1,2}=\frac{16\sqrt{2}{B}_r}{{\tau }^3·n\pi }\left(e^{\frac{H_{\mathrm{m}}}{2}\tau }-e^{-\frac{H_{\mathrm{m}}}{2}\tau }\right)\left(e^{\frac{H_{\mathrm{c}}}{2}\tau }-e^{-\frac{H_{\mathrm{c}}}{2}\tau }\right)·\hfill \\ \hfill & \frac{\sin \left(\tau \frac{R_{\mathrm{c}}}{2}\cos \gamma \right)·\sin \left(\tau \frac{L_{\mathrm{c}}}{2}\sin \gamma \right)·\sin \left(\tau \frac{W_{\mathrm{c}}+R_{\mathrm{c}}}{2}\cos \gamma \right)}{\cos \gamma ·\sin \gamma }\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill & K_{yr1,2}=\frac{16\sqrt{2}{B}_r}{{\tau }^3·n\pi }\left(e^{\frac{H_{\mathrm{m}}}{2}\tau }-e^{-\frac{H_{\mathrm{m}}}{2}\tau }\right)\left(e^{\frac{H_{\mathrm{c}}}{2}\tau }-e^{-\frac{H_{\mathrm{c}}}{2}\tau }\right)·\hfill \\ \hfill & \frac{\sin \left(\tau \frac{R_{\mathrm{c}}}{2}\sin \gamma \right)·\sin \left(\tau \frac{L_{\mathrm{c}}}{2}\cos \gamma \right)·\sin \left(\tau \frac{W_{\mathrm{c}}+R_{\mathrm{c}}}{2}\sin \gamma \right)}{\cos \gamma ·\sin \gamma }\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill & K_{yr3,4}=\frac{16\sqrt{2}{B}_r}{{\tau }^3·n\pi }\left(e^{\frac{H_{\mathrm{m}}}{2}\tau }-e^{-\frac{H_{\mathrm{m}}}{2}\tau }\right)\left(e^{\frac{H_{\mathrm{c}}}{2}\tau }-e^{-\frac{H_{\mathrm{c}}}{2}\tau }\right)·\hfill \\ \hfill & \frac{\sin \left(\tau \frac{R_{\mathrm{c}}}{2}\cos \gamma \right)·\sin \left(\tau \frac{W_{\mathrm{c}}}{2}\sin \gamma \right)·\sin \left(\tau \frac{L_{\mathrm{c}}+R_{\mathrm{c}}}{2}\cos \gamma \right)}{\cos \gamma ·\sin \gamma }\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill & K_{xr3,4}=\frac{16\sqrt{2}{B}_r}{{\tau }^3·n\pi }\left(e^{\frac{H_{\mathrm{m}}}{2}\tau }-e^{-\frac{H_{\mathrm{m}}}{2}\tau }\right)\left(e^{\frac{H_{\mathrm{c}}}{2}\tau }-e^{-\frac{H_{\mathrm{c}}}{2}\tau }\right)·\hfill \\ \hfill & \frac{\sin \left(\tau \frac{R_{\mathrm{c}}}{2}\sin \gamma \right)·\sin \left(\tau \frac{W_{\mathrm{c}}}{2}\cos \gamma \right)·\sin \left(\tau \frac{L_{\mathrm{c}}+R_{\mathrm{c}}}{2}\sin \gamma \right)}{\cos \gamma ·\sin \gamma }.\hfill \end{array}$$

(35)

4. Validation of the Mathematical Model and the Simulation Model

4.1. Experimental Facility

The mathematical model proposed in this paper was verified by measuring the force and torque acting on the translator through a 6-DoF load cell (ATI mini40) located on a 4-DoF platform. The co-simulation model was validated via the measurement date of the translator position by the position sensors including eddy sensors and laser sensors. The measurement system consisted of the 6-DoF magnetically levitated actuator, the displacement detection system, the 6-DoF load cell, and the 4-DoF motion platform, as shown in Figure 6. In the structure of the magnetically levitated actuator, every two of the eight racetrack coils were fixed side by side on the stator, and four magnetic arrays were glued under the backplane of the translator. Each magnetic array corresponds to the two coils. The specification of each coil is AWG#22 (300 windings), and each Halbach array was made up of twelve NdFeB permanent magnets ($B_0$ = 1.28 T). The parameters of the magnet and coil are given in

Table 1. Parameters of the magnet and coil.

${\mathit{L}}_{\mathit{c}}$	${\mathit{W}}_{\mathit{c}}$	${\mathit{H}}_{\mathit{c}}$	${\mathit{R}}_{\mathit{c}}$	${\mathit{L}}_{\mathit{m}}$	${\mathit{W}}_{\mathit{m}}$	${\mathit{H}}_{\mathit{m}}$
40 mm	10 mm	10 mm	10 mm	60 mm	10 mm	10 mm

. In the displacement detection system, three laser sensors were used to measure the horizontal displacement of the backplane, and the vertical displacement was measured by the eddy sensors. The 4-DoF motion platform can control the translator with four Halbach arrays rotating round the z-axis and moving along the x-axis, y-axis, and z-axis. The yaw angle of the translator can be calculated by the measurement data of the sensors.

4.2. Validation of the Mathematical Model

The purpose of this experiment is to test the accuracy of this mathematical model. For the mathematical model of the 1D magnetic field, the magnetic force and torque were measured at three different rotation angles around the z-axis under static conditions. The total force and torque on the translator are the sums of the four units’ contributions where each unit is controlled separately. Therefore, measuring the contribution of only one unit can verify the accuracy of this mathematical model. In a unit, the two coils were controlled separately as well. Therefore, it was only necessary to measure the force generated by a coil on the magnetic array in a unit to verify the accuracy of the wrench model.

In this experiment, Coil 1 was energized shown in Figure 6. The position of Halbach array can be described with ${}^c\overrightarrow{p}$ and $\gamma$. At each rotation angle, the array starts at ${}^c\overrightarrow{p}=(0\mathrm{mm},0\mathrm{mm},13\mathrm{mm})$, and the force and torque are measured at 21 points along the $x-$ axis where $x_p\in [-10\mathrm{mm},10\mathrm{mm}]$. Coil 1 is excited by a current of 200 mA. Figure 7a–c shows the results of the measurement and computation at all points with three different yaw angles. Meanwhile, the force and torque values calculated by the boundary element model (BEM) are also shown in this figure to verify the mathematical model. The BEM is also a numerical method that has great calculation accuracy. It has a faster calculation speed than the FEM and has been used to calculate the magnetic force [30] via the software packages Mathematica and Radia. The force and torque curves illustrate that the mathematical model and the BEM almost coincide. Viewing the measurement results as a benchmark, the errors of the force and torque are also shown in Figure 7a–c. At three different rotation angles, the mathematical errors and the BEM errors are close compared with the measurement results. The maximum errors of force and torque between the mathematical model and BEM are not beyond 0.03 N and 0.007 N.m. Based on the measurement data, the force average error percentage of the proposed mathematical model is 9.07%, and that of the BEM is 8.95%. The average error percentages of torque are respectively 27.28% and 30.51%. These results verify the static accuracy of this mathematical model considering the rotation angles around the z-axis. Although the force and torque for the yaw angles around the x-axis and y-axis are not verified and considered, the yaw angle around the other two axes can be taken into consideration in the same way.

The validation result of the 2D mathematical model is shown in Figure 8. This figure shows the solution of the 2D mathematical model and the numerical solution of the boundary element method. This comparison illustrates that this 2D mathematical model has high calculation accuracy. In fact, for this superimposed 2D magnetic field, the magnetic field is the vector sum of two 1D magnetic fields, which means that the resultant magnetic force generated by this 2D field is also the vector sum of the magnetic force separately generated by two 1D magnetic fields, and the mathematical model for the 1D magnetic field is verified. In this case, although the magnetic force is not measured, we believe that the accuracy of the 2D mathematical model can still be explained by comparing the calculation results with the numerical solution.

In this experiment, the calculation time at one point was also measured. The calculation time of the proposed mathematical model is 0.07122 s, and that of the BEM is about 6 s. This result proves that the proposed mathematical is suitable for a real-time co-simulation model.

4.3. Validation of the Simulation Model

To verify the consistency of the simulation model results and the real control results, the co-simulation results were compared with the traditional simulation method and the measurement position data. The measurement data obtained by deploying the tested control strategy into a real MLA system are shown in Figure 2a. The traditional simulation method was also implemented based on MATLAB/Simulink, where the output of the control strategy is viewed as the real force and torque on the magnetic array. By integrating the output twice, the results of the traditional simulation method can be obtained.

In this case, the discrete-time PID controller was used to control the MLA to achieve the desired motion. However, a more complex controller also can be used for the control of the MLA system, which is not the main point of this paper. The design of the PID controller was done with the help of MATLAB SISOTOOL. To realize a fast response and aggressive transient behavior, the coefficients of this PID controller were set as $K_p=76.53$, $K_i=1151.8046$, and $K_d=1$, and the frequency was 4000 Hz.

In this experiment, the PID controller was given the desired position information where the ${}^c{\overrightarrow{p}}_{desire}=(0\mathrm{mm},0\mathrm{mm},13\mathrm{mm})$ and the rotation angles around the z-axis were $\gamma =5\mathrm{mrad}$ and $\gamma =50\mathrm{mrad}$. The results of this experiment are shown in Figure 9. It can be seen that compared with the traditional simulation model, the resulting curve of this simulation model proposed in this paper for the MLA is closer to the real control result curve. This phenomenon shows that the simulation model proposed in this paper has a high consistency with the real control results, and this simulation model for the MLA is reliable. Although the consistency of the co-simulation model with the 2D magnetic array cannot be verified by the measurement data, we believe it also has a high consistency with the real control system.

To further validate the accuracy of this proposed co-simulation model, the proposed mathematical model and the numerical model of the MLA were respectively utilized to play as the real plant in the simulator. In this test, the MLA whose starting position was ${}^s\overrightarrow{p}=(0\mathrm{mm},0\mathrm{mm},0\mathrm{mm},0\mathrm{rad},0\mathrm{rad},0\mathrm{rad})$ was controlled to track a trajectory, and the expressions of the trajectory in the ${}^s{x}^s{y}^{s}z$ plane shown in Figure 6 are given as:

$${}^{s}X=5sin\left(5\pi t\right),$$

(36)

$${}^{s}Y=5sin(5\pi t+\pi /2),$$

(37)

$${}^{s}Z=3,$$

(38)

while the trajectory of the rotation angle around ${}^{s}z$ are, respectively:

$$\gamma =\pi /180sin\left(10\pi t\right)$$

(39)

and:

$$\gamma =\pi /60sin\left(10\pi t\right)$$

(40)

Observing the result of this test shown in Figure 10, the error of the mathematical model and the numerical is almost zero. The average percentage deviation between the simulation results respectively based on the numerical model and mathematical model is shown in

Table 2. Average percentage deviation between the simulation results respectively based on numerical model and mathematical model.

Rotation Amplitude (rad)	x-axis Motion	y-axis Motion	z-axis Motion	z-axis Rotation
$\frac{\pi }{180}$	0.045%	0.012%	0.072%	0.15%
$\frac{\pi }{60}$	0.013%	0.067%	0.075%	0.066%

. This result highlights that the proposed simulation model based on the proposed mathematical model for the MLA has a high accuracy of the simulation result. The time cost of the co-simulation model based on the mathematical model is less than one second, and the numerical model needs about fifteen minutes, which can prove that the simulation model for the MLA based on the mathematical model is computationally-efficient. On the other hand, rotation around the z-axis will disturb the motion along the z-axis and the rotation around the x- and y-axes. As the rotation amplitude increases, the disturbance to other axes also increases. This phenomenon is caused by the coupling of the magnetic force and torque and can be improved by a more advanced controller.

5. Conclusions

This paper proposes a new co-simulation model for a 6-DoF magnetically levitated actuator. In this co-simulation model, the physical process of the MLA is described by the proposed mathematical model, and it can be used for 1D and 2D magnetic fields. In this mathematical model, the harmonic model is used to calculate the magnetic flux density, and the force and torque are calculated by the Lorenz integral law. Utilizing the coordinate transformation, the rotation of the moving part around the z-axis is taken into consideration when calculating the force and torque. This mathematical model is verified by comparing the calculation results of force and torque with the measurement results and numerical simulation software. Comparing the results shows that this model has good accuracy when considering the rotation. The performance of the co-simulation model for MLA is verified by comparing it with the traditional simulation model and measurements. The comparison result indicates that this simulation model has better performance and higher consistency with the real control system than the traditional simulation model for MLA. The co-simulation model based on the mathematical model is further validated by comparing it with the simulation model based on the numerical simulation. The comparison result shows that the mathematical model has not only high accuracy, but also high efficiency.