Hang Zhao

Chunhua Liu

Zaixin Song

Senyi Liu

Thierry Lubin

Analytical Model for Magnetic-Geared Double-Rotor Machines and Its d-q Axis Determination

This paper presents an analytical subdomain model for the prediction of various electromagnetic parameters of magnetic-geared double-rotor machine (MGDRM). By dividing the MGDRM into six subdomains and solving the Maxwell equations in polar coordinates for each region, the vector magnetic potential distribution can be derived. Subsequently, the flux density distribution, back electromotive force (EMF), and output torque, can be calculated. Furthermore, based on the proposed model, the equivalent d-q axis is elaborated. Also, the corresponding d and q inductance can be obtained. It proves that the MGDRM can be regarded as a non-salient-pole synchronous machine. Besides, the power factor of the MGDRM under i d =0 control is deduced and optimized for the different thicknesses of modulators. In addition, the demagnetization capability is analysed by adopting the subdomain model. Finally, the accuracy of proposed subdomain model for MGDRMs is validated via the finite-element method (FEM).

Nomenclature

A z Vector magnetic potential. B r , B θ Radial and tangential component of magnetic flux density. M , J Magnetization vector of magnets, Current density vector. μ 0 Magnetic permeability of vacuum. P, Q Number of stator slots and modulator pieces. P i Pole pair Number of inner rotor permanent magnets (PMs). β, γ, δ Modulator opening angle, stator slot width angle and slot opening angle. φ 0 , θ 0 Initial circumferential position of the inner rotor and modulator. θ i , θ j , J l Circumferential positions of PMs, i th modulator opening and j th slot opening, current density in the l th slot. R 1 , R 2 , R 3 , R 4 Radii of inner rotor yoke, inner rotor outer surface, modulator inner surface, modulator outer surface. R 5 , R 6 , R 7 Radii of stator bore, slot top, slot bottom. L Axial length of the machine. D Root mean square (RMS) current density in windings. B r , B' r PM remanence at normal temperature and at 120 °C. N, K, T, M Harmonic order used in air gap and PM subdomain, modulator subdomain, slot opening subdomain and slot subdomain. φ l , ψ, ψ d , ψ q Flux linkages of one coil side, certain phase, d axis and q axis. L xx , M xy Self inductance of phase x, mutual inductance between phase x and y. L d , L q Inductance of d axis and q axis. E, f Back electromagnetic force (EMF) and current frequency of windings. α, φ Electric angle of the equivalent d axis, power factor angle of the machine. Ω in , Ω out , Ω d Rotating speed of inner rotor, modulator and equivalent d axis. T rot , T mod Electromagnetic torque of inner rotor and modulator.

Introduction

Due to the high efficiency, long range and smooth start, hybrid electric vehicles (HEVs) are taking up the more share in the whole vehicle market [START_REF] Emadi | Power Electronics and Motor Drives in Electric, Hybrid Electric, and Plug-in Hybrid Electric Vehicles[END_REF]. One of the key components in HEV is the power split device (PSD), which allows the internal combustion engine (ICE) and the electric motor to provide the drive power at different ratios [START_REF] Mashadi | Dual-Mode Power-Split Transmission for Hybrid Electric Vehicles[END_REF]. The existing PSDs in HEVs, such as Toyota Prius and Honda Insight, have adopted the gearbox to realize the power distribution, which inevitably decreases the efficiency and introduces the vibration and noise.

The magnetic-geared double-rotor machine (MGDRM) has a great potential for serving as the PSD in HEVs [START_REF] Liu | Emerging Electric Machines and Drivesan Overview[END_REF], [START_REF] Zhao | Analytical Modeling and Comparison of Two Consequent-Pole Magnetic-Geared Machines for Hybrid Electric Vehicles[END_REF], since it utilizes the electromagnetic force to transmit the power among different moving components. Thus, it can eliminate the disadvantages caused by the physical contact in mechanical gearboxes. The operating principle and structure requirements, such as the pole pair combination and winding configuration of MGDRMs, are investigated in [START_REF] Bai | Design and Analysis of a Magnetic-Field Modulated Brushless Double-Rotor Machine-Part I: Pole Pair Combination of Stator, Pm Rotor and Magnetic Blocks[END_REF]. Also, a step-skew modulator structure is introduced to solve the asymmetrical flux path problem, whose cogging torque can be significantly reduced [START_REF] Sun | Analysis of a Novel Magnetic-Geared Dual-Rotor Motor with Complementary Structure[END_REF]. In addition, the AlNiCo material and magnetizing windings are adopted to realize the gear-ratio variation in [START_REF] Liu | Electromagnetic Design of a New Electrically Controlled Magnetic Variable-Speed Gearing Machine[END_REF]. Besides, to enhance the flux-weakening ability of MGDRM, a consequent-pole and inserted PM structure are discussed in [START_REF] Zhao | Design of an Effective Double-Rotor Machine with Robust Mechanical Structure[END_REF]. Moreover, PMs are inserted into the gap of modulator to increase the torque density or interact with extra windings in [START_REF] Wang | Electrical-Continuously Variable Transmission System Based on Doubly Fed Flux-Bidirectional Modulation[END_REF], [START_REF] Ren | A Brushless Dual-Mechanical-Port Dual-Electrical-Port Machine with Spoke Array Magnets in Flux Modulator[END_REF].

Particularly, an accurate expression of the air-gap magnetic field distribution of MGDRMs is necessary to evaluate the performance and executing optimization at the initial design stage. The finite element method (FEM) can assist the designers to obtain the magnetic field distribution within a machine with the high accuracy. However, the mesh density of FEM should be high enough to acquire the electromagnetic parameters as accurately as possible, which further increases the computation time. Nevertheless, the numerical result in each mesh can be obtained, which would blurry the understanding of the fundamental principle of a certain machine. This is exactly true for MGDRMs, since it utilizes the magnetic-gearing effect. The magnetic field distribution is complex with the introduction of modulator [START_REF] Atallah | Design, Analysis and Realisation of a High-Performance Magnetic Gear[END_REF].

Another approach of machine design is the magnetic equivalent circuits (MECs), which have been developed for decades [START_REF] Severson | Magnetic Equivalent Circuit Modeling of the Ac Homopolar Machine for Flywheel Energy Storage[END_REF], [START_REF] Amrhein | 3-D Magnetic Equivalent Circuit Framework for Modeling Electromechanical Devices[END_REF]. However, the nodes in MEC must be selected delicately to reflect the flux paths within the machine. The motion of rotor makes the modeling even complicated [START_REF] Bash | Incorporating Motion in Mesh-Based Magnetic Equivalent Circuits[END_REF]. Since MGDRMs have two rotating components, MEC is unsuitable for the design and optimization. Besides, the analytical model based on conformal transformation is very useful in predicting the airgap magnetic flux distribution by taking into account of the slotting effect [START_REF] O'connell | A Schwarz-Christoffel-Based Analytical Method for Electric Machine Field Analysis[END_REF][START_REF] Boughrara | Magnetic Field Analysis of Inset and Surface-Mounted Permanent-Magnet Synchronous Motors Using Schwarz-Christoffel Transformation[END_REF][START_REF] Jalali | Analytical Modeling of Slotless Eccentric Surface-Mounted Pm Machines Using a Conformal Transformation[END_REF]. This method cannot be applied in analyzing the magnetic field of MGDRMs, because the modulation poles on the modulator cannot be regarded as equipotential.

