Magnetic field modeling and validation for a spherical actuator with cylindrical permanent magnets

Magnetic field modeling and validation for a spherical actuator with cylindrical permanent magnets

Accepted Manuscript

Magnetic field modeling and validation for a spherical actuator with cylindrical permanent magnets Jingmeng Liu, Xuerong Li, Weihai Chen, Lu Liu, Shaoping Bai PII: DOI: Article Number: Reference:

S1569-190X(19)30087-5 https://doi.org/10.1016/j.simpat.2019.101954 101954 SIMPAT 101954

To appear in:

Simulation Modelling Practice and Theory

Received date: Revised date: Accepted date:

19 March 2019 2 July 2019 22 July 2019

Please cite this article as: Jingmeng Liu, Xuerong Li, Weihai Chen, Lu Liu, Shaoping Bai, Magnetic field modeling and validation for a spherical actuator with cylindrical permanent magnets, Simulation Modelling Practice and Theory (2019), doi: https://doi.org/10.1016/j.simpat.2019.101954

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A permanent magnet spherical actuator (PMSA) which can generate 3-

CR IP T

DOF rotations in one driving unit is designed. • The complex magnetic field distribution of PMSA is formulated analytically with current model by toroidal functions.

• The influence of PM structure is investigated using the analytical model.

• An experimental platform is designed to measure the 3-D magnetic field

AC

CE

PT

ED

M

AN US

of PMSA for validation.

1

ACCEPTED MANUSCRIPT

Magnetic field modeling and validation for a spherical actuator with cylindrical permanent magnets

CR IP T

Jingmeng Liua , Xuerong Lia,b , Weihai Chena,∗, Lu Liua , Shaoping Baib a School

of Automation Science and Electrical Engineering, Beihang University, 100191 Beijing, China b Department of Materials and Production, Aalborg University, 9220 Aalborg, Denmark

Abstract

AN US

A permanent magnet spherical actuator is an electromagnetic device character-

ized by achieving multi-degree-of-freedom rotations in one driving unit. Analytical magnetic field modeling of the rotor array is needed for the design, performance analysis and optimization of the spherical actuator. This paper presents an analytical method to formulate the complex magnetic field distribution of

M

rotor array with cylindrical permanent magnets in 3D spherical space. This method combines the current model by toroidal functions and transformation method, thus the magnetic field can be calculated efficiently with reasonable

ED

accuracy in order for further design optimization and real-time control. The analytical magnetic model of the permanent magnet rotor array is validated with the numerical finite element method by Ansoft Maxwell and experimental

PT

measurement.

Keywords: Permanent magnet spherical actuator, magnetic modeling,

CE

cylindrical permanent magnet, current model, toroidal functions.

AC

1. Introduction A spherical actuator is a device which can provide multi-degree-of-freedom

(DOF) rotations in one driving unit [1, 2]. Compared with spherical parallel ∗ Corresponding

author. Tel.: +86 10 8231 5560; fax: +86 10 8231 5920. Email address: [email protected] (Weihai Chen)

Preprint submitted to Simulation Modelling Practice and Theory

July 22, 2019

ACCEPTED MANUSCRIPT

manipulators [3, 4, 5], spherical actuators are advantageous due to simple me5

chanical structure and fast response, which can find many potential applications in robot manipulators [6, 7], aerospace [8, 9] and camera systems [10]. So far,

CR IP T

various designs and working principles of 3-DOF spherical actuators have been proposed such as spherical ultrasonic motors [11], spherical induction actuators

[12, 13] and permanent magnet spherical actuators [14]. Among them, perma10

nent magnet spherical actuators are widely investigated. A permanent magnet

spherical actuator mainly consists of a rotor constructed by spherical permanent magnet array and a set of multi-layer stator windings distributed in a

AN US

spherical shape. By energizing the stator coils with appropriate current control

strategies, the rotor can generate specified 3-DOF motion on the sphere. The 15

performance of the electromagnetic spherical actuator is greatly influenced by the air-gap magnetic flux density generated by rotor permanent magnet array [15]. Therefore, an effective analytical magnetic model of the rotor is required for design and control of the spherical actuators.

20

M

Up to now, most rotor arrays are mainly constructed by dihedral poles and cylindrical permanent magnet poles. Several works on magnetic modeling with

ED

both types of rotor arrays have been proposed. Magnetic equivalent circuit method is used to analyze the magnetic field of spherical actuators constructed with dihedral permanent magnets in 3D space [16]. The geometry of the device

25

PT

is meshed into elements and a magnetic circuit is employed based on different magnetic boundary conditions. However, this method is a kind of lumped-

CE

parameters analysis with which requirements of accuracy and computation time cannot be met simultaneously. Lee and Son et al. developed a distributed multipole model which was originated from the concept of a magnetic dipole to

AC

calculate the magnetic field of cylindrical permanent magnets with axial mag-

30

netization [17, 18]. In this method, the permanent magnets are replaced with magnetic dipoles, and then the magnetic scalar potential at any point can be obtained by the superposition of all the magnetic dipoles. Spherical harmonics is an analytical method for the design and analysis of spherical actuators with dihedral permanent magnets [19, 20, 21]. Spherical harmonics method can be 3

ACCEPTED MANUSCRIPT

35

also used to analyze approximately the magnetic field distribution of cylindrical permanent magnets [22]. In [23], the analytical torque model of spherical actuators with cylindrical permanent magnets is obtained by adjusting the an-

CR IP T

alytical terms with the measured experimental results of the actual magnetic field. The charge model and current models are other alternative choices to 40

calculate the air-gap magnetic flux density of permanent magnets with better accuracy. In the models, the permanent magnets are reduced to a distribution of equivalent ”magnetic charges” or equivalent currents [24, 25, 26]. The ge-

ometry is taken into consideration and the analytical model is obtained using

45

AN US

the volume and surface integration. However, geometry-dependent integration increases the computation complexity [27, 28]. To simplify the computation problem, the magnetostatic potential expressed in terms of elliptic integrals is used for uniformly magnetized cylindrical permanent magnets [29, 30]. In [31], the toroidal harmonics is proposed to compute the 3D magnetic field from a bipolar cylindrical permanent magnet. It is completely an analytical method which is useful for uniform magnetization vectors [32].

M

50

Compared with dihedral poles, a design of cylindrical permanent magnet has

