Friction torque field distribution of a permanent ‐ magnet spherical motor based on multi ‐ physical field coupling analysis

Non ‐ linear friction between the rotor and the rotor support mandrels would occur during the movement of the permanent ‐ magnet spherical motor (PMSpM). The friction which varies with the position causes PMSpM jitter or even stop it when running at low speeds, and the system operation accuracy is reduced. To solve these problems, this work proposes a parameter identification method of the Coulomb friction model. First, a dynamics model with the Coulomb friction of PMSpM is established by using generalised Euler angles and the Lagrangian energy method. Then, a parameter identification method of the Coulomb friction model based on COMSOL Multi physics coupled simulation and MEMS experiment is proposed to obtain the friction torque data of the PMSpM. Finally, the friction torque compensation in a closed ‐ loop control experiment is completed by using the friction data. The results show that the accuracy


| INTRODUCTION
With the development of modern industrial technologies, highprecision multi-degree-of-freedom (Multi-DOF) motion actuators that perform complex motions are being widely used. Most of the existing Multi-DOF motion actuators are composed of single-degree-of-freedom (Single-DOF) motion actuators through a complex transmission mechanism, with complex mechanical structures, huge volume and weight, and a large amount of non-linear friction. Therefore, the Multi-DOF spherical motor with a simple transmission structure and direct drive has attracted the attention of scholars at home and abroad. [1,2]. In recent years, the developed spherical motor mainly include types such as the variable reluctance type [3], permanentmagnet type [4][5][6], induction type [7,8], ultrasonic type [9] etc.
The permanent-magnet spherical motor (PMSpM) has a unique spherical rotor structure, and the output shaft can complete Multi-DOF movement within a certain range. The dynamic stability of the output torque and the trajectory tracking control accuracy are key technical parameters for evaluating the performance of the motor and for realising industrial applications. The above parameters are all affected by the friction torque. Hence, studying the friction torque field distribution of the motor is of great significance for friction torque compensation.
The contact between the rotor and the rotor support mandrels is the main cause of non-linear friction. Researchers try to change the support structure to reduce the non-linear friction term. The PMSpM designed by Kok-Meng Lee and others from the Georgia Institute of Technology in the United States adopts an SRJ ball bearing support structure [10,11]. In order to improve the load capacity of the motor, Gan Lei from the Harbin Institute of Technology and others improved the motor to a joint bearing support structure [4,12]. The outer rotor-type spherical motor developed by Masahito Tsukano and others from Osaka University in Japan adopted a universal joint support structure [13,14]. Yan Liang from the Beijing Aviation University and others designed the Halbach PMSpM with a roller support structure, imitating the support of a rolling ball in a mechanical mouse [15,16]. In the above research works, there is a lack of research on friction torque, resulting in errors between the simulated value of the output torque and the experimental value, and the accuracy of the trajectory tracking control is limited [17]. If the friction model is ignored and the maximum static friction torque is directly used for compensation, it will cause obvious overcompensation [18]. Therefore, the Multi-DOF motors represented by the PMSpM lacks a theoretical model of friction torque and an accurate compensation method, and they cannot provide data support for the design optimisation of the motor and the improvement of the trajectory tracking control accuracy.
In view of the above problems, this work uses the Coulomb friction model to theoretically model the friction torque of the PMSpM. Fit the dynamic data of the COMSOL Multi-physics (COMSOL Co, Ltd.) coupling simulation and MEMS experiment and identify the dynamic friction factor of the Coulomb friction model. Based on the identified parameters, the distribution of the friction torque field in the rotor motion range is obtained by simulation, which provides theoretical and data support for the design optimisation and friction compensation of the PMSpM.