The third approach of machine design is the subdomain model, which is an effective analytical method for solving the magnetic field of various electromagnetic devices. By dividing the machine into several subdomain regions, solving Maxwell equations within each region, and applying Fourier series to satisfy the boundary conditions, the magnetic field distribution can be achieved [START_REF] Zhu | Instantaneous Magnetic Field Distribution in Brushless Permanent Magnet Dc Motors. I. Open-Circuit Field[END_REF], [START_REF] Gysen | General Formulation of the Electromagnetic Field Distribution in Machines and Devices Using Fourier Analysis[END_REF]. In fact, various machine topologies have been investigated by using the subdomain model method [START_REF] Rahideh | Two-Dimensional Analytical Investigation of the Parameters and the Effects of Magnetisation Patterns on the Performance of Coaxial Magnetic Gears[END_REF][START_REF] Zhu | An Accurate Subdomain Model for Magnetic Field Computation in Slotted Surface-Mounted Permanent-Magnet Machines[END_REF][START_REF] Boughrara | General Subdomain Model for Predicting Magnetic Field in Internal and External Rotor Multiphase Flux-Switching Machines Topologies[END_REF][START_REF] Oner | Analytical on-Load Subdomain Field Model of Permanent-Magnet Vernier Machines[END_REF][START_REF] Lubin | 2-D Exact Analytical Model for Surface-Mounted Permanent-Magnet Motors with Semi-Closed Slots[END_REF][START_REF] Shakibapour | 2d Analytical Model for Heteropolar Active Magnetic Bearings Considering Eccentricity[END_REF]. In fact, the calculation time of FEM increases immensely with the size of electric machine. Yet, the calculation time of subdomain model is only related to the number of harmonic orders for a given electric machine. Therefore, the subdomain model could offer the convenience to engineers at the initial stage of motor design.

Although the analytical model of magnetic-geared machines with the single outer rotor is given in [START_REF] Zhang | Analytical Investigation on the Power Factor of a Flux-Modulated Permanent-Magnet Synchronous Machine[END_REF], the procedure of derivation is not given in detail. In particular, the rotation of modulator is not discussed. To the authors' knowledge, the comprehensive characteristics, operating principle and performances of MGDRMs have not been well investigated with an accurate analytical model in previous studies.

This paper proposes a subdomain model for MGDRMs, which can effectively predict the magnetic field distribution and secondary parameters in high accuracy. First, the modeling process of MGDRM is developed in Section II. Then, the calculated magnetic field distribution, back EMF, and torque-angle relation are verified by FEM. In addition, the equivalent d-q axis of MGDRM is initially acquired by using the subdomain model, and i d =0 control is applied to get the maximal output torque. Moreover, the power factor and demagnetization capability of MGDRM are also analyzed.

Subdomain Model

Fig. 2. The symbol and definitions of subdomains

A fourteen-pole/twenty-four-slot MGDRM with nine ferromagnetic bars and single-layer winding is investigated in this paper to verify the proposed method. It should be noted that the proposed model is also suitable for doublelayer-winding machine by changing the expression of current density in each slot. Since there are currents flowing through the windings, vector magnetic potential is adopted to calculate the magnetic field distribution within the MGDRM. The geometry of the investigated MGDRM is shown in Fig. 1: the modulator is composed of equally distributed iron bars on circumferential direction, serving as the outer rotor of the MGDRM. φ 0 , θ 0 are the initial circumferential position of the inner rotor and modulator, respectively, they are defined as zero where the positions of inner rotor and modulator are as shown in Fig. 1. θ i , θ j can be further given by: 0 2 + with 1

2 i i iQ Q   = -+   (1)
2 with 1 2

j j jP Q   = -+   (2)
Assumptions are made to simplify the calculation as:

(1) The magnetic field on the axial direction is regarded as zero, and the end effect is ignored [START_REF] Sun | Analysis and Control of Complementary Magnetic-Geared Dual-Rotor Motor[END_REF];

(2) The permeability of soft-magnetic material is assumed to be infinite. Thus, the magnetic field strength direction should be vertical to the soft-magnetic boundary, namely it does not have tangential direction value. Also, the magnetic field strength is zero within the iron;

(3) The relative permeability of permanent magnets (PMs) and coil is assumed to be 1;

(4) The geometrical structure side of MGDRM is parallel to either θ or r direction in polar coordinates.

As shown in Fig. 2, the MGDRM can be divided into several subdomains: PM (Region I), inner air gap (Region II), modulator opening (Region III), outer air gap (Region IV), slot opening (Region V) and stator slot (Region VI). The connections among different subdomains can be expressed as boundary conditions. The magnetostatic Maxwell equation is satisfied in each subdomain, which can be further reduced to Poisson's equation or Laplace's equations under aforementioned assumptions. Then, by utilizing the method of separation of variables in polar coordinates, the general expression of vector magnetic potential in each subdomain can be obtained. Finally, the Fourier coefficients within these expressions can be achieved by solving the boundary condition equations. To simplify the equation form, the following notations are adopted in this paper: ( , ) and ( , )

zz zz zz x y x y P x y E x y y x y x         = + = -                
(3)

Boundary Condition

In polar coordinates, the flux density can be decomposed into the radial and tangential component B r , B θ . Their relations with the vector magnetic potential are as:

For simplification, the following boundary conditions are all expressed by using the vector magnetic potential.

Actually, the boundary conditions are inserted at the boundaries of different regions, which can be divided into three forms for the proposed MGDRM.

(1) Periodic condition: In tangential direction, the magnetic field distribution has a period of 2π, thus:

( , ) ( , 2 ) ij A r A r    =+ (5)
(2) Continuous condition: The vector magnetic potential should be continuous at the interface between two regions:

( , ) ( , ) Boundary is radial, is a constant ( , ) ( , ) Boundary is tangential, is a constant i j i j A r A r A R A R R   =      =   (6)
(3) Neumann condition: Generally, the tangential component of the magnetic field intensity H should be continuous at the interface. Therefore, air region and coil region have:

Boundary is in radial direction Boundary is in tangential direction j i j i rR rR A A A A rr    = = = =   =        =    (7)
I rR A r =  =  (8) [0,2 ]   IV, V 5 5
()

0 j IV rR rR A A r h r  = =     ==     (9) [ , ] jj elsewhere     +    I, II 22 ( , ) ( , )
I II A R A R  = (10) [0,2 ]   IV, V 55 ( , ) ( , ) j IV A R A R  = (11)
[0,2 ]  

I, II 22 II I r R r R AA rr ==  =  (12) [0,2 ]   V 0 0 j j j j A A       = =+  =       =    (13) 56 [ , ] r R R  II, III 3 3 () 0 i II rR rR A A r f r  = =     ==     (14)
[ , ]

ii elsewhere

    +    V, VI 66 ( , ) ( , ) jl A R A R  = (15) [0,2 ]   II, III 33 ( , ) ( , ) i II A R A R  = (16) [0,2 ]   V, VI 6 6