ED

the advantages of easy fabrication and magnetization, and its use in the actuator can reduce the inertial moment of the rotor in order to improve the dynamic response [23, 33]. Therefore, cylindrical permanent magnets are proposed as the rotor poles in this study. In our previous work, the 3-DOF permanent magnet

PT

55

spherical actuator consists of 8 cylindrical permanent magnets in one layer and

CE

30 stator coils in three layers [34]. The magnetic field is analyzed numerically with finite element method, and then the torque is obtained by curve fitting method [35]. The finite element method can yield good accuracy with consid-

eration of nonlinear parameters and complicated geometry shapes. However,

AC

60

it is time-consuming, especially for design and parameters optimization of the structure [36], and it isn’t favourable for the real-time control of the spherical actuator. This paper proposes a new design of spherical actuator that is constructed

65

by a two-layer rotor array with 16 permanent magnets in total, and 36 sur4

ACCEPTED MANUSCRIPT

rounding stator coils to enhance the output torque. The objective of this paper is to develop an analytical method for the analysis of the magnetic field distribution of the 3D magnetic array with cylindrical permanent magnets. The

70

CR IP T

proposed method integrates the current model by toroidal functions and the transformation method together. The magnetic field model of the single cylindrical permanent magnet is first formulated analytically using the current model

by toroidal functions in cylindrical coordinates. Utilizing coordinate transformation method, the magnetic field model is achieved in spherical coordinate

system, and then the complete magnetic field distribution of the rotor array is obtained by superposition principle. Compared with the conventional current

AN US

75

model, the toroidal analysis reduces the computational complexity of the integral function. Moreover, the introduction of the transformation method helps to formulate model of cylindrical permanent magnet at any position in the global coordinate system without the loss of accuracy.

The rest of the paper is organized as follows. Section 2 presents the design

80

M

and working principle of the spherical actuator. Section 3 formulates the analytical model of the magnetic field distribution of a single permanent magnet

ED

in cylindrical coordinate. Section 4 presents the magnetic field distribution of the whole rotor structure in spherical coordinate system. Section 5 illustrates 85

the simulation results, and further analyzes the effect of the permanent magnet

PT

structure parameters on magnetic field. Section 6 presents the experiments.

CE

Finally, the paper is concluded in Section 7.

2. Design and working principle The new design of the permanent magnet spherical actuator is shown in

Fig. 1. It consists of a rotor and a shell-shaped stator. This rotor can generate

AC 90

continuous spinning motion about its output shaft, and incline about its equatorial plane to a position that the axes of a permanent magnet pole and the coil are aligned. The maximum tilting angle is ±15◦ . The specification of the spherical actuator is listed in Table 1.

5

ACCEPTED MANUSCRIPT

Shaft

Stator

Coil

CR IP T

Rotor

PM

AN US

(a) Prototype

Iron hoop

Rotor shell

z

x

y

ED

M

PM

(b) CAD model of rotor structure

PT

Figure 1: Design of the permanent magnet (PM) spherical actuator

Parameter

Value

Outer/inner stator radius

Ros = 98 mm, Ris = 84 mm

Outer/inner rotor radius

h1 = 50 mm, h2 = 34 mm

Cylindrical permanent magnet

R = 10 mm, l = 16 mm

Number of coils

Nc = 36/3 layers

Number of permanent magnet poles

Np = 16/2 layers

AC

CE

Table 1: Design parameters of the spherical actuator

6

ACCEPTED MANUSCRIPT

There are in total 36 air-core coils mounted on the outer stator shell in three

95

layers and all the stator coils point to the center of the sphere. The specific structure of the rotor is shown in Fig. 1(b). The permanent magnets are made

CR IP T

of NdFeB-N35UH which is a kind of rare earth material with high coercive force. The cylindrical permanent magnets are radial-symmetric distributed on 100

the rotor surface.

The working principle is illustrated in Fig. 2. The rotor motion is generated by the repelling and attracting electromagnetic forces between the rotor permanent magnets and stator coils. The rotor can produce corresponding motion

105

AN US

within the workspace by activating the input currents of stator coils in a specific

approach. The spinning motion is governed by energizing all the circumferential coils of both upper and lower layers in order as shown in Fig. 2(a). With coils activated in longitudinal directions, the rotor can realize tilting motion as illustrated in Fig. 2(b). To tilt the end-effector of the rotor from the initial balanced state, the upper and lower coils are energized first. As the the rotor inclines to a certain angle, the middle-layer coil begins to involve in the actuation. When

M

110

the upper and middle coils are aligned with the double-layer permanent magnet

ED

poles, respectively, the end-effector then reaches its maximum tilt motion. In our previous work, the rotor permanent magnet array was constructed with eight cylindrical permanent magnets in one layer [37]. For the spinning motion, only middle-layer stator coils are employed, and the tilting motion is

PT

115

implemented by the permanent magnet interacted with the upper and lower

CE

coils. Compared to spherical actuator with the one-layer rotor array, the newly designed spherical actuator can make more permanent magnets and stator coils involved in the electromagnetic driving operation effectively, which can improve the torque output to implement the 3-DOF spherical rotations.

AC

120

7

Spinning motion

PM

CR IP T

ACCEPTED MANUSCRIPT

Attraction

N

Attraction

S

S

PM z

x N

S

N

S

N

Coil

N

y

S

y S

z

AN US

N

x

N

N

N

S

N

Coil S

Repulsion

Repulsion

PM

Attraction z

N

S

N

S

y S

CE Figure 2:

z

Attraction

Coil x

y

N

Repulsion

PM

z

Coil

x

S

N

N

Tilting motion

Tilting motion

Attraction

PM

PT

x N

ED

Tilting motion

M

(a)

y

Coil

Repulsion

Repulsion

(b)

Working principle of the permanent magnet spherical actuator in generating (a)

AC

spinning motion and (b) tilting motion

8

ACCEPTED MANUSCRIPT

z0(z’)

0

yk(y’)

0

0

θk

r(kr)

0

k

CR IP T

zk

o

y0

0

ϕk

x’

AN US

xk

x0 Figure 3: kth

P(kP)

Coordinate rotations from global spherical frame to the local cylindrical frame of

permanent magnet

M

3. Magnetic field distribution of a cylinder permanent magnet in cylindrical coordinate system