| Basic structure
The research group designed and manufactured the PMSpM which is shown in Figure 1 [6]. The PMSpM adopts a mechanical structure with a stator spherical shell surrounding the rotor sphere.
The internal structure of the PMSpM after removing the stator shell is shown in Figure 2. The 40 sintered neodymium iron boron permanent magnets embedded on the rotor are equally divided into four layers, with an angle between layers of 30°and an angle between rows of 36°. Similarly, the 24 coils installed on the stator are equally divided into two layers; the angle between the layers is 45°, the angle between the rows is 30°and each coil is wound with 1200 turns. The polarities of any adjacent permanent magnets are opposite, and the coils are wound with 1200 turns. Each layer of permanent magnets and coils is equally spaced and symmetrically distributed on both sides of the equatorial plane. Based on this structure, the stator and rotor of the PMSpM have the same spherical centre, small rotational inertia and air gap, which can further improve the magnetic flux density in the air gap and the output torque [17]. Figure 3 shows two different configurations of the rotor support mandrel. The PMSpM in Figure 1 adopted the support mandrels shown in Figure 3(a). It is generally considered that the dynamic friction factor is a quantitative F I G U R E 1 Structure of the permanent-magnet spherical motor F I G U R E 2 Internal mechanism of the permanent-magnet spherical motor after removing the stator F I G U R E 3 Support mandrels with two different structures. (a) Smallball support mandrel (b) Big-ball support mandrel parameter related to the material and structure of the contact surface [19]. Under the premise of not changing the rotor structure, the friction torque optimisation can be achieved by optimising the material and structure of the support mandrels, which is an important direction for motor optimisation.

| Dynamic model
The rotor of the PMSpM performs Multi-DOF movement around the centre of rotation in the stator shell. In order to facilitate the description of the relative motion between the rotor and the stator, a fixed coordinate system (O-XYZ) is defined for the stator and a motion coordinate system (o-xyz) is defined for the rotor. Establish the coordinate system as shown in Figure 4, assuming that the origin O of the O-XYZ and the origin o of the o-xyz always coincide and the rotor is a rigid body a with constant centre of rotation.
Each motion of the rotor can be decomposed into three rotations from the O-XYZ to the o-xyz. The angles of the three rotations can be expressed as the tilt angle α, pitch angle β, and rotate angle γ. First, the rotor arrives at the x 1 y 1 z 1 coordinate system by rotating the α angle around the X-axis. Then, the coordinate system x 2 y 2 z 2 and xyz can also be achieved by rotating the β angle around the y 1 -axis and rotating the γ angle around the z 2 -axis to get the o-xyz coordinate system, respectively.
The above three rotation coordinate transformations can be expressed in the following matrix form: ð2Þ The transformation matrix from O-XYZ to o-xyz can be obtained by multiplying the above three matrices according to the transformation order, which is shown as The angular velocity in O-XYZ is denoted as ω i ¼ ½ω x ω y ω z � Τ , as shown in Figure 4. The rotor position at any time is described as θ i ¼ ½α β γ� Τ by the generalised Euler angle, and the angular velocity in the o-xyz is expressed as Denoting the element at the corresponding position in P O−o ∈ ℝ 3�3 as r ij ði; j ¼ 1; 2; 3Þ can obtain.
Substitute Equation (4) into Equation (5). Then the matrix can be expressed as cos β cos γ sin γ 0 −cos β sin γ cos γ 0 sin β 0 1 The centroid of the rotor coincides with the centre of gravity, so the rotor has no potential energy [20]. Based on the Lagrangian energy method [21], the dynamic equation of the rotor is established as where E is the total kinetic energy. During rotor movement, the control torque Q i ∈ ℝ 3 received by the rotor at the angular position θ i ∈ ℝ 3 can be expressed by Equation (7). Since the structure of the rotor is symmetrical and the mass of the output shaft and the MEMS can be ignored relative to the rotor, the rotational inertia components of the rotor have a relationship as follows: and the total kinetic energy of the rotor can be expressed as [22].
where _ ω x , _ ω y and _ ω z are obtained by taking the first derivative of the component of ω i ∈ ℝ 3 . Therefore, according to Equations (2)-(6), a dynamic model of the PMSpM is established as Among them, T c ∈ ℝ 3 is the friction torque described by the Coulomb model. M ðθÞ ∈ ℝ 3�3 and Cðθ; θ Þ ∈ ℝ 3�3 are the inertial matrix and Coriolis force matrix of the rotor, respectively [23].