()

0 j l rR rR A A r s r  = =     ==     (17) [ , ] jj elsewhere     +    III 0 0 i i i i A A       = =+  =       =    (18) 23 [ , ] r R R  VI 1 () 2 1 () 2 0 ( ) 0 j j l l A lj A           = + - = + +  =    =    =    (19) 67 [ , ] r R R  III, IV 44 ( , ) ( , ) i IV A R A R  = (20) [0,2 ]   VI 7 0 l rR A r =  =  (21) 1 [ ( ), 2 1 ( )] 2 j j         + - ++ III, IV 4 4 () 0 i IV rR rR A A r g r  = =     ==     (22)
[ , ]

ii elsewhere

    +   
Since H=0 within the soft-magnetic material, then the boundary condition for the interface with soft-magnetic material can be expressed as:

0 Boundary is in radial direction 0 Boundary is in tangential direction i i rR A A r   = =  =       =    (23)
By adopting the different boundary condition forms in the different regions of MGDRM, Table 1 can be obtained.

Besides, from [START_REF] Lubin | Analytical Computation of the Magnetic Field Distribution in a Magnetic Gear[END_REF], the Region III and Region V have Neumann boundary conditions. Extra constraints should be exerted to guarantee the existence of solution [START_REF] Farlow | Partial Differential Equations for Scientists and Engineers[END_REF]. 00 11 0 and 0

Q P ij ij BB == ==  (24)
By using (68) to replace 0 j B with Jl, one can get:

1 0 P l l J = =  (25)
where ( 25) can be satisfied automatically, if the windings are energized with three-phase symmetrical currents.

General Solution of Poisson's and Laplace's

Equations Generally, if Coulomb gauge is adopted, the vector magnetic potential for every region should satisfy:

2 00 () A M J   = - - (26)
In regions where M and J are zero, ( 26) is reduced to Laplace's equation: Under polar coordinate, ( 28), ( 29), [START_REF] Bai | Characteristic Analysis and Verification of the Magnetic-Field-Modulated Brushless Double-Rotor Machine[END_REF] can be written as:

2 2 2 11 0 with =II, III, IV, V i i i A A A i r r r r     + + =    (31) 0 2 2 2
11

I I I A A A M r r r r r       + + =     (32) 0 2 2 2
11 with 1

j j j j A A A J j P r r r r      + + = -      (33)
By taking into account the boundary conditions in Table 1, the expressions of vector magnetic potential can be obtained, as shown in Table 2. Besides, the detailed Fourier coefficients are given in Appendix.

Finally, by solving the system of first-order linear equation, the Fourier coefficients of each expression can be obtained. So, the exact expression of vector magnetic potential in each subdomain is acquired.

Flux Density and Torque Calculation

By utilizing (39) and ( 4), the radial and tangential components of flux density in outer air gap can be derived as: ( , ) ( , ) ( ( , ) ( , )

( , ) )sin( ) ( , ) ( , ) ( , ) ( )cos( ) ( , ) ( , )
IV n IVr n n n IV n n n IV IV nn nn n nn P r R R B r A r E R R R r R Bn r E R R P r R R P r R R C D n r E R R r E R R     =  = =- + ++   (34)
( , ) ) cos( ) ( , ) ( , ) ( ( , )
)sin( ) ( , ) ( , ) ( , ) 34) and (35) [START_REF] Lubin | Analytical Computation of the Magnetic Field Distribution in a Magnetic Gear[END_REF]. According to the Maxwell stress tensor expression, the electromagnetic torque of inner air gap T ein and outer air gap T eout can be expressed as:

IV n IV n n n IV n n IV IV nn nn n nn E r R R B r A r E R R R E r R Bn r E R R E r R R E r R R C D n r E R R r E R R      =  = =- + -+   (35)
I I I nn I n n i n n i nn nn P r R P r R A r A A X r n n C X r n n P R R P R R       == = + + + +  (36) II 3 3 2 3 3 2 2 2 0 1 1 2 3 3 2 2 3 3 2 ( , ) ( , ) ( , ) ( , ) ( , ) ( )cos( ) ( )sin( ) ( , ) ( , ) ( , ) ( , )
II II II II II n n n n II n n n n n n n n n n P r R R P r R P r R R P r R R R A r A A B n C D n n E R R n E R R n E R R n E R R      = = = + + + +   (37) III / 4 / 3 00 1 / 3 4 / 3 4 ( , ) ( , ) ( , ) ln ( )cos( ( )) ( , ) ( ,
)
kk i i i i i k k i k kk E r R E r R k A r A B r A B E R R E R R               = = + + - -  (38)
( , ) ( , ) ( , ) ( , ) ( )cos( ) ( )sin( ) ( , ) ( , ) ( , ) ( , )
IV IV IV IV IV n n n n IV n n n n n n n n n n P r R R P r R P r R R P r R R R A r A A B n C D n n E R R n E R R n E R R n E R R      = = = + + + +   (39)
j j i i tt j t t j t tt E r R E r R t A r A B r A B E R R E R R               = = + + - -  (40) VI /7 22 6 0 0 7 1 / 6 7 ( , ) 1 1 1 ( , ) ( ln ) cos( ( ( ))) 2 2 ( , ) 2
m ll l l m j m m P r R R m A r A J R r r A m E R R             = = + - + -- -  (
2 2 0 0 ( , ) ( , ) i ein IIr i II i LR T B R B R d       =  (42) 2 2 0 0 ( , ) ( , ) o eout IVr o IV o LR T B R B R d       =  (43)
where R i , R o are the average radii of inner air gap and outer air gap, respectively. R i , R o are given by:

2 3 4 5 , 2 2 i o R R R R R R + + = = (44)
For the inner rotor, only T ein is exerted to drive it. However, the modulator is driven by T ein and T eout simultaneously. Hence, the resultant electromagnetic torque on the inner rotor T rot and modulator T mod can be expressed as:

, rot ein mod eout ein T T T T T = = - (45)

Flux Linkage and Back EMF Calculation

The flux linkage of each coil can be calculated as:

7 6 1 () 2 1 () 2 j j R turn ll R slot N L A rd dr A        ++ +- =  (46)
where A slot is the cross-sectional area of each slot as given by:

22 76 () 2 slot A R R  =- (47)
The coil connection of each phase is given by a matrix as: 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 0 1 1 0 0 1 1 0 0 0 0 1 1 0 0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0 0 1 1 0 0 0 0 1 1 0 0 1 1 0 0 0 0 1 1 0 0 0

C --   = --   --  --   --   --  (48)
Therefore, the flux linkage of each phase can be derived by:   1 2 24 , ,... ,...

A T B l C C          ==    ψ ( 49
)
The three-phase back EMF can be calculated as:

,,

i in i out i i A tt B C E d E dt t E         +  +   -  = - =     (50)

Analytical Prediction and FEM Validation

To validate the effectiveness of the proposed subdomain model, the analytical prediction is compared with the 2-D FEM result and considered the saturation effect of softmagnetic material from the tool of JMAG Designer. Table 3 gives the parameters of the investigated 14-pole/24-slot MGDRM. The FEM model has 80952 elements, 50581 nodes; the element size near the air gap is set as 0.5 mm to maintain the calculation accuracy. The calculation time of one step by FEM is 16.98 s, while it only takes 7.65 s to run a single step by the subdomain model. So, the proposed method could help designers to save the time at primary stage of motor design. 34567show the comparisons between the analytical prediction and FEM calculation results, where φ 0 =0, θ 0 =0. Figure 3 and4 show the no-load the radial and tangential flux density distribution for inner air gap and outer air gap; Figure 5 and6 show the radial and tangential flux density distribution for inner air gap and outer air gap when the windings are electrified with AC current. It can be seen that the outer air-gap flux density distribution changes greatly after the windings are electrified, while that of inner air gap changes slightly. Figure 7 shows the back-EMF under noload condition, where the rotating speed of inner rotor and modulator are 600 r/min and 400/3 r/min, respectively. A good agreement between the analytical prediction given by the subdomain model and FEM simulation can be seen for all these figures. However, due to the magnetic saturation within the modulator pieces and slot teeth of MGDRM, a small difference of magnetic flux density can still be observed, as shown in Fig. 4 and Fig. 6. Besides, the magnetic saturation becomes severe as the current increases in the stator windings, since the magnetic flux density difference is larger in Fig. 6 compared to that in Fig. 4.

Fig. 3. Analytically predicted and FEM calculated results with no-load inner air-gap flux density distribution

Fig. 6. Analytically predicted and FEM calculated results with on-load outer air-gap flux density distribution

Fig. 7. Analytically predicted and FEM calculated results with three-phase back EMFs

Electromagnetic Performance Analysis

Equivalent d-q Axis Determination

When the MGDRM is electrified with three-phase AC currents, the magnetic field generated by stator windings interacts with that produced by PMs. So, the MGDRM can be regarded as a PMSM, when observed from the outer airgap. However, the equivalent d-q axis of MGDRM is much more difficult to obtain compared with that of conventional PMSM, because the magnetic field distribution is changed by the modulation effect of ferromagnetic bars.

By adopting the subdomain model of MGDRM, the equivalent d-q axis of MGDRM can be gained. For the 14pole MGDRM with nine ferromagnetic bars, the fundamental harmonic of windings is 2 nd . Therefore, the radial-direction magnetic field distribution of 2 nd harmonics on the middle of outer air gap can be obtained by substituting n=2 into (34): _ ( , ) sin( 2) cos( 2)

outer gap B R G H    =+ (51)
IV IV IV IV P R R R R P R R G A B R E R R R E R R with P R R R R P R R H C D R E R R R E R R  = - -     =+   (52) where 2 IV A , 2 IV B , 2 IV C , 2
IV D are the second order Fourier coefficients of vector magnetic potential in (66), and R =(R 4 +R 5 )/2. Then, (51) can be further rewritten as: 22 _ ( , ) sin( 2)

outer gap B R G H    = + +
with  =arctan(H/G). So, the fundamental harmonic observed from stator windings is a sinusoidal wave with two pole pairs. If the d-q axis theory in traditional PMSM is applied to the MGDRM, the electrical angle of the equivalent d axis α on the outer air gap produced by PM can be calculated by:

if G>0 2 if G<0 2       -   =   --   (54)
After obtaining the d axis angle for a certain modulator and rotor position, the variation rule among α, modulator and rotor position can be explored. Since the modulator pole structure has a period of 2π/Q and the inner rotor PMs has a period of 2π/P i in circumferential direction, namely Figure 9 is the corresponding contour map of Fig. 8. To illustrate the relation among Ω in Ω out and the angular position of equivalent d axis, different rotating speed combinations are set as following: line 1 represents Ω in = 900 r/min, Ω out = 700 r/min; line 2 represents Ω in = 600 r/min, Ω out = 2000/3 r/min; line 3 represents Ω in = 600 r/min, Ω out = 200/3 r/min; line 4 represents Ω in = -600 r/min, Ω out = -2000/3 r/min. Figure 10 shows the comparison of the variation of the angular position of d axis by analytical method and FEM. For the analytical result, the gradient of the angular position of the d axis is a constant under four different conditions, and the value of that constant changes with different Ω in and Ω out . Actually, it can be generalized that the rotating speed of the d axis Ω d satisfies [START_REF] Bai | Characteristic Analysis and Verification of the Magnetic-Field-Modulated Brushless Double-Rotor Machine[END_REF]:

d i in out PQ  =  - (55)
However, the d axis variation obtained by FEM does not strictly follow these lines, which is caused by magnetic saturation. The overall change rule of Ω d by taking account of magnetic saturation still follows (55).

Hence, to produce a steady torque, the axis rotating magnetic field generated by the stator winding should keep relative still with the equivalent d axis. Hence, the frequency of stator winding should satisfy:

/ 60 d f = (56)
Thus, the speed, torque and power relationship among the three rotating parts of MGDRM are similar to those of planetary gear.

Since the d axis can be calculated by using the angular position information of inner rotor and modulator in real time, i d =0 control can be carried out for the MGDRM. Fig. 11 shows the output torque of inner rotor and modulator variation in one electrical period under i d =0 control with D=5 A/mm 2 , where T out is set as the absolute value for a clear illustration. The torque difference becomes larger with the increase of RMS current density in the stator windings, that is to say, the magnetic saturation of the MGDRM becomes severe as the current in stator windings increases.

d-q Axis Inductance Analysis and id=0 Control

The flux linkage in (49) can be divided into two parts, one is produced by the PMs, named as f  ; another is produced by the stator windings. As for the later one, it can be written in an inductance form, so the flux linkage can be expressed as:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) A fA AA AB AC A B fB BA BB BC B C fC CA CB CC C L M M i M L M i M M L i                                    = +                       (57)
where θ is the angle between the axis of A phase and d axis. L ii is the self inductance of phase i, M ij is the mutual inductance between phase i, j. The self and mutual inductance can be obtained by only letting one phase be electrified and the remanence of PMs be zero. For instance, when the phase A is electrified with direct current I A with an amplitude of 0.1 A, while I B = I C = 0, and B r = 0. The flux linkage matrix ψ can be obtained via (49). Then, the self inductance and mutual inductance can be derived as:

, , I I I C AB AA AB AC A A A L M M   = = = (58)
Figure 13 shows the comparison of L AA , M AB and M AC by analytical model and FEM within an electric period. The saliency effect of MGDRMs is caused by the modulator, whose rotating speed is different from that of the equivalent d axis. So the rule that L AA contains a second harmonic component with respect to electrical angle for conventional PMSMs is false for MGDRMs, as can be seen in Fig. [START_REF] Amrhein | 3-D Magnetic Equivalent Circuit Framework for Modeling Electromechanical Devices[END_REF]. In addition, it can be obtained that the change of L AA is only 0.93 percent due to the saliency effect of the modulator pieces, which means the air-gap permeance observed from stator can be regarded as a constant. Hence, the equivalent synchronous machine is a non-salient-pole one, namely L d = L q , where L d , L q are the d-axis and q-axis inductance, respectively. Hence, under i d =0 control, the output torque of MGDRM can reach its maximum. Besides, L d , L q can be obtained by:

q dq q LL I  == (59)
where ψ q , I q is the q-axis flux and q-axis current, respectively. They can be calculated by using Clark's Transformation and park's Transformation to find the projection of three-phase flux and current on the q axis [START_REF] Jahns | Flux-Weakening Regime Operation of an Interior Permanent-Magnet Synchronous Motor Drive[END_REF], where the q axis is obtained by using the method in section IV part 1.

The values of L d , L q keep constant for a certain MGDRM with the condition on how Ω in and Ω out changing, which simplifies the control of MGDRM after d-q axis decoupling.

Under i d =0 control, the phasor diagram can be plotted as Fig. 14, where the resistance stator winding is ignored. Thus, the power factor of MGDRM can be obtained by:

0 22 0 cos( ) () qq E E L I   = + (60)
where 0 E is the back EMF of one phase obtained from (50), ω is the angular speed of d axis, and ω=2πΩ d /60. The manufacture of double-rotor structure is a tough problem. The air-gap length directly influences the electromagnetic performance of MGDRM. Figure 15 shows the relation of outer rotor torque, q-axis inductance and power factor with the air gap length. It can be seen that the thicker the air-gap length is, the less the outer rotor torque will be. However, a very small air gap requires a high manufacture accuracy. Hence, it is a trade off between the electromagnetic performance and price. Besides, it can be seen from Fig. 15 that the power factor of the MGDRM is very low, that is a common drawback for flux-modulation electric machines due to the high self inductance and low PM flux linkage on the stator windings [START_REF] Oner | Analytical on-Load Subdomain Field Model of Permanent-Magnet Vernier Machines[END_REF].

Demagnetization Capability Analysis

The permanent magnet suffers from the irreversible demagnetization under high-temperature working condition. It is assumed that the ambient temperature for the motors in EV to be 120 °C [START_REF] Pellegrino | Comparison of Induction and Pm Synchronous Motor Drives for Ev Application Including Design Examples[END_REF], and then the demagnetization risk of PMs at different current levels (in RMS value) are analyzed under i d = 0 control. The PM material N35SHX is selected, and its remanence B' r =1.08 T, while at the knee point, B d =0.22 T. Thus, the absolute value of magnetic flux density distribution within PMs should be above B d by considering different magnetization direction, otherwise irreversible demagnetization would happen. Since the outer surface of PMs is the closest to the armature windings, where the PMs have the highest risk of demagnetization, Fig. 16 shows the flux density distribution in the magnetization direction on the outer surface of PMs along the circumferential direction under different armature current. From Fig. 16, the flux density is below 0.22 T when the current density is 8 A/mm 2 . Thus, the maximum safe operation current density of MGDRM is 7 A/mm 2 at 120 °C. Also, a comparison between the analytical prediction and FEM where D= 5 A/mm 2 , temperature is 20 °C and D= 7 A/mm 2 , temperature is 120 °C is shown in Fig. 17-18. It can be observed that the larger the armature current and the higher the temperature are, the smaller of minimum flux density on some places of outer surface of PMs is. Thus, the demagnetization risk goes higher with the increase of armature current and temperature. Besides, a good agreement can be seen between the analytical prediction and FEM.

Conclusion

In this paper, a subdomain model is developed for predicting the magnetic field distribution and electromagnetic performances of MGDRMs. The analytical model can give the accurate results with a significant time saving, when it is compared with the FEM. Moreover, based on the proposed model, the equivalent d-q axis is derived. And the MGDRM is proved to be a non-salient-pole synchronous machine, where the id=0 control can be exerted, and maximal output torque is obtained. Furthermore, the power factor and the demagnetization capability under id=0 control are analysed. Finally, the analytical model is validated by using the FEM. It tells that this model can be a useful tool for the primary stage design and optimization of MGDRMs.

Fig. 1 .

 1 Fig. 1. Geometry of the investigated MGDRM

  Thus, the vector magnetic potential in different regions should follow: 2 0 in region II, III, IV, V A =

  can be calculated by replacing R 4 , R 5 with R 2 and R 3 , and II n

3 Fig. 4 .Fig. 5 .

 345 Fig. 4. Analytically predicted and FEM calculated results with no-load outer air-gap flux density distribution

.

  Fig.8shows the predicted d axis variation with i  and i  .

Fig. 8 . 3 2 1Fig. 9 .Fig. 10 .Fig. 11 .

 83291011 Fig. 8. Variation of d axis with modulator angular position and inner rotor angular position

Fig. 12 .

 12 Fig. 12. Analytically predicted and FEM calculated electromagnetic torque of outer rotor at different current densities

Fig. 13 .Fig. 14 .Fig. 15 .

 131415 Fig. 13. Self inductance and mutual inductance of one phase

Fig. 16 .Fig. 17 .

 1617 Fig. 16. Analytically predicted flux density distribution in the magnetization direction on the surface of PMs at 120 °C at different current levels

Table 1

 1 Boundary Conditions of Each Subdomain

	Region s	Boundary condition	Applied range	Regio ns	Boundary condition	Applied range
	I	0
		1

Table 2

 2 Expressions of Vector Magnetic Potentials in Each Subdomain

	Regions			Corresponding vector potential expression
	I	( , )	0	11 11 2 1 2 1 ( , ) ( , ) ( ( )cos( ))cos( ) ( ( , ) ( , )	( )sin( ))sin( )

Table 3

 3 Geometrical Parameters of MGDRM

	Symbol	Value	Symbol	Value
	R 1	40 mm	P i	7
	R 2	50 mm	P	24
	R 3	52 mm	Q	9
	R 4	62 mm	B r	1.21 T
	R 5	64 mm	r B	1.08 T
	R 6	66 mm	D	5 A/mm 2
	R 7	95 mm	N	100
	L	110 mm	K	100
	β	π/9 rad	T	100
	δ	π/60 rad	M	100
	γ	π/24 rad	-	-
	Figures

Acknowledgement

This work was supported in part by a grant (Project No. 51677159) from the Natural Science Foundation of China (NSFC), China. Also, the work was supported by a grant (Project No. JCYJ20180307123918658) from the Science Technology and Innovation Committee of Shenzhen Municipality, Shenzhen, China. Moreover, it was supported in part by a grant (Project No. CityU 21201216) from the Research Grants Council of HKSAR, China.

Appendix

Table 4 shows the expressions of the Fourier coefficients of each subdomains when boundary conditions are applied.

( , )sin( )

( , )cos( ( ))

( , )cos( ( ))

( , )cos( ( ))


Thierry Lubin
email: thierry.lubin@univ-lorraine.fr

Alexandre Colle

Simulation Analysis and Experimental Evaluation of the Transient Behavior of a Reluctance Magnetic Coupling

The aim of this paper is to present the dynamic behavior of two co-axial rotating shafts connected by a reluctance magnetic coupling. The analysis is based on nonlinear differential equations of motion. The reluctance magnetic coupling is modeled by an equivalent torsional stiffness coefficient and an additional viscous damping coefficient to account for the effects of induced currents in the salient-pole during the transient. These coefficients are obtained from 3D finite element electromagnetic models. In order to evaluate the accuracy of the proposed model, several transient responses for the shaft line are considered. The simulation results are compared with those obtained from measurements.