ED

The rotor array is constructed by cylindrical permanent magnets that are distributed equally in two layers. The permanent magnets on the same latitude are uniformly magnetized with alternative variation of N and S poles, and the

PT

permanent magnet poles on the same longitude are with the same parallel magnetization. The specific magnetization axis of k th permanent magnet pole in

CE

the global rotor spherical frame is given by the unit vector 0 ˆrk which is defined

AC

as



  cos(0 φk ) sin(0 θk )         0 ˆrk =  0 yk  = (−1)k  sin(0 φk ) sin(0 θk )      0 cos(0 θk ) zk 0

xk



where k = 1, 2, ..., Np presents the k th permanent magnet, 0 φk = and 0 θk =

125

π 2

−

π 12

· sgn(

Np +1 2

(1)

π 8

+ π4 (k − 1)

− k) are the azimuth and zenith angles of the k th

permanent magnet in global spherical coordinate. Every permanent magnet pole

9

ACCEPTED MANUSCRIPT

is modeled separately in its cylindrical coordinate system. To analyze the magnetic distribution at point 0 P(0ρ ,0 θ ,0 φ ) in global spherical coordinate system, P needs to be transformed into each local cylindrical coordinate system as

k

P(kρ ,k φ ,k z). By transforming further the magnetic field distribution k B from

CR IP T

130

0

the local frame of every permanent magnet pole to global spherical frame 0 Bk , the magnetic field in the global frame is obtained. The coordinate frame rotations between the global and local frame are shown as Fig. 3. 3.1. The current model

R

AN US

zk

A

k

ρ'

k

z'

ED

h1

k

P(k²ˈk¶ˈkz)

P'

M

r'

ρ

r

k

z

CE

PT

h2

yk

¶'

k

¶

xk

AC

Figure 4: A single permanent magnet pole in kth local cylindrical coordinate

135

The magnetic field of a single permanent magnet pole is analyzed using

the current model. The cylindrical permanent magnet in its local cylindrical coordinate is shown in Fig. 4. The derivation of the current model starts with

10

ACCEPTED MANUSCRIPT

the magnetostatic field equations for current-free regions, (2)

∇ · kB = 0

(3)

CR IP T

∇ × kH = 0

where k presents the local coordinate system of the k th permanent magnet,

k = 1, 2, ..., Np , k H is the magnetic field strength and k B is the magnetic flux density in the k th local coordinate. The magnetic flux density, k B in Eq. (3), can be expressed as the curl of the magnetic vector potential field k A B = ∇ × kA

AN US

k

(4)

where k A is guaranteed by the coulomb gauge condition ∇ · kA = 0

(5)

B = µ0 (k H + k M)

(6)

The constitutive law is

M

k

where k M is the magnetization and µ0 is the permeability of the permanent

ED

magnet pole. Substituting Eqs. (4) and (6) into Eq. (2) ∇2 · k A = −µ0 (∇ × k M)

(7)

CE

tion

PT

Equation (7) can be written in integral form using the free-space Green’s funcG(k r, r0 ) = −

Then, we find that k

µ0 A( r) = 4π k

Z

1 4π|k r − r0 |

(8)

∇ × kM 0 dV |k r − r0 |

(9)

AC

where r0 is the position vector of the source point P0 and k r represents the position vector of the observation point k P. If the magnetization k M is confined to a volume V , and falls abruptly to zero outside of V , Eq. (9) becomes Z k I k µ0 Jm (r0 ) µ0 jm (r0 ) 0 k 0 A(k r) = dV + ds 4π V |k r − r0 | 4π S |k r − r0 | 11

(10)

ACCEPTED MANUSCRIPT

where S is the surface of the magnet, and Jm and jm are equivalent volume and surface current densities, respectively, which are defined as J m = ∇0 × k M

(A/m2 )

k

ˆ0 jm = k M × n

(A/m)

(11)

CR IP T

k

ˆ 0 is unit normal vector of the permanent magnet surfaces. The cylindrical where n permanent magnet is polarized along its axis with a uniform magnetization k

ˆ M = Ms k z

(12)

AN US

The equivalent current densities are computed using Eq. (11). We readily determine that k

J m = ∇0 × k M = 0

M

To evaluate the surface term, we first evaluate the unit surface normals    ˆ0 , z 0 = h1 z   0 ˆ = n ρˆ0 , ρ0 = R     z ˆ0 , z 0 = h

ED

and then compute

k

(13)

(14)

2

ˆ 0 , (ρ0 = R) ˆ 0 = Msφ jm = k M × n

(15)

PT

Equation (10) becomes k

A(k r) =

µ0 Ms 4π

I

s

ˆ0 φ ds0 , (ρ0 = R) |k r − r0 |

(16)

CE

3.2. Magnetic field analysis by toroidal functions In cylindrical coordinate, the reciprocal distance between the source point

0

AC

P and the observation point k P is given by 1 =q |k r − r 0 | k

1 02

(17)

ρ2 + ρ + (k z − z 0 )2 + 2 · k ρρ0 cos(k φ − φ0 )

The Green’s function expansion of Eq. (17) is given by ∞ X 1 1 p = m Qm− 12 (k β) cos[m(k φ − φ0 )] |k r − r0 | π k ρρ0 m=0

12

(18)

ACCEPTED MANUSCRIPT

2

2

where k β = (k ρ + ρ0 + (k z − z 0 )2 )/(2 · k ρρ0 ) > 1. The coefficient m is called

the Neumann factor, which is equal to 1 for m = 0 and 2 for m ≥ 1. Qm− 21 (k β)

is called half-integral Legendre function of the second kind. It is also referred

as Qm− 21 (k β) =

π (2 · k β)

m+ 21

2m

CR IP T

to as Q-function or a toroidal function of zeroth order which can be represented ∞ X (4n + 2m − 1)!! 1 2n (n + m)!(n)! (2 · k β)2n 2 n=0

(19)

where (4n + 2m − 1)!! = 1 · 3 · 5 · · · (4n + 2m − 1) for all m, n ≥ 0. With Eq. (18), Eq. (16) is written as, A(k r) =

µ0 Ms ˆ 0 φ 4π ×

∞ X

m=0

Z

h1

h2

Z

2π

π

0

1 p

k ρR

AN US

k

k

k

0

0