| Friction torque modelling
The simplified support structure of the PMSpM is shown in Figure 5. The distance from the top of the output shaft to the centre of the rotor is R, and the radius of the rotor is r. The support mandrels are divided into two parts. The upper five support mandrels are symmetrically distributed on the equatorial plane of the motor, and their function is to prevent direct contact between the rotor and the stator. Meanwhile, the six support mandrels in the lower layer are in contact with the rotor and jointly support the rotor to complete the movement. The six contact surfaces between the lower support mandrels and the rotor are defined as a whole, and a comprehensive dynamic friction factor μ is defined to complete the Coulomb friction modelling.
In the dynamics model expressed in Equation (7), the control torque matrix Q i ∈ ℝ 3 of the rotor is calculated to obtain the electromagnetic force matrix The operation is as follows: The total mass of the rotor is m, the gravity and the electromagnetic forces received by the rotor finally act on the lower support mandrels in the form of pressure. Therefore, the Coulomb friction torque T c ∈ R 3 can be expressed as Generally, the change of the friction torque is complicated in the process of the rotor moving from rest to motion. Stribeck et al. found that the non-linear friction torque can be expressed as a combination of three processes through experiments. First, the maximum static friction torque needs to be overcome at the moment of PMSpM stars. Next, initiate the Coulomb friction torque phase during low speed. Finally, when the system runs at a high-speed, it exhibits a viscous friction torque that is positively correlated with the rotational angular velocity. All the above stages are continuous, which can be expressed as a function of angular velocity ω. The simplified Stribeck model is shown in Figure 6. In the low-speed motion stage, the friction torque decreases with increasing angular velocity from ω th1 to ω th2 , F I G U R E 5 Simplified support structure of the permanent-magnet spherical motor, including upper frictionless support mandrels and lower frictional support mandrels F I G U R E 6 Three stages of the simplified Stribeck friction torque model which is called the Stribeck effect [24]. From this, the friction torque expression can be established as where δ is the empirical constant, T brk is the maximum static friction torque, T c is the Coulomb friction torque, and B is the viscous friction coefficient. When the maximum static friction torque is reached, the angular velocity is ω th1 and ω th1 is infinitely close to 0. Due to the non-linearity of the friction torque when the PMSpM moves, the friction torque will hardly reach the linear stage at high speed [24]. Therefore, when the motor is moving, the angular velocity of the PMSpM is maintained around ω th2 , with the friction torque T f ≈ T c . As mentioned in Equation (11), T f is only the Coulomb friction torque.

| IDENTIFICATION OF THE DYNAMIC FRICTION FACTOR AND THE FRICTION TORQUE FIELD DISTRIBUTION
In this section, the dynamics data of the tilting motion are obtained through simulation and experiment respectively, where the AC total current of 3A is used to drive the PMSpM. In the simulation and experiment, the same energisation strategies are adopted to ensure that the trajectories in the simulation are consistent with that in the experiment. Input the pre-determined tilting motion trajectories' parameters into the upper computer, and analyse the energisation strategies of 24 coils by using the pseudo-inverse matrix [20]. Figure 7 shows the tilting motion trajectories. After completing a tilting motion along the previous trajectory, step 3.6°in the direction of rotation motion to complete the next tilting motion. Limited by the structure of the PMSpM, the theoretical maximum tilting angle of the output shaft is 37.5°.
Hence, the angle of each tilting motion from the initial position is 37.5°.