Introduction

Magnetic couplings are generally used to transmit a rotational movement between two shafts through a sealed wall in order to avoid any contact with the product to be treated. The torque transmission between the two shafts is carried out without contact by the interaction between magnetic fields in the air-gap. There are two main types of magnetic couplings: synchronous and asynchronous (eddycurrent), with radial and axial flux topologies. The synchronous magnetic coupling [1-5] consists of rare-earth permanent magnets (PMs) placed on both rotors that move at the same speed. It has a pull-out torque, which is the maximum torque that can be transmitted by the coupling before the stall. For eddy-current couplings [6-9], one of the two rotors is made with a conductive plate, usually copper. The torque transmission is due to the induced currents which depend on the relative speed between the two rotors (slip). This results in Joule losses, and adequate cooling may be necessary to avoid overheating.

.

. To the best knowledge of the authors, no studies have been conducted to analyze the transient behavior of a reluctance magnetic coupling.

For synchronous magnetic couplings or magnetic gears with permanent magnets placed on both rotors [12],

 [START_REF] Desvaux | Analysis of the dynamic behaviour of magnetic gear with nonlinear modelling for large wind turbines[END_REF][START_REF] Liu | Investigation of the dynamic characteristics of a coaxial magnetic gear under loading condition based on analytical model[END_REF], it has been shown that a stiffness coefficient is sufficient in the transient model to correctly predict the speed response during start-up. Indeed, the impact of induced currents in the magnets during the transient can be neglected due to the low value of their electrical conductivity. An important experimental result given in [START_REF] Lubin | Experimental and theoretical analyses of axial magnetic coupling under steady-state and transient operation[END_REF] shows that the synchronism can be lost for an abrupt start-up of the motor which cause the load to stall. This will be no longer the case for a reluctance magnetic coupling due to the supplementary torque caused by induced currents in the salient-pole during the transient, as shown later in the paper.

For eddy-current magnetic couplings, the dynamic behavior is well predicted by considering only the damping coefficient in the transient model [START_REF] Davis | Transient performance of eddy-current couplings[END_REF][START_REF] Amati | Dynamic behaviour of torsional eddy-current dampers: sensitivity of the design parameters[END_REF][START_REF] Aberoomand | Design optimization of double-sided permanent-magnet axial eddy-current couplers for use in dynamic applications[END_REF]. It is shown that no oscillation appears on the speed response during start-up [START_REF] Lubin | Steady-state and transient performance of axial-field eddy-current coupling[END_REF], as is the case for synchronous magnetic couplings with permanent magnets on both sides [START_REF] Lubin | Experimental and theoretical analyses of axial magnetic coupling under steady-state and transient operation[END_REF]. A comparison of the transient performances between eddy-current coupling and PMs synchronous coupling is given in [START_REF] Lubin | Comparison of transient performances for synchronous and eddy-current torque couplers[END_REF].

For the reluctance magnetic coupling studied in this paper, the analysis is more complicated in terms of modeling because of the induced currents in the salient-pole during the transients that cannot be overlooked to correctly predict the dynamic performances. To obtain good prediction for the dynamic behavior, it will be shown in the next parts that it is necessary to take into account both stiffness and damping coefficients in the transient model.

The paper is organized as follows. The description of the studied reluctance magnetic coupling is presented in section 2. In section 3, the torsional spring stiffness and the viscous damping coefficients representing the reluctance magnetic coupling will be determined by using 3-D finite element simulations. Based on the knowledge of these coefficients, a mathematical model for the dynamic behavior of the shaft line will be developed in section 4 and used to study several transient regimes such as starting. It will be shown in section 5 that the synchronization between the two shafts can be lost in the event of a sudden start if the air-gap is too large. The positive impact of induced currents in the salient-pole for the stability will also be demonstrated.

The simulation results will be compared to experimental tests to demonstrate the accuracy of the proposed mathematical model

Fig. 1. Geometry of the studied reluctance magnetic coupling

Fig. 2. Reluctance magnetic coupling prototype placed on the test bench between two electrical machines

Fig. 3. Equivalent mechanical system of the shaft line with the magnetic coupling

Description of the Axial-Field Reluctance Magnetic Coupling

As shown in Fig. 1, an axial-field reluctance magnetic coupling is composed of two discs facing each other and separated by a small air-gap. One disc is made with permanent magnets (PMs) regularly distributed along the circumferential direction in order to obtain north and south poles. The second disc consists of ferromagnetic salient poles, their number being the same as PMs. If one of the rotors is rotated, there is a change in energy stored in the magnetic circuit. This change is resisted by a restoring torque toward the aligned position. The angular shift between the two discs is directly related to the transmitted torque value. In steady state, the two rotors rotate in synchronism.

The geometrical parameters of the studied reluctance magnetic coupling are given in Fig. 1 and defined in Table 1. The prototype, which is placed between two electric motors, is shown in Fig. 2. The parameters given in Table 1 have been obtained after an optimization procedure but this is not within the scope of this paper. The optimization goals were to obtain a pull-out torque of 15 Nm for an air-gap length of 3 mm for a external radius set at R 2 = 60 mm and the minimum volume for the magnets.

3-D Finite Element Modelling

The first step before studying the dynamic behavior of the shaft line is to determine the mechanical parameters that represent the magnetic coupling as well as possible. The mechanical model of the system is shown in Fig. 3. It is composed of a load connected to a motor (drive torque) via the reluctance coupler and is equivalent to the test bench shown in Fig. 2 where J 1 and B 1 are the moment of inertia and the friction coefficient of the drive rotor, respectively; J 2 and B 2 denote the moment of inertia and the friction coefficient of the load, respectively. T 1 is the drive torque (dc motor), and T 2 is the load torque (dc generator).

As shown in Fig. 3, the reluctance magnetic coupling can be seen as a rotary spring with torsional stiffness K and a viscous rotary damper with a damping factor B connected in parallel. The spring stiffness is directly related to the synchronous torque between the magnets and the salient pole. The damping effect is due to the eddy current loss in the salient-pole that only exists during the transient.

Fig. 4. B-H curve of the salient poles used in 3D FEM

Fig. 5. Reluctance torque as a function of the load angle for 3 values of the air-gap length

Therefore, the torque transmitted through the magnetic coupling can be written as the sum of two terms

    1 2 1 2 c T K B         (1)
where  1 and  2 are the drive and the load rotor position, respectively, and Ω 1 =dθ 1 /dt, Ω 2 =dθ 2 /dt represent the speeds of each of the rotors. In steady state, the two rotors rotate in synchronism ( 1 =  2 ) and the transmitted torque only depends on the load angle defined as

δ =  1 - 2 .
Here, we have assumed that the spring torque is directly proportional to the relative position of its two ends, as well as for the damper torque but for the relative angular speed. We will see in the next section that the reluctance torque coupler has a non-linear behavior and the associate characteristics should be linearised to define the coefficients K and B in [START_REF] Furlani | Analysis and optimization of synchronous couplings[END_REF].