m Qm− 12 ( β) cos[m( φ − φ )]Rdφ dz

(20) 0

ED

M

ˆ 0 = − sin(φ0 )k x ˆ . Equation (20) is further written as ˆ + cos(φ0 )k y where φ s Z h1 ∞ µ0 Ms R X k k m Qm− 21 (k β)dz 0 A( r) = kρ 4π 2 h 2 m=0 Z 2π ˆ + cos(φ0 )k y ˆ ] cos[m(k φ − φ0 )]dφ0 × [− sin(φ0 )k x 0 | {z }

(21)

kζ

ˆ . Substituting Eq. (19) where k ζ is not zero only when m = 1, and k ζ = π · kφ into Eq. (21), we obtain

A(k r) =

CE

PT

k

∞

µ0 Ms k ˆ X R2n+2 (4n + 1)!! φ 4 22n (n + 1)!n! n=0 Z h1 k 2n+1 ρ × dz 0 k 2 2 k 0 2 2n+ 23 h2 [ ρ + R + ( z − z ) ]

(22)

According to Eq. (4), the magnetic flux density of k th permanent magnet in its local cylindrical coordinate system is denoted as   (n) ∞ k k X ∂ k Aφ  k k  1 ∂ Az − ∂ Aφ   ˆ Bρ = kρˆ  ρ = −  kρ ∂kφ ∂kz ∂kz n=0    k k  ∞ k (n) k (n)  k X  A A ∂ ∂( ρ A ) ∂ A 1     φ ρ φ φ k   ˆ ˆ Bz = k z − k    k = kz   ∂kρ + kρ   ρ ∂kρ ∂ φ n=0  k  ∂ Aρ ∂ k Az  k ˆ  Bφ = kφ  k − k  =0 ∂ z ∂ ρ

AC

140

13

(23)

(24) (25)

ACCEPTED MANUSCRIPT

So far, we have developed the model for the magnetic field distribution of a permanent magnet in the cylindrical coordinate system by Eqs. (23)-(25). The model allows magnetic field analysis of the whole permanent magnet rotor array

145

CR IP T

in spherical coordinate system.

4. Magnetic field distribution of rotor structure in spherical coordinate system

The complete magnetic field distribution of the 3D magnetic array in global spherical coordinate is obtained by superposition of all poles that are trans-

150

AN US

formed from corresponding local cylindrical coordinate system, following the steps described presently.

Firstly, the magnetic field at point 0 P created by the k th permanent magnet

M

pole is calculated. The point 0 P in the global Cartesian coordinate system is     0 sin(0 θ) cos(0 φ) x         0 (26) r = 0 y  = 0 ρ  sin(0 θ) sin(0 φ)      0 0 cos( φ) z

ED

As shown in Fig. 3, 0 P is reformulated as k P in the local coordinate system

of the k th permanent magnet by rotation matrix. The rotation matrix k R0 is

PT

given by

k

R0 = Ry (0 θk )Rz (0 φk )

AC

CE

Then, the point k P is obtained   k x   k  k k r =  y  = R0 · 0 r = Ry (0 θk )Rz (0 φk ) · 0 r   k z     cos(0 θk ) 0 − sin(0 θk ) cos(0 φk ) sin(0 φk ) 0 0 x         = − sin(0 φk ) cos(0 φk ) 00 y  0 1 0     0 0 1 0z sin(0 θk ) 0 cos(0 θk )

(27)

14

(28)

ACCEPTED MANUSCRIPT

(29)

the current model in the k th local cylindrical coordinate system     P∞ ∂ k Aφ(n) k − n=0 ∂ k z B   ρ   k   k k   0 B( r) =  Bφ  =   k (n)   P  k (n)  A ∞  ∂ Aφ k  ∂ k ρ + k φρ   Bz n=0

(30)

CR IP T

Denoting the point k P in k th local cylindrical coordinate system,    p k k x2 + k y 2 ρ     k     k r = k φ =  tan−1 ( k xy )      k k z z

AN US

the magnetic field distribution of the point k P is obtained by Eqs. (23)-(25) of

By using the transformation matrix

C

Mc from cylindrical coordinate system

to Cartesian coordinate system, the magnetic field distribution of point k P is restructured

   cos(k φ) 0  k  kB    B ρ ρ    k  = B(k r) = k By  = C Mc (k φ)   sin(k φ) 0 k  k Bz    Bz k 0 1 Bz k

Bx





M





(31)

ED

The magnetic field distribution of the k th permanent magnet pole is now transformed from the k th local coordinate to the global Cartesian coordinate system 0

Bkx

    0 T 0 k Bk =0 Bky  = 0 Rk · k B = RT z ( φk ) · Ry ( θk ) B   0 Bkz     cos(0 φk ) − sin(0 φk ) 0 cos(0 θk ) 0 sin(0 θk ) k Bx     k    = sin(0 φk ) cos(0 φk ) 0  By  0 1 0     0 0 1 − sin(0 θk ) 0 cos(0 θk ) k Bz

AC

CE

0



PT



(32)

To obtain the magnetic field of the k th permanent magnet distributed in spher-

15

ACCEPTED MANUSCRIPT

CR IP T

ical components, 0 Bk is reformulated as     0 0 Bkx Bkρ         0 Bk (0 r) = 0 Bkθ  = s MC (0 ρ, 0 θ, 0 φ)0 Bky      0 0 Bkz Bkφ     sin(0 θ) cos(0 φ) sin(0 θ) sin(0 φ) cos(0 θ) 0 Bkx      0   =cos(0 θ) cos(0 φ) cos(0 θ) sin(0 φ) − sin(0 θ)  Bky      0 Bkz − sin(0 φ) cos(0 φ) 0

(33)

where s MC is the transformation matrix from the Cartesian coordinate system

AN US

to spherical coordinate system.

Finally, the magnetic field distribution at the point 0 P by all the permanent magnets in global spherical coordinate system is computed by superposition 155

principle 0

Np X

Bρ = 0ρˆ

0

Bkρ

(34)

0

Bkθ

(35)

M

k=1

ED

0

0

Bθ = 0θˆ

Np X

k=1

ˆ Bφ = 0φ

Np X

0

Bkφ

(36)

k=1

Above all, the complete magnetic field distribution of the spherical actuator

PT

is obtained as a system of equations (34)-(36), which can be used for further design optimization of the spherical actuator and theoretical analysis of the ac-

CE

tuator performance. A flowchart is illustrated in Fig. 5 to summarize all steps

160

of calculation.

AC

5. Numerical validation of magnetic field and structure parameter analysis

5.1. Numerical model with finite element method Finite element approach is employed to validate the analytical magnetic

165

model. As the azimuthal flux component is almost equal to zero, this implies that 16

ACCEPTED MANUSCRIPT

Start

k=1 Yes

AN US

k<=16

CR IP T

Magnetic field 0B( 0Br,0Bf,0Bq ) of point 0 0 0 0 P( x, y, z) in spherical coordinate system

k

P(kx,ky,kz) = kR0 0P

Current model B( kBr,kBf,kBz )

M

k

C

ED

B( kBx,kBy,kBz )

Mc

k

Bk( kBkx,kBky,kBkz ) = 0Rk kB

AC

CE

PT

0

0

s

Bk( 0Bkr,0Bkf,0Bkq ) 0

MC

B=Ʃ0Bk k=k+1

No

End Figure 5: Flowchart of 3D magnetic field calculation

17

ACCEPTED MANUSCRIPT

only axial and radial flux components need to be analyzed. A single permanent magnet pole is modeled using the software Maxwell Ansoft, and the air-gap magnetic field distribution of a single permanent magnet is obtained. For both

170

CR IP T

analytical and numerical results, the magnetic field distribution is generated for a surface A as presented in Fig. 4.

Figures 6 and 7 present the axial and radial components of magnetic flux density along yk -direction, on the lines of 0.5 mm (zk = 50.5 mm) and 1.5 mm

(zk = 51.5 mm) away from the upper surface of a single permanent magnet

pole in the local coordinate system, respectively. Both analytical model and the numerical result are presented for comparison.

AN US

175

From Figs. 6(a) and 7(a), it’s seen that the value of the axial flux component is relatively uniform and high at the permanent magnet’s circular surface (−10 mm< yk < 10 mm). It decreases significantly close to the edge of the permanent magnet (yk = ±10 mm). As the distance from the permanent magnet 180

surface becomes larger (yk > 10 mm and yk < −10 mm), the slope at the edge

M

of permanent magnet tends to be less steep. From Figs. 6(b) and 7(b), it is seen that the radial flux component is zero at the center of permanent magnet

ED

pole (yk = 0 mm) and gradually increases to the highest value at about the permanent magnet fringe (yk = ±10 mm), and then decreases again. Comparison of the results shows that the developed analytical model fits

PT

well with the finite element results. We further evaluate the performance of the analytical model using two methods. Firstly, the normalized root-mean-square

CE

deviation (NRMSD) is applied to estimate the accuracy of the model, which is

AC

defined as

185

1 NRMSD = |umax − umin |

s

PN

j=1 (uj

N

− vj )

× 100%

(37)

where u is the value from the numerical model, v is the value of the analytical model, and N is the number of the data. The NRMSD values of the axial flux components from Figs. 6(a) and 7(a) are 7.63% and 5.94%, respectively. The NRMSD values of the radial flux components from Figs. 6(b) and 7(b) are 6.88% and 5.90%, respectively. This result shows that the analytical model is 18

ACCEPTED MANUSCRIPT

0.5 AM FEM

CR IP T

0.4

0.2

k

Bz (7)

0.3

0

-0.1 -25

-20

-15

-10

AN US

0.1

-5

0

5

10

15

20

25

yk (mm)

(a) Axial magnetic flux density k Bz along yk -direction

0.25

B (7)

ED

0.2

k

AM FEM

M

0.3

0.15

PT

0.1

AC

CE

0.05

Figure 6:

0 -25

-20

-15

-10

-5

0

5

10

15

20

25

yk (mm)

(b) Radial magnetic flux density k Bρ along yk -direction

Comparison of the magnetic field distribution along yk -direction on the line 0.5

mm away from the upper surface of a permanent magnet pole in local cylindrical coordinate (zk = 50.5 mm) between the analytical method (AM) and finite element method (FEM)

19

ACCEPTED MANUSCRIPT

0.45 AM FEM

CR IP T

0.4 0.35 0.3

0.2

k

B z (T)

0.25

0.15 0.1

0 -0.05 -25

-20

-15

-10

AN US

0.05

-5

0

5

10

15

20

25

yk (mm)

(a) Axial magnetic flux density k Bz along yk -direction

0.25

B (T)

ED

0.2

k

AM FEM

M

0.3

0.15

PT

0.1

AC

CE

0.05

Figure 7:

0 -25

-20

-15

-10

-5

0

5

10

15

20

25

yk (mm)

(b) Radial magnetic flux density k Bρ along yk -direction

Comparison of the magnetic field distribution along yk -direction on the line 1.5

mm away from the upper surface of a permanent magnet pole in local cylindrical coordinate (zk = 51.5 mm) between the analytical method and finite element method

20

ACCEPTED MANUSCRIPT

190

effective with a reasonable accuracy less than 10%. Secondly, Pearson correlation coefficient (PCC) is introduced to assess the correlation between the analytical model and numerical model. The correlation

CR IP T

coefficient ranges from -1 to 1. As the value of PCC closes to 1 or -1, it implies that the correlation becomes higher between two models. The PCC is defined as

PN

− u)(vj − v) qP N 2 2 (u − u) j j=1 j=1 (vj − v)