| Tilting motion experiment with MEMS
The experimental equipment is mainly composed of the PMSpM, upper computer, current source, control circuit and the MEMS sensor, as shown in Figure 8. The sensor is mounted on the top of the output shaft [25], which can obtain the angular displacement θ, angular velocity _ θ and the angular acceleration € θ components of the rotor under a given energisation strategy in real time. The acquired data is transmitted to the upper computer via a wireless data transmission module.
The sampling frequency of the sensor is 100 Hz, and the sampling period is T period ¼ 0:01 s. Therefore, the experimental data can be expressed as a discrete function of time, and the forward difference operation can be performed on this discrete function. The sequence of the sampling points is numbered with i by time. The increment of the function in each sampling period is Δf ðiÞ ¼ f ði þ 1Þ − f ðiÞ, and the first-order difference equation of the discrete function can be obtained as Therefore, based on the angular displacement θ collected by the sensor, the linear velocity v, linear acceleration α and angular acceleration _ ω of the rotor in O-XYZ can be obtained by vector operation and difference operation respectively [18]. In fact, v, α and _ ω are all vectors, but in the subsequent parameter identification process, these vectors need to be decomposed into mutually independent component forms. The components of the above vectors are all discrete functions of time, so they can be drawn as time-varying curves.
At the beginning of this chapter, the tilting motion trajectory specified in the experiment as shown in Figure 7, is introduced. At the same time, it introduces the method of using the upper computer to analyse the energisation strategy corresponding to a given trajectory. Based on the above preparation work, a total of 100 tilting motion experiments were completed. Before each experiment starts, adjust the initial position of the tilting motion through the vertical restraint device of the output shaft, and the o-xyz is also adjusted to coincide with the O-XYZ. In the experiment, a total of 100 sets of the deflection motion experimental data were collected.

F I G U R E 7
Multiple tilting motion trajectories in simulation and experiment LI ET AL.

| Tilting motion simulation with COMSOL
The advantage of COMSOL is that it can complete the multiphysics coupling simulation. In order to simulate the actual working state of the PMSpM, the coupling simulation of the magnetic field and the mechanical motion are simultaneously completed in the AC/DC module and the multi-body dynamics module.
Select a tilting motion trajectory from Figure 7 and analyse the corresponding energisation strategy. As in the experiment, a total of 3 A current needs to be passed through two coils with a difference of 180°in longitude. The electromagnetic force is generated by the permanent magnet under the action of the coil, and the rotor overcomes the non-linear friction between the support mandrels and relatively rotates.
At the same time, use spherical joint plug-in to characterise the nonlinear friction torque experienced by the rotor. There is a fixed axis e and a moving axis e d in the spherical joint plug-in. The supporting mandrels corresponds to the axis e, and the rotor corresponds to the axis e d that can complete the Multi-DOF space motion. Under this definition, the friction torque can be controlled by adjusting the dynamic friction factor μ between the source axis and the target axis of the sphere pair.
In COMSOL, the transient calculation is performed in the multi-physics coupling simulation according to the previous settings. During the simulation, the rotor completes the tilting motion under the combined action of the electromagnetic force and friction force, the two modules run simultaneously and there is data interaction. As a result, the magnetic field coupling between the permanent magnet and the coil at different times is shown in Figure 9. The vector's magnetic potential line in the figure tends to find the shortest closed path, thereby generating electromagnetic torque. From the beginning of the rotor movement in Figure 9(a) until the maximum deflection angle is reached in Figure 9(c), it took less than 0.4 s. Follow the same energisation strategy shown in Figure 9. The tilting motion trajectory of the PMSpM is shown in Figure 10. In this process, the electromagnetic torque and friction torque of each movement point can be collected, and the dynamic data v, a and _ ω can also be obtained. Similar to the processing method of the experimental data in Section 3.1, the dynamic data obtained by COMSOL is drawn as a timevarying curve for subsequent parameter identification.
As the energisation strategy changes, the PMSpM will complete the tilting motion along different trajectories, resulting in different friction torque. Thus, based on the same energisation strategies as the MEMS experiment, the multiphysical coupling simulation of the PMSpM is completed to obtain the dynamics data.
After completing the basic settings of multi-physics, encapsulate the COMSOL model as an APP, which can significantly improve the simulation speed. In the APP, only 25 variable entries are reserved of which 24 are used to input the coil current matrix and the remaining one is used to input the dynamic friction factor. Thus, the multiple tilting motion simulation can be completed quickly.
The main moment of inertia is needed in some calculations, which can also be done in COMSOL. In order to ensure calculation accuracy, it is necessary to accurately establish the geometric model of the PMSpM and add material properties. The calculated component of the moment of inertia J x , J y and J z is shown in Table 1.