Three dimensional finite element simulations (3-D FEM) are performed using COMSOL Multiphysics® to define K and B. The geometrical parameters given in Table 1 and the actual geometry of the coupling (Fig. 1) are considered for the simulations. The non-linear magnetic property of the iron parts (B-H curve) has been taken into account, Fig. 4.

In order to show the effect of the air-gap length, three values for the air-gap were considered (b = 3 mm, b = 5 mm, and b = 7mm). For the determination of the torsional stiffness coefficient, a 3-D magnetostatic formulation was used while a fully 3-D transient electromagnetic model was required to obtain the viscous damping coefficient. In each case, the mesh size was refined until converging results were obtained [START_REF] Wu | Analysis and design optimisation of magnetic couplings using 3D finite element modelling[END_REF]- [START_REF] Chan | Performance of an Axial-Flux Permanent Magnet Synchronous Generator From 3-D Finite-Element Analysis[END_REF].

Torsional Stiffness Coefficient

The torque as a function of the load angle δ is shown in Fig. 5. The results have been obtained for three values of the air-gap length. As can be seen, the torque varies almost sinusoidally with the load angle δ with a period equal to /p and a peak value T max that corresponds to the pull-out torque. This is a classical result for reluctance-type actuators [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF].

Of course, the pull-out torque decreases as the air-gap length increases. From Fig. 5, the reluctance torque can be approximated by a sinusoidal function as follows

    max 1 2 sin 2 c T T p   (2)
In order to obtain the torsional stiffness expression defined in (1), we have to make a linear approximation of (2) at a given operating load angle δ = δ 0 using the Taylor theorem. The torsional stiffness coefficient can then be written as follows

  max 0 2 cos 2 K pT p  with 0 04 p   (3)
It should be noted that the torsional stiffness depends directly on the pull-out torque and is a multiple of twice the number of pole pairs. Fig. 6 shows the torsional stiffness as a function of the air-gap length. It can be observed that the stiffness decreases when the air-gap increases, and reaches a value of around K = 60 N.m/rad for b = 7mm, which is much lower than the stiffness of mechanical couplings and can be the cause of large torsional oscillation between the two rotors during the transients, as will be seen later in the paper

Viscous Damping Coefficient

The viscous damping coefficient B is much more difficult to determine because 3D FE simulations must be performed by considering transient electromagnetic modeling with motion to predict the induced current in the salient pole due to the relative velocity between the two rotors. The simulations have been conducted using the Rotating Machinery, Magnetic interface in COMSOL Multiphysics.

Fig. 7 shows the eddy current distribution in a salient-pole at a given time and for a relative speed of 10 rad/s between the two rotors. It can be observed that induced currents are mainly distributed on the edge of the teeth due to the skin effect. It needs a fine mesh to converge towards the solution. The simulation times takes about 10 hours to obtain the result given in Fig. 7. From the knowledge of the current density distribution in steady-state condition, the eddy-current losses in salient-poles and back-iron is computed as follows 2 1

( , , )

eddy vol P J x y z dV    ( 4
)
where J is the current density, V is the volume, and  is the iron conductivity given in table 1. Fig. 8 shows the average eddy-current losses as a function of the relative speed between the two rotors for three values of the air-gap. It can be seen that for the considered velocities, the eddy-current losses can be approximated as being proportional to the square of the relative speed and can then be written as follows

  2 12 eddy PB   ( 5
)
where B is the damping coefficient to be determine.

From the eddy-current losses given in Fig. 8, we obtain the damping coefficient as a function of the air-gap length as shown in Fig. 9. As expected, the damping coefficient decreases with the air-gap and reaches a value of around B = 0.01 N.m/(rad/s) for b = 7mm.

Equations of Motion

Nonlinear Model

The governing equations of motion for the mechanical system of Fig. 3 are derived from Newton's second law for rigid body. If we consider a nonlinear spring element for the magnetic coupling as defined in (2), the mathematical model consists of two coupled nonlinear differential equations for the angular shafts position  1 and  2

2 11 1 1 1 2 22 2 2 2 c c dd J B T T dt dt dd J B T T dt dt         (6)
with

    12 max 1 2 sin 2 c dd T T p B dt dt          (7)
The nonlinear model described by equations ( 6) and ( 7) will be simulated in the MATLAB/Simulink environment. The Simulink scheme is shown in Fig. 10.

Linearised Model

If we consider the linearised model for the spring stiffness as defined in (1), we obtain the following coupled linear differential equations for  1 and  2

    2 1 1 2 1 1 1 1 2 2 2 2 1 2 2 2 1 2 d d d J B B K T K B dt dt dt d d d J B B K K B T dt dt dt                     (8)

Fig. 12. Experimental setup

From these linear equations, it is possible to identify the damping factor and the oscillation period of both rotors, as discussed later in the document. In addition, most control design techniques are based on a linear model of the controlled system. Equations ( 8) can therefore be used to define a transfer function between the drive torque and the load rotor position, but this will not be developed here.

Experimental Results

For the experimental validation, we have manufactured an axial-flux reluctance coupler prototype. The prototype, which is placed between two electric motors, is shown in Fig. 2. The geometrical parameters of the prototype are given in Table 1. These geometrical parameters have been obtained after an optimization procedure but this is not within the scope of this paper. The optimization goals were to obtain a pull-out torque of 15 N.m for an air-gap length of 3 mm for a external radius set at R 2 = 60 mm and the minimum volume for the magnets.

The layout of the test bench is shown in Fig. 11 where the reluctance coupler is placed between two dc machines. One of the machines is used as a motor to impose the drive torque T 1 (main motor). It is powered by a dc chopper to control the armature current with a PI controller. A dc generator (secondary motor) is used to apply a sudden load torque by connecting a resistance to its terminals armature winding. Two high-resolution encoders (4096 pulses/revolution) are placed on the test bench to measure the relative angular position between the two rotors and to measure the speed on both sides of the coupling during the transient. The rotor of the dc generator can be locked if necessary. A photograph of the setup is given in Fig. 12.

The parameters of the drive system under consideration are given in Table 2. These parameters have been obtained by classical tests. The armature resistance R a and the armature inductance L a of the dc motor (main motor of Fig. 12) have been identified by the time response of the armature current to a step change in voltage (5V) as shown in Fig. 13. For this test, the rotor was locked.

The friction coefficients B 1 and B 2 were obtained from the no-load test for each dc motor. The motor's moment of inertia J 1 (the same for J 2 ) was measured from the velocity decay rate once current to the dc motor has been cut. For these tests, the two motors were mechanical decoupled, each with a part of the magnetic coupling.

The torsional stiffness and the damping coefficient for the reluctance coupler are given in Fig. 6 and Fig. 9, respectively.

In the following sections, we study two types of transient responses and we compare the simulation results with tests. For the simulation results, we will consider the model with or without damping effect to highlight the differences. The first test concerns the angular rotor position response to a step input in the drive torque T 1 when the other side of the magnetic coupling is locked ( 2 = 0°). As explained previously, the control of the drive torque is done by regulating the armature current of the dc motor.