PCC = qP N

j=1 (uj

(38)

where u and v are the mean values of the numerical and analytical results,

AN US

respectively. The PCC values of the axial flux components from Figs. 6(a) and

7(a) are 0.9864 and 0.9910, respectively. The PCC values of the radial flux components from Figs. 6(b) and 7(b) are 0.9877 and 0.9876, respectively. It is 195

seen that there is quite high correlation between the analytical and numerical models. Additionally, compared to the numerical model with a computing time of about 2.03 minutes, the running time of analytical model is only about 0.013

M

seconds, which greatly improves the efficiency of calculation. 5.2. Analysis of permanent magnet structure

ED

The analytical model is finally used to investigate the effect of permanent

200

magnet structure on the magnetic flux density. The radius and the height of

PT

permanent magnet are main parameters of the permanent magnet. By varying the radius and height simultaneously, the magnetic field distribution of the permanent magnet is obtained by analytical model. Two representative points are chosen on the line 0.5 mm away from the surface of the permanent magnet

CE

205

as shown in Fig. 8. The axial magnetic flux density is comparatively high along

AC

the central axis of permanent magnet, therefore point k PA at the center is used

to analyze the effect of permanent magnet size on the axial flux component. The radial magnetic flux density reaches the maximum at the edge of the permanent

210

magnet, and then point k PB is selected to analyze the effect of permanent magnet size on the radial flux component.

21

ACCEPTED MANUSCRIPT

zk R k

PB

k

PA

CR IP T

d

AN US

h1

h2

yk

O

M

Figure 8: The position of the calculated points (d = 0.5 mm)

Figure 9 shows that when the permanent magnet length is fixed, both mag-

ED

netic flux components decrease with the increasing permanent magnet radius. However, as the permanent magnet length increases, the decreasing trend grad215

ually levels off. Conversely, when the permanent magnet radius is set, both flux

PT

components will increase with the increasing permanent magnet length. As the permanent magnet radius becomes larger, the change is more obvious.

CE

In summary, the increase of magnetic flux density is not always proportional to the length or the radius of the permanent magnet structure. The results from

220

the analytical model can be used as a reference to design the structure param-

AC

eters of the permanent magnet poles, which can avoid the waste of materials and unnecessary volume of the mechanism. Taking the cylindrical permanent

magnet design reported in [37] as example, the magnet has dimensions of radius R = 12.5 mm, height l = 20 mm and the volume being equal to 9.817 cm2 . By

225

looking the contour plot of Fig. 9, we can find that a design with parameters of radius R = 10 mm and height l = 16 mm can yield almost the same magnetic 22

CR IP T

ACCEPTED MANUSCRIPT

(12.5, 20, 0.4910)

AN US

(10, 16, 0.4870)

(a) Axial magnetic flux density Bz at the point

kP A

on the line 0.5

M

mm away from the surface of the permanent magnet

ED

(12.5, 20, 0.2496 )

AC

CE

PT

(10, 16, 0.2516)

Figure 9:

(b) Radial magnetic flux density Bρ at the point k PB on the line 0.5 mm away from the surface of the permanent magnet surface surface Magnetic flux density at the points 0.5 mm away from the permanent magnet

surface

23

ACCEPTED MANUSCRIPT

flux density, but the volume is reduced to 5.027 cm2 . By this way, the rotor structure can thus be improved in a more compact structure. The outer rotor

230

CR IP T

shell can correspondingly be reduced from 58.5 mm to 50 mm (Table 1). 6. Experiments

Upon our numerical validation of the analytical model of a single cylindri-

cal permanent magnet pole, experiments were conducted to verify the model. As the rotor array is symmetric in the spherical coordinate, only two adjacent

permanent magnet poles with opposite directions of magnetization on the same latitude are required to be measured. 6.1. Measurement platform

M

1

4

5

3

2

AC

CE

PT

ED

6

AN US

235

Figure 10: Experimental bench, (1) rotor, (2) markers, (3) angle indicator platform, (4) hall probe, (5) supporting platform, (6) Gauss meter

The experiments were conducted with a prototype of the permanent magnet spherical actuator as shown in Fig. 1(a), with specification presented in Table 1. 24

ACCEPTED MANUSCRIPT

In order to measure the air-gap magnetic field distribution of the 3D permanent 240

magnet array, a measuring platform was designed as shown in Fig. 10. The rotor is installed on an angle indicator platform so that it can rotate around

CR IP T

the platform by a certain angle. A series of points are marked on the measuring part of rotor surface. A hall probe is fixed by a supporting platform above the

measuring points. The probe is connected to a Gauss meter. The Gauss meter 245

(TM-801) is produced by KANETEC with a measuring range of 0-3000 mT and a resolution of 0.01 mT. The measured magnetic flux density can be directly

AN US

read and displayed by the computer. 6.2. Experimental results

The results of 3D magnetic field distribution in spherical coordinate system 250

are shown in Figs. 11 - 13, displayed together with analytical and numerical results for comparison. Good agreement of three results can be observed. From Fig. 11, it can be seen that the results have two peaks consistent with

M

the permanent magnet arrangement at the same latitude, and Bρ reaches the peak value at the magnetization axis of the permanent magnet pole in the global rotor frame. When the measured points are along the longitudinal direction

ED

255

passing by the center of the permanent magnet pole, i.e., φ =

3π 8

or

5π 8

with

varying θ, Bφ is equal to zero. As the measured points become far from the

PT

permanent magnet center, the value of Bφ gradually grows larger or smaller along the latitudinal direction, and reaches its peak value close to the edge of 260

the permanent magnet pole as shown in Fig. 12. Similarly, when the measured

CE

points are along the latitudinal direction passing by the center of the permanent magnet pole, i.e., θ =

5π 12

with varying φ, Bθ is always equal to zero. As the

AC

measured points become far from the permanent magnet center, the value of Bθ gets larger or smaller along the longitudinal direction, and reaches its peak

265

value close to the edge of the permanent magnet pole as shown in Fig. 13. In terms of efficiency, the numerical model takes more than 10 minutes for

every calculation, whereas the analytical model only takes about 0.35 seconds for the same results, which indicates a high efficiency with the analytical model 25

ED

M

AN US

(a) Analytical results

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

(b) Simulation results

(c) Experimental results Figure 11: Comparison of the radial flux component Bρ among analytical method, simulation model and experimental results at ρ = 50.5 mm

26

ED

M

AN US

(a) Analytical results

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

(b) Simulation results

(c) Experimental results Figure 12:

Comparison of the azimuthal flux component Bφ among analytical method,

simulation model and experimental results at ρ = 50.5 mm

27

ED

M

AN US

(a) Analytical results

CR IP T

ACCEPTED MANUSCRIPT

AC

CE

PT

(b) Simulation results

(c) Experimental results Figure 13: Comparison of the zenith flux component Bθ among analytical method, simulation model and experimental results at ρ = 50.5 mm

28

ACCEPTED MANUSCRIPT

developed.

270

7. Conclusion

CR IP T

This paper presents an analytical model to formulate the complex 3D magnetic field distribution of the spherical actuators in global spherical coordinate.

The model is developed for a new permanent magnet spherical actuator designed with a two-layer permanent magnet rotor array and 3-layer stator coils which 275

potentially improve the system performance with higher torque output. The model of magnetic field is developed analytically by combining of the current

AN US

model by toroidal functions and transformation method. Both simulations and experiments are carried out to show the validity of the analytical model.

A contribution of this paper is the development of the analytical magnetic 280

model which integrates the current model by toroidal functions and transformation method. The proposed method is characterized by efficient computation

M

with reasonable accuracy and versatility for cylindrical permanent magnet poles in any magnetostatic field. Using the analytical magnetic model, design parameters of the permanent magnet can be optimized. A case of design is included, which leads to a more compact structure of the spherical actuator.

ED

285

In future work, a more systematic analysis of the spherical actuator will be

PT

carried out. The analytical model will be used for in-depth design analysis and further investigation of the analytical torque model, which can be utilized to reveal the influence of the electrical parameters and the running conditions of the spherical actuator. Study of real-time control implementation is also a topic

CE

290

for future research.

AC

Acknowledgment This work is partially supported by National Natural Science Foundation of

China under Grant No. 51475033.

29

ACCEPTED MANUSCRIPT

295

References [1] K.-M. Lee, C.-K. Kwan, Design concept development of a spherical stepper

7 (1) (1991) 175–181. doi:10.1109/70.68082.

CR IP T

for robotic applications, IEEE Transactions on Robotics and Automation

[2] J. Wang, G. W. Jewell, D. Howe, A novel spherical actuator: design and control, IEEE Transactions on Magnetics 33 (5) (1997) 4209–4211. doi:

300

10.1109/20.619712.

[3] S. Bai, X. Li, J. Angeles, A review of spherical motion generation using

AN US

either spherical parallel manipulators or spherical motors, Mechanism and

Machine Theory 140 (2019) 377 – 388. doi:https://doi.org/10.1016/ j.mechmachtheory.2019.06.012.

305

[4] C. M. Gosselin, J.-F. Hamel, The agile eye: a high-performance threedegree-of-freedom camera-orienting device, in: Proceedings of the 1994

M

IEEE International Conference on Robotics and Automation, Vol. 1, 1994, pp. 781–786. doi:10.1109/ROBOT.1994.351393. [5] S. Bai, Optimum design of spherical parallel manipulators for a prescribed

ED

310

workspace, Mechanism and Machine Theory 45 (2) (2010) 200 – 211. doi:

PT

https://doi.org/10.1016/j.mechmachtheory.2009.06.007. [6] D. Kang, J. Lee, Analysis of electric machine charateristics for robot eyes using analytical electromagnetic field computation method, IEEE Trans-

CE

actions on Magnetics 50 (2) (2014) 785–788. doi:10.1109/TMAG.2013.

315

2278937.

AC

[7] S. Ahmadi, J. S. Moghani, M. Mirsalim, Simulation and analysis of a

320

novel PM spherical 3-DOF actuator with E-shaped stator and bladeshaped rotor structure, in: Power Electronics, Drives Systems and Technologies Conference (PEDSTC), 2018 9th Annual, IEEE, 2018, pp. 59–64. doi:10.1109/PEDSTC.2018.8343772.

30

ACCEPTED MANUSCRIPT

[8] J. Chabot, H. Schaub, Spherical magnetic dipole actuator for spacecraft attitude control, Journal of Guidance, Control, and Dynamics (2016) 911– 915doi:10.2514/1.G001471. [9] D.-K. Kim, H. Yoon, W.-Y. Kang, Y.-B. Kim, H.-T. Choi, Development

CR IP T

325

of a spherical reaction wheel actuator using electromagnetic induction,

Aerospace Science and Technology 39 (2014) 86–94. doi:10.1016/j.ast. 2014.09.004.

[10] J.-S. Lee, D.-K. Kim, S.-W. Baek, S.-H. Rhyu, B.-I. Kwon, Newly struc-

AN US

tured double excited two-degree-of-freedom motor for security camera,

330

IEEE Transactions on Magnetics 44 (11) (2008) 4041–4044. doi:10.1109/ TMAG.2008.2002802.

[11] S. Toyama, A. Kobayashi, Development of spherical ultrasonic motor, CIRP Annals-Manufacturing Technology 45 (1) (1996) 27–30. doi:10. 1016/S0007-8506(07)63010-8.

M

335

[12] J. F. P. Fernandes, V. M. Machado, P. J. C. Branco, Magnetic field analysis

ED

in shell-like spherical induction machines with zenithal traveling waves, IEEE Transactions on Energy Conversion 32 (3) (2017) 1081–1089. doi: 10.1109/TEC.2017.2696164. [13] J. F. P. Fernandes, S. M. Vieira, P. J. C. Branco, Multiobjective op-

PT

340

timization of a shell-like induction spherical motor for a power-assisted

CE

wheelchair, IEEE Transactions on Energy Conversion 33 (2) (2018) 660– 669. doi:10.1109/TEC.2017.2761983.

AC

[14] W. Chen, L. Zhang, L. Yan, J. Liu, Design and control of a three degree-

345

of-freedom permanent magnet spherical actuator, Sensors and Actuators A: Physical 180 (2012) 75 – 86. doi:https://doi.org/10.1016/j.sna. 2012.04.010.

[15] X. Chen, J. Hu, K. Chen, Z. Peng, Modeling of electromagnetic torque considering saturation and magnetic field harmonics in permanent magnet 31

ACCEPTED MANUSCRIPT

synchronous motor for hev, Simulation Modelling Practice and Theory 66

350

(2016) 212 – 225. doi:https://doi.org/10.1016/j.simpat.2016.02. 012.

CR IP T

[16] B. Li, G. D. Li, H. F. Li, Magnetic field analysis of 3-DOF permanent

magnetic spherical motor using magnetic equivalent circuit method, IEEE Transactions on Magnetics 47 (8) (2011) 2127–2133. doi:10.1109/TMAG.

355

2011.2123102.

[17] K.-M. Lee, H. Son, Distributed multipole model for design of permanent-

AN US

magnet-based actuators, IEEE Transactions on Magnetics 43 (10) (2007) 3904–3913. doi:10.1109/tmag.2007.904709. 360

[18] H. Son, K.-M. Lee, Two-DOF magnetic orientation sensor using distributed multipole models for spherical wheel motor, Mechatronics 21 (1) (2011) 156–165. doi:10.1016/j.mechatronics.2010.10.001.

M

[19] L. Yan, I.-M. Chen, G. Yang, K.-M. Lee, Analytical and experimental investigation on the magnetic field and torque of a permanent magnet spherical actuator, IEEE/ASME Transactions on Mechatronics 11 (4) (2006) 409–

365

ED

419. doi:10.1109/TMECH.2006.878545. [20] C. Xia, J. Xin, H. Li, T. Shi, Design and analysis of a variable arc perma-

PT

nent magnet array for spherical motor, IEEE Transactions on Magnetics 49 (4) (2013) 1470–1478. doi:10.1109/TMAG.2012.2231092. [21] Z. Li, L. Zhang, W. He, Y. Zhang, Q. Wang, Magnetic field analysis of

CE

370

hybrid driven permanent magnet multi-dof motor, 2017 20th International

AC

Conference on Electrical Machines and Systems (ICEMS) (2017) 1–6.

[22] L. Yan, F. Liang, Z. Jiao, T. Wang, Magnetic field analysis of novel spher-

375

ical actuators with three-dimensional pole arrays, Review of Scientific Instruments 87 (6) (2016) 33–175. doi:10.1063/1.4953920.

[23] L. Yan, I.-M. Chen, C. K. Lim, G. Yang, W. Lin, K.-M. Lee, Hybrid torque modeling of spherical actuators with cylindrical-shaped magnet 32

ACCEPTED MANUSCRIPT

poles, Mechatronics 21 (1) (2011) 85–91. doi:10.1016/j.mechatronics. 2010.08.009. 380

[24] H. Y. Kim, H. C. Kim, D. G. Gweon, Magnetic field analysis of a VCM

doi:10.1016/j.sna.2013.02.024.

CR IP T

spherical actuator, Sensors & Actuators A Physical 195 (6) (2013) 38–49.

[25] X. Li, J. Liu, W. Chen, S. Bai, Integrated design, modeling and analysis

of a novel spherical motion generator driven by electromagnetic principle, Robotics and Autonomous Systems 106 (2018) 69 – 81. doi:https://doi.

385

AN US

org/10.1016/j.robot.2018.04.006.

[26] X. Li, S. Bai, W. Chen, J. Liu, Torque modelling and current optimization of a spherical actuator built as an electro-magnets driven spherical parallel manipulator, in: 2017 IEEE International Conference on Cybernetics and Intelligent Systems (CIS) and IEEE Conference on Robotics, Automation

390

M

and Mechatronics (RAM), 2017, pp. 626–631. doi:10.1109/ICCIS.2017. 8274850.

ED

[27] J. Li, E. S. Barjuei, G. Ciuti, Y. Hao, P. Zhang, A. Menciassi, Q. Huang, P. Dario, Magnetically-driven medical robots: An analytical magnetic model for endoscopic capsules design, Journal of Magnetism and Magnetic

395

PT

Materials 452 (2018) 278 – 287. doi:https://doi.org/10.1016/j.jmmm. 2017.12.085.

CE

[28] L. Yan, Z. Wu, Z. Jiao, C. Y. Chen, I.-M. Chen, Equivalent energized coil model for magnetic field of permanent-magnet spherical actuators, Sensors & Actuators A Physical 229 (2015) 68–76. doi:10.1016/j.sna.2015.03.

AC

400

016.

[29] T. Taniguchi, An analytical computation of magnetic field generated from a cylinder ferromagnet, Journal of Magnetism and Magnetic Materials 452 (2018) 464 – 472. doi:https://doi.org/10.1016/j.jmmm.2017.11.078.

33

ACCEPTED MANUSCRIPT

405

[30] A. Caciagli, R. J. Baars, A. P. Philipse, B. W. Kuipers, Exact expression for the magnetic field of a finite cylinder with arbitrary uniform magnetization, Journal of Magnetism and Magnetic Materials 456 (2018) 423 – 432. doi:

CR IP T

https://doi.org/10.1016/j.jmmm.2018.02.003. [31] Z. Qian, Q. Wang, G. Li, X. Guo, C. Hu, H. Yan, Design and analysis

of permanent magnetic spherical motor with cylindrical poles, in: 2013

410

International Conference on Electrical Machines and Systems (ICEMS), 2013, pp. 644–649. doi:10.1109/ICEMS.2013.6754492.

AN US

[32] J. P. Selvaggi, S. Salon, O. M. Kwon, M. V. K. Chari, Computation of

the three-dimensional magnetic field from solid permanent-magnet bipolar cylinders by employing toroidal harmonics, IEEE Transactions on Magnet-

415

ics 43 (10) (2007) 3833–3839. doi:10.1109/TMAG.2007.902995. [33] Q. Wang, Z. Qian, G. Li, Vision based orientation detection method and

M

control of a spherical motor, in: 2010 53rd IEEE International Midwest Symposium on Circuits and Systems, 2010, pp. 1145–1148. doi:10.1109/ MWSCAS.2010.5548861.

ED

420

[34] L. Zhang, W. Chen, J. Liu, X. Wu, Analysis and decoupling control of a permanent magnet spherical actuator, Review of Scientific Instruments

PT

84 (12) (2013) 125001. doi:10.1063/1.4833681. [35] J. Liu, H. Deng, C. Hu, Z. Hua, W. Chen, Adaptive backstepping sliding mode control for 3-DOF permanent magnet spherical actuator, Aerospace

CE

425

Science and Technology 67 (2017) 62 – 71. doi:https://doi.org/10.

AC

1016/j.ast.2017.03.032.

[36] A. Kumar, S. Marwaha, A. Marwaha, N. Kalsi, Magnetic field analysis

430

of induction motor for optimal cooling duct design, Simulation Modelling Practice and Theory 18 (2) (2010) 157 – 164. doi:https://doi.org/10. 1016/j.simpat.2009.10.002.

34

ACCEPTED MANUSCRIPT

[37] L. Zhang, W. Chen, J. Liu, C. Wen, A robust adaptive iterative learning control for trajectory tracking of permanent-magnet spherical actuator, IEEE Transactions on Industrial Electronics 63 (1) (2016) 291–301. doi: 10.1109/TIE.2015.2464186.

AC

CE

PT

ED

M

AN US

CR IP T

435

35