| Data-based parameter identification
By adopting the same energisation strategy, the PMSpM can be driven to make a tilting motion along the same trajectory in the experiment and simulation. Therefore, the MEMS sensor and the COMSOL software can obtain similar dynamic data, respectively, under the same motion trajectory. In actual operation, the ultimate goal of the simulation is to obtain dynamic data that are extremely similar to those in the experiment. To achieve this goal, it is necessary to constantly adjust the dynamic friction factor μ between the rotor and the supporting mandrels in the simulation to change the dynamic data.
In this process, the simulation and experimental data are drawn as time-varying curves, and the simulation curve approaches the experimental curve gradually. Therefore, when the aforementioned approximation reaches a certain level, it can be considered that μ in the simulation is the identified dynamic friction factor of the motor. The flow chart of parameter identification is shown in Figure 11.
A reliable initial value can significantly reduce the number of simulations in parameter identification; thus, the selection of the initial value μ 0 of the dynamic friction factor in the simulation is important. Aiming at the situation where a hard sphere rolls on a soft surface, article [26] proposes a method for estimating the dynamic friction factor, which is typically about 0.1 or less. Furthermore, the maximum static friction torque T brk of the PMSpM was measured with a spring dynamometer, which was about 0:1 N ⋅ m [18]. The estimation formula of the PMSpM parameter μ 0 is deduced as the initial value of the dynamic friction factor of PMSpM μ 0 ≈ 0:025.
Extract a set of linear velocity component data of the tilting motion experiment, and draw the time-varying curve of the X-axis as shown in Figure 12. Input the initial value μ ¼ 0:025 of the dynamic friction factor in the simulation, and adopt the same energisation strategy as the experiment. Because the maximum static friction force F brk is greater than the dynamic friction force F f , the given initial value μ is greater than the actual value. Based on the above analysis, the value of μ needs to be gradually reduced in the following simulation.
Therefore, through the adjustment of μ, the simulation timevarying curves of the linear velocity component are infinitely close to the experimental value curve. The process of the simulation curve gradually approaching the experimental curve is shown in Figure 12; the result of this identification is μ = 0.0193.
Similar to the above process, continue to perform parameter identification on the components of v, α and _ ω obtained by simulation and experiment, and each component can get an identification value. The specific identification process is shown in Figure 13, from (a) to (c). In the above three sets of simulated data the components approach the experimental data components. nine The dynamic friction factor μ can be obtained as shown in Table 2.
To verify the rationality of each identification result, take the experimental data as the benchmark, define the difference between the simulation data and the benchmark at the same time as the residual data to observe the effect of each identification. Meanwhile, define −10% and 10% of the peak  Figure 13(d). When most of the residual data points fall within the confidence interval between the residual observation lines, it can be considered that the identification process has obtained an accurate dynamic friction factor μ. If the residual distribution is not reasonable, the dynamic friction factor can be adjusted according to the trend of the simulation value curves until the residual distribution is reasonable. According to the above conditions, the identification results in Table 2 are valid and credible. Based on all the identified results, a reasonable dynamic friction factor μ = 0.0194 is determined by the root mean square of all the identified values.
Fix the dynamic friction factor of the Coulomb friction model in COMSOL to μ = 0.0194 and control the PMSpM to complete the tilting motion along the trajectory in Figure 7. As in the experiment, the energisation strategy in the simulation is analysed by the upper computer. Extract the friction torque T f corresponding to each trajectory in COMSOL, and obtain the distribution of the friction torque field of the PMSpM, as shown in Figure 14.