The experimental results are given in Figs. [START_REF] Lubin | Comparison of transient performances for synchronous and eddy-current torque couplers[END_REF]. As shown in fig. 14(a), the dc motor applies a torque which is zero for t < 0s and increases to a step of 4 N.m at t = 0s. The transient responses for three values of the rotor position are given in Fig. 14(b). We can see large oscillations in the rotor position during the transient. The oscillation period T increases with the air-gap value. This is due to the decrease of the torsional stiffness as shown in Fig. 6.

As the tests were performed for small variation in the rotor position, we can use the linearised model defined in [START_REF] Wang | Analytical Prediction of Torque Characteristics of Eddy Current Couplings Having a Quasi-Halbach Magnet Structure[END_REF] to analyze the results obtained in Figs. [START_REF] Lubin | Comparison of transient performances for synchronous and eddy-current torque couplers[END_REF]. From ( 8) and considering  2 = 0, we obtain the following second-order differential equation between the rotor position and the drive torque.

 

2 11 1 1 1 1 dd J B B K T dt dt       (9)
If we consider an under-damped system that corresponds to the experimental results of Fig. 14(b), the damping ratio ζ and the natural frequency ω n are given by:

1 1 2 BB KJ    1 n K J   (10)
From [START_REF] Nagrial | Analysis and performance of variable reluctance (V. R.) torque coupler[END_REF], we obtain the classical expression for the oscillation period for a second-order system

In order to show the accuracy of the proposed model, Table 3 compares the oscillation period T for simulation and tests. In all cases, the error is less than 10%.

Transient Performance During Start-Up

For the second test, we study the transient performance during start-up by considering two values for the air gap length: 5 mm and 7 mm. The results are given in Fig. 15 and Fig. 16, respectively.

The transient performances are evaluated by running the dc motor from the stop (at t = 0 s) at a speed of 400 rpm under no load condition. A step voltage U dc = 58 V is applied to the dc motor at t = 0 s. Because these tests correspond to large variations in the rotors positions, the simulations were performed using the non-linear model of Fig. 10, adding the dc motor model whose parameters are given in Table 2.

The results obtained for b = 5 mm are given in Fig. 15. It can be observed speed oscillations between the two rotors before the synchronization. This is due to the low stiffness specific to synchronous magnetic couplings [START_REF] Lubin | Experimental and theoretical analyses of axial magnetic coupling under steady-state and transient operation[END_REF]. The simulation result presented in Fig. 15 (dotted line) shows a fairly good agreement with the test. The difference between test and simulation can be explained by the armature inductance variation caused by the magnetic saturation due to the field current. The armature inductance was measured without the field current (Fig. 13). For an airgap of 5 mm, one can observe that the speed of the two rotors remains close to each other during the starting. This is no longer the case if we consider a higher air-gap value, as shown in Fig. 16.

The test result given in Fig. 16(a) corresponds to an air-gap of 7 mm. It shows that for this large air-gap, the start-up is now a two-step process. At the beginning (t < 1.2 s), the speed difference between the two rotors is important.

For this stage, the torque transmitted through the magnetic coupling is mainly due to the damping effect (asynchronous torque). Then, synchronization between the two rotors occurs with large oscillations at t  1.2s. This is due to the synchronous torque between the magnets and the salient poles which becomes predominant. Finally, the two rotors pull into synchronism at around t = 2s. This result is very different from the one we had obtained for a synchronous magnetic coupling with permanent magnets placed on the 2 rotors [START_REF] Lubin | Experimental and theoretical analyses of axial magnetic coupling under steady-state and transient operation[END_REF]. This can be explained from the electrical conductivity of the salient-pole which is around ten times greater than the electrical conductivity of the NdFeB magnets and causes a significant asynchronous torque at start-up.

As shown in Fig. 16(b), the transient behaviour of the shaft-line during start-up is well predicted if the damping effect is taken into account in the simulation. This is no longer the case if we consider the magnetic coupling without damper (B = 0), as can be seen in Fig. 16(c) where the synchronism between the two rotors is lost at the beginning of the start, that does not correspond to the test. This last result clearly shows that it is necessary to represent the reluctance magnetic coupling by the equivalent mechanical model illustrated in Fig. 3, with the two components: spring and damper. It can be noted that the transient behaviour of a reluctance magnetic coupling during the start-up with asynchronous and synchronous torque is similar to that of a line start synchronous reluctance machine with a damper cage winding, although the analytical models for their study are not the same [START_REF] Ogunjuyigbe | Dynamic and transient behavior of a line start, capacitance compensated synchronous reluctance machine[END_REF][START_REF] Mingardi | Line-start PM-assisted synchronous motor design, optimization, and tests[END_REF].

It is worth noting that oscillations between the two rotors during starting (as shown in Fig. 15) can be problematic to achieve high performance for speed or position control in servomechanism applications [START_REF] Montague | Characterisation and modelling of magnetic couplings and gears for servo control systems[END_REF]. If a reluctance magnetic coupling is used to transmit the torque between two shafts, its transfer function with damping and stiffness coefficients (8) must be taken into account when designing the PID controller.

Conclusion

In this paper, we have developed a new tool to effectively predict the transient characteristics of a shaft line including a reluctance magnetic coupling. We have shown that the reluctance coupler causes speed oscillations between the two shafts when starting. One of the most important contributions of this paper is the consideration of the damping effect due to the induced currents in the salient poles during the transient which cannot be neglected in order to obtain good predictions. This damping effect is advantageous because it prevents the stall between the two shafts when starting. Comparisons between experimental and simulation results have clearly demonstrated the value of the proposed model.

Fig. 6 .

 6 Fig. 6. Torsional stiffness as a function of the air-gap length for δ 0 = 0°.

Fig. 7 .

 7 Fig. 7. Eddy current distribution (module) in the salientpole for a relative speed of 10 rad/s between the two rotors.

Fig. 9 .

 9 Fig. 9. Damping coefficient as a function of the air-gap length.

Fig. 8 .

 8 Fig. 8. Eddy-current losses in the salient-pole versus the relative speed.

Fig. 10 .

 10 Fig. 10. Simulink block-scheme for the nonlinear model.

Fig. 11 .Table 2 :Fig. 13 .

 11213 Fig. 11. Test bench arrangement for the transient performance measurement.

Fig. 14 .

 14 Fig. 14. Transient responses to a step input in the drive torque, (a) Drive torque T 1 , (b) Rotor position for b = 3 mm, b = 5 mm, and b = 7 mm, Experimental results.

Table 3 :

 3 Comparison between the model and the test for the oscillation period5.1. Rotor Position Response to a Step Input onthe Drive Torque T 1

Fig. 15 .Fig. 16 .

 1516 Fig. 15. Start-up of the shaft line for an air-gap b = 5mm: ( ──── ) Experimental result, (------) Simulation.

Table 1 :

 1 Geometrical parameters of the magnetic coupling

	Symbol	Quantity	value
	R1	Inner radius of the magnets	30 mm
	R2	Outer radius of the magnets	60 mm
	a	Magnets thickness	10 mm
	b	Air-gap length	variable
	c	Salient-pole thickness	15 mm
	α	PMs pole-arc to pole-pitch ratio	0.83
	β	Slot-opening to pole-pitch ratio	0.77
	p	Pole-pairs number	5
	Br	Remanence of the magnets	1.25 T
		Conductivity of the salient-pole	6 MS/m