| TRAJECTORY TRACKING VERIFICATION EXPERIMENT
In the third section, on the basis of the parameter μ = 0.0194, the friction torque field distribution of the PMSpM in Figure 14 is obtained. To verify the validity of the field distribution results, design a closed-loop trajectory tracking experiment with friction torque compensation. Each position θ i ¼ ½α β γ� Τ in Figure 14 corresponds to a uniquely determined friction torque value T f . Therefore, the compensation of the friction torque is determined by the position of the rotor in trajectory tracking. First, set the expected trajectory in the experiment as x d . The actual trajectory obtained by tracking is x s , and the angle error δ between x d and x s is obtained. Second, the PD controller calculated the control torque T d , which is required to correct the tracking error. Finally, the vector sum of the friction torque compensation amount T f and the control torque T d is solved and substituted into the pseudo-inverse matrix to obtain the required energisation strategy. Based on the above analysis, the system block diagram is shown in Figure 15.
Define the desired trajectory in the upper computer, which can be expressed as To verify the trajectory tracking effect after friction compensation, adjusted the initial position of the output shaft before the experiment to make it deviate from the desired trajectory. The initial position is After completing the above settings, the PMSpM is controlled to complete trajectory tracking under the premise of a total current of 3 A.
Decomposing the tracking trajectory into the three Euler angle directions of α, β and γ, which can more intuitively show the compensation effect of the friction torque field distribution. Figure 16 shows the trajectory tracking curves in three directions. In the above figures, the tracking effect after friction compensation is significantly better than that before compensation, mainly reflected in the tracking speed and accuracy. To analyse the above two indicators, the angle errors of the three Euler Angle directions in the trajectory tracking process were obtained through calculation, as shown in Figure 17.
Perform quantitative analysis on Figures 16 and 17, and get the response time and tracking error in the α, β and γ direction as shown in Table 3. Compared with the previous compensation, the response time of α, β and γ directions entering stable tracking after friction compensation is reduced by 44.83%, 44.01% and 43.48% respectively. At the same time, during the tracking time of t ∈ ½0; 8� seconds, the tracking error after friction compensation was reduced by 52.35%, 52.87% and 42.75%, respectively.
Because the error is artificially set at the starting point of trajectory tracking, and the errors of the three Euler angle directions are 10°, 18°and 30°respectively, the error is large in the time period t ∈ ½0; 2� seconds. It can be seen from Figure 16 that the trajectory tracking becomes stable in the α, β and γ directions after 2 s. Analyse the tracking angle error in the period t ∈ ½2; 8� seconds in Figure 17(b) and drop to 1.17°, 1.29°and 0.30°, respectively.
The experimental results show that the friction torque field is used to determine the compensation value in the process of trajectory tracking, and the tracking trajectory convergence speed and accuracy are significantly improved. In summary, the method of using COMSOL simulation and MEMS experiment to identify the dynamic friction factor and obtain the friction torque field distribution is effective.

| CONCLUSION
This work presents a method to obtain the dynamic friction factor of the contact surface between the rotor and the support mandrels by combining the experiment and simulation and then obtaining the distribution of the friction torque field. The accuracy and validity of the friction torque field distribution are proved, which provides the reference data for precise tracking control. In terms of structural optimisation design, the distribution diagram of the friction torque field can intuitively show the position where the friction torque is larger, so it can provide empirical guidance for the arrangement number and position of the support mandrels. Besides, the coil and the permanent magnet can be optimised under the premise that the rotor and support mandrels' structures remain unchanged. Based on the above conditions, the dynamic friction factor μ is fixed; it can be coupled by the multi-physics simulation to optimise the effect of the magnetic circuit.
The steps to obtain the distribution of the friction torque field through multi-physical field coupling simulation and experiment are complex and slow, so better means and methods are needed to improve efficiency. The main workload of obtaining the friction torque field distribution is in the simulation, so some methods are needed to improve the efficiency of the simulation. On the one hand, based on the structural symmetry characteristics of the PMSpM, the simulation model is simplified through the sector symmetry function in COMSOL, which can reduce simulation time. On the other hand, encapsulating the simulation model into an APP and cooperating with the energisation strategy of the upper computer can achieve efficient simulation.