Dynamic analysis of Halbach coaxial magnetic gears based on magnetic equivalent circuit modelling

H. Feshki Farahani, Department of Electrical Engineering, Ashtian Branch, Islamic Azad University, Ashtian, Iran. Email: hfeshki@aiau.ac.ir Abstract Halbach permanent magnet arrays (HPMAs) are an attractive configuration that can be applied in PM electric machines such as magnetic gears. These arrays have some merits, namely near‐sinusoidal distribution of flux density along air gap, low torque ripple and iron losses. This article proposes a two‐dimensional magnetic equivalent circuit method to determine the magnetic field distribution of HPMA based on coaxial magnetic gears (HCMG). Furthermore, regarding the proposed model, distribution of radial and tangential components of flux densities is extracted within outer and inner rotors. Moreover, the dynamic model of HCMG is extracted and its pull‐out torque and output speed are calculated using it. Finally, in order to validate the effectiveness of the proposed model, the obtained results for a sample HCMG are compared with finite element analysis method results.


NOMENCLATURE cross section of integral element surface for inner and outer rotor in radial direction
A ragi ; A rago cross section of integral element surface for inner and outer rotor in tangential direction A θagi ; A θago radial and tangential flux density  inner and outer rotors angular positions (degree) θ i , θ o air permeability coefficient μ 0 relative permeability μ r network loops flux vector φ net inner and outer rotor radial air gap flux φ ragi , φ rago inner and outer rotor tangential air gap flux φ θagi , φ θago inner and outer rotor angular velocity ω i ; ω o number of magnet pole pairs for inner and outer rotors p i ; p o

| INTRODUCTION
Gearbox is one of the essential components of industrial applications such as wind power industry, crashers, electro-pumps and electrical vehicles which are required to change the speed and torque. However, mechanical gears play an important role in the aforementioned applications, but they have some disadvantages such as acoustic noises and vibration maintenance problems as well as the need for periodical lubrication [1][2][3].
In order to overcome the mechanical gears deficiencies, a new kind of gear is introduced known as magnetic gears (MGs) [4]. An MG includes two permanent magnet (PM) arrays mounted on the inner and outer rotors. So, it can transform torque and speed via interacting with the magnetic fields of the rotors. This characteristic of the MGs brings out some special advantages including high reliability, low noise, protection against overload, isolation of input and output shafts as well as reducing the requirement of maintenance [2,5,6]. Various kinds of MGs are introduced in the literature such as linear MGs [7], trans-rotary axial MGs [8,9], and coaxial MGs [10]. Each kind has advantages and limitations.
One of the important problems in MG design is maximizing the transmitted torque; so, various literature are published related to MGs structures and their PM arrays. For example some structures like spur type [11,12], worm type [13,14], bevel type [13,15,16] have been proposed, but they have low torque transmission capability because of engaging a few number of magnet poles [15][16][17]. Since, in order to overcome this drawback, the improved topologies are introduced for the MGs, namely the coaxial MGs [4], axial MGs [18] and so on. These structures have some limitations too, such as high torque ripple, non-sinusoidal flux density within the air gaps. One solution to improve their performance is applying Halbach PM Array which can provide higher torque density compared to conventional coaxial MGs [19][20][21][22][23]. As it was mentioned, the Halbach coaxial magnetic gears (HCMG) have some intrinsic capabilities [24] and advantages over conventional coaxial magnetic gears (CMGs), worth considering in MGs design and selection for an application. In the literature, there wasn't found any magnetic equivalent circuit (MEC) modelling study on these new and effective structures. So based on the low computational burden of MEC modelling approach, the present work was carried out to provide a dynamic reluctance network analysis (RNA) model based on MEC theory for HCMGs, usable for fast design purposes.
Various analytical methods are proposed related to MGs modelling such as analytical method [6,[25][26][27][28][29][30], electrical machine theory [31], finite element analysis (FEA) [32][33][34][35] and RNA [36][37][38][39]. The analytical method is based on Maxwell equations [40]. In this method, the Laplace equation is applied for air gap regions and the Poisson equation for the PM sections. In FEA method, the flux density distribution in the different regions of the MG is calculated based on Poisson and Laplace equations [41]. This method is usually used to model the electrical machine, however, it is so time-consuming for the preliminary design stage. In other words, this method can be employed to verify a designed MG and refine it.
The reference [42], adopted a play model to model the magnetic hysteresis behaviour of the motor core in order to calculate the iron loss of the machine. As the authors claimed, this model shortens the calculation time compared with the previous methods, because it needs no convergence calculation. But this model needs a large number of dc hysteresis loops with different maximum magnetic flux densities which can be derived from one or two measurements. So considering the hysteresis behaviour is only reasonable for modelling the machines with high core losses due to the high variation rates of the flux densities, otherwise, it causes an extra computational burden that is not really necessary, especially for the design purposes.
The overall reluctance network-based dynamic analysis in magnetic machines is provided in [43], which finds an equivalent circuit in order to calculate the distributed flux densities and calculate the torque. But in this work because of considering the core nonlinearity, solving the system equations was complicated. So the authors employed the analytical Runge-Kutta's method which needs a powerful simulator like SPICE to find the solution. BEIRAMI ET AL.
In the MGs with only PM magnetic sources (no electric windings), there is no wide range of the core BH operating point, because there is no control over the magneto motive force after the construction. So these devices are designed to work within the linear area of the core BH curve in order to maximize use of the magnets power to produce the maximum torque. But in the Halbach structure the behaviour is some different. In this structure, the flux in the inner rotor and outer rotor cores is weakened and in the pole-piece area, it is strengthened. Also, it must be said the flux in the inner rotor and outer rotor cores due to PMs is generally constant. So modelling hysteresis behaviour for this area is not needed. For the pole pieces, although the flux direction varies with the rotation of the rotors but considering hysteresis loop is not helpful for low-frequency applications. The simulated model in [42] is implemented with a 100-Hz source with high-frequency PWM which really affects the results.
A reluctance network model is employed in [44] for coaxial MG considering very small reluctance elements on the tangential surface of the pole pieces which adds to the accuracy of flux estimation but if the number of circumferential elements are taken enough, this approach won't be very useful but makes code implementation a bit tricky. On the other hand, the authors in the employed reluctance network, in the middle of the pole pieces ring, have just taken one magnetic node which really affects the flux paths reluctances and is not accurate. Also in this model, the cores reluctance in both outer and inner rotors, have been ignored.
The dynamic reluctance network employed in the present work, is very useful for dynamic analysis of Halbach MG, because of the symmetry in the model makes it very easy for code implementation and on the other hand there is no change in the reluctance matrix and for the time simulation, it must be inversed just once at the time zero. Only the magnetic motive forces (MMF) matrix is changing over the time with changing the angular position of the elements and to find fluxes distribution over the network in each time step, there is just need for once matrix multiplication.
The RNA is a much popular modelling method to design of a MG and it helps to save time in preliminary design stage. In this method a MEC is used to calculate flux flow in different points of the MG. Finally, the inner and outer rotor torques could be obtained by some mathematical calculations using flux densities in the air gaps.
To the best of our knowledge, the main contribution of this paper is presenting a model to predict the magnetic components (i.e. radial and tangential flux densities) that are applied to calculate the pull-out torques and speeds within inner and outer rotors. Another novelty of this article is the dynamic analysis of HCMG based on proposed model and extracting torque-angle curve of the inner and outer rotors. Furthermore, the performance of proposed model has been verified with comparing with FEA method results.
The study presented in [45], employed a MEC model for modelling HCMGs similar to the present work, but it has not investigated the performance of this model for dynamic analysis of HCMGs which is the main contribution of the present work. In this article, there is modelled a reluctance network in a way which is very easy for code implementation due to its symmetry and the reluctances are invariant with time or angular position of the rotors as [45]. But in this work, the invariant characteristic of an MG reluctance matrix is used for fast time-domain analysis. So there is just need for one matrix inversion (at time zero) which is the main challenge of MEC model solving in electrical machines. So only the MMF matrix is changing with time or angular position and the inverted matrix is multiplied by the new MMF matrix at instant to find the instantaneous flux distribution.
The article is organized as follows: In Section 2, the principle of HCMG is described. Then extracting of MEC model for these configurations has been presented in Section 3. The obtained results from proposed model and their validation with Maxwell software FEA method results are provided and discussed in Section 4. Finally, conclusions of the present work are discussed in Section 5. Figure 1 shows a conventional CMG and an HCMG respectively the main difference between them are PM array magnetization directions. In other words, in a CMG, the inner and outer rotor PMs are magnetized only in the radial directions, while in an HCMG, those are magnetized in both radial and tangential directions and in some cases in diagonal directions. As shown in Figure 2, the magnetic fields are weakened on one side of the magnetic ring (the yoke side) and it is amplified on the other side (the air gap side) [19]. This property leads to concentrate the fluxes at one side of the magnet ring and increases the flux density. In an HCMG, the magnetic fluxes of the rotors are focused toward pole pieces. This structure leads to utilize the PMs Array effectively. In this configuration, the torque ripples and iron losses are significantly reduced because of the low flux linkage in the iron path. Also, the sinusoidal and amplified flux is generated in the air gap which provides a smoother and greater torque at the output.

| CONVENTIONAL AND HALBACH MAGNETIC GEARS MECHANISM
The number of PM segments per pole (NSP) is usually considered 2 and 3. However selecting higher NSP leads to achieve more sinusoidal air gap flux as well as suppressing the cogging torque, but it increases the fabrication cost. Moreover, the constraints of Halbach cylinders such as the radius and thickness limit the possible value of NSP [19].
In this work, as it is discussed in [19], the number of inner and outer rotors segments per pole are considered 3 and 2 respectively as shown in Figure 2.

| MEC MODELLING OF AN HCMG
The prediction of the magnetic fluxes distribution within the air gaps is a prerequisite to obtain transmitted torque of a MG. MEC analysis models the MG mechanism with lumped circuit parameters. So by a circuit analysis, the magnetic flux flows will be determined in various parts of the gear especially within the air gaps.
The MEC model employed in this work, is a twodimensional model that the HCMG is divided into N segments in tangential direction and 7 layers (6 mesh levels) in radial direction of cylindrical coordinate. This model is shown in Figure 3 in Cartesian coordinate.

| HCMG equivalent net reluctances calculations
Based on Figure 3 and previous descriptions, the magnetic reluctances could be given by (1): where, L stands for the flux flow path length, A stands for the area of cross-section of the path and μ stands for magnetic permeability of the material.
The partial permeance of a unit part in the tangential direction (dP θ ) could be derived from (2) considering Figure 4b: The magnetic reluctance and permeance of the different parts of the equivalent circuit including yokes, pole pieces, magnets and air gaps are calculated by considering Figure 4 and Equations (3) to (5).
In this modelling, it is assumed that the reluctances are approximately linear.

| HCMG equivalent sources modelling
One of the important differences between CMG and HCMG is the tangential components. Due to magnetization direction of HPMAs in the tangential direction, the magneto motive force should be considered in the model and its value depends on in the tangential direction.
Considering Figure 5, the magnetization vectors of each magnet in the Cartesian [20] and Cylindrical coordinate are represented as (6) and (7), respectively: where the sign "-" is for the inner rotor (NSP = 3) and the sign "+" is for the outer rotor (NSP = 2). Also p is the number of pole pairs, θ is the angle between the x-axis and the centre line of each PM and finally M _r is amplitude of remnant magnetization which is calculated as follow where the H c is the field intensity of the PMs and μ 0 is permeability of free space: When a ferromagnetic material is magnetized in one direction, it will not relax back to zero magnetization when the imposed magnetizing field is removed. Every PM material has two kinds of the hysteresis loop, the magnetic induction B curve versus the field intensity H, called extrinsic or normal curve and the magnetization M curve versus H, called intrinsic curve as shown in Figure 6, in electromagnetic machines, where one is concerned with determining the amount of the flux a magnet is capable of producing, the normal demagnetization (B-H) curve is used which illustrates the relationship between the magnetic induction, B, and the applied field strength, H. The intrinsic curve (M-H or J-H) is more useful when working with a PM's reaction to an external magnetic field, which describes the relationship between the magnetization/polarization M/J and H. So in order to determine the PMs Magneto motive forces in a Magnetic circuit, the normal B-H curve is considered. As shown in Figure 6, the operating point of a PM in a magnetic circuit will be in its 2 nd quadrant B-H curve. The magnets used in electrical machines like NdFeB or SmCo usually have a high H ci (intrinsic coercivity) and a linear extrinsic B-H curve in the 2 nd quadrant which makes them ideal for this application. This linear area is affected highly by operating temperature. Increasing the temperature reduces the linear range, bringing the knee-point toward the 2 nd quadrant (usually in temperatures above 120°).
In cases where no demagnetization field is applied, the output magneto motive force of the PM can be written as: The intersection of the load B-H line (based on the circuit shown in Figure 7  Ignoring the leakage fluxes and fringing effects, the circuit flux can be written as: And by setting the MMF of the PM and equivalent circuit equal: By dividing the two sides of Equations (10) and (11), and considering a constant B m /H m ratio (the circuit elements are air and linear core) as a permeability μ d , the load line equation can be obtained as: Obtaining the operating points and the load lines for the complicated structures containing several PMs and multi flux passing loops, is required for solving several simultaneous equations which is presented in this article.
As it is shown in Figure 8, to simplify the equations and making the reluctances matrix symmetrical, the reluctances and sources are divided into two equal sections.
The MMF of the inner and outer rotors magnets are modelled as voltage sources as it is shown in Figure 8. The MMF of the PM section can be written as: where H ck is magnetic coercivity of the magnets, l mk is the magnetization length of the magnet elements, the subscript k stands for inner and outer rotors (k ∈ fi; og) and μ r is relative permeability of the magnets. So, for a unit magnet element as Figure 8a, the radial and tangential sources are calculated as follow: where l k and w k are average tangential and radial length of a PM segment respectively (Figure 8b). The MMF of each PM element is divided into two equal parts in the radial direction to increase the precision of the model. Also, the air gaps, magnets and pole pieces reluctances are divided into two equal parts only in the radial direction.

| HCMG MEC solution
Based on the RNA method, the HCMG could be divided into N sections as shown in Figure 3 that N is an integral number. In other words, the angle of each section would be 360/N. By dividing the CMG into seven radial layers, including (1) outer rotor yoke, (2) outer rotor magnets, (3) outer air gap, (4) pole pieces stationary ring, (5) inner air gap, (6) inner rotor magnets and (7) inner rotor yoke layer ( Figure 3). There will be 6�N circuit meshes that the simultaneous mesh equations (based on Kirchhoff's voltage law) need to develop and solve. According to Figure 3, to find the mesh fluxes, the 6�N equations are provided as follow: where, R Net and F Net are network reluctances and sources matrixes. The fluxes matrix can be presented as: That the flux distribution in each part of the circuit can be given by solving (18): After achieving the flux matrix, tangential and radial components of the flux within the inner (φ ragij and φ θagij ) and outer air gaps (φ ragoj and φ θagoj ) can be calculated.
The tangential and radial flux density components at air gaps can be specified as follows:

| HCMGS TORQUE AND SPEED CALCULATION
The torques of the gear can be determined by the following function based on magnetic energy [19]: where B θ and B r are the components of the air gap flux density and r is the radius of the air gap. The integration of (20) is implemented for whole tangential elements and the torque is calculated for each air gap. Equation (20) is called stress tensor equation which is applied for rotors surface elements [19]. So, the surface flux densities components are considered as shown in Figure 9.  Discretizing the integration (20), the output and inner rotors torques can be calculated as The non-zero average torque of the outer rotor is generated due to the interaction of inner and outer rotors fluxes. The rotor dynamic (force-movement) equation is expressed as where T Load and J o are load torque and moment inertia of the outer rotor, respectively. The speed of the outer rotor (ω o ) and its position (θ o ) at each time interval (dt) can be calculated by (25) and (26): The flowchart of HCMG dynamic model is depicted in Figure 10. In this model, at each instant, the torques are extracted based on HCMG geometric parameters such as

| SIMULATION RESULTS
To validate the model effectiveness, this model has been implemented using MATLAB script and the results have been compared with finite element method (FEM) results. The parameters of HCMG are listed in Table 1. The numbers of pole pairs of the inner and outer rotors are 4 and 17, respectively which leads to gain 1:4.25 gear ratio. The number of pole pieces is equal to sum of the inner and outer pole pairs (17 + 4 = 21). The material used in the MG is Nd-Fe-B that its residual flux density is 1.25 T and coercive magnetic force is 975 kA/m. The core material of the pole pieces and the rotors yoke is non-oriented silicon steel.
In order to validate the proposed analytical model for HCMG, some simulation scenarios have been considered with the proposed model, compared with those obtained with twodimensional (2D) finite element simulation using Ansoft Maxwell transient modelling. The main advantages of FEA are its high accuracy and reliability. This method has been employed to verify the analytical MEC approach. Both models are performed with the same assumptions as shown in Table 1.
In Figure 11, the waveforms of radial and tangential components of the flux densities in inner and outer air gaps are shown. The results which are extracted by taking the number of MG tangential divisions (N ) equal to 450 in MEC model show that there is a good agreement between the analytical and FEA results. Therefore, the proposed 2D equivalent magnetic circuit model of Halbach MGs has the capability of predicting the distribution of the flux densities within the air gaps of this type of MGs. Figure 12 shows flux density distribution and magnetic flux lines of the simulated HCMG, respectively. It is obviously observed that in Halbach MG, the flux density in the inner and outer rotor yokes are much lower than the main torque producing parts and there are low flux lines.
By holding the outer rotor in the stationary condition and rotating the inner rotor for one period of electrical degree, the torque capability (static torque curve) of the gear will be achieved that is depicted in Figure 13.
It can be seen that the obtained torque curve exactly tracks the one obtained from FEM simulation. The maximum torque calculated using MEC model on the outer and inner rotors is 174.4 and 40.9 N.m, respectively. So, the obtained ratio between two torques is 4.26 (174.4/40.9) that has only 0.3 percent tolerance compared with the expected ratio 4.25 (17/4).
From the MGs theory, when the pole pieces number is selected equal to sum of the rotors pole pair's numbers, the rotors will rotate in opposite direction, when the pole pieces ring is held stationary [46].  The torque and speed behaviour of HCMG is obtained from MEC-based dynamic analysis which is plotted in Figure 14. This behaviour is extracted when the load torque is gradually increased in three steps while the inner rotor rotates with constant speed (1000 rpm). The results indicate that the proposed model properly calculates the transmitted torque and speed with the ratio 1:4.25 before the load torque exceeds the gear maximum achievable torque. It is observed that the proposed model considers the gear's loss of synchronism when the gear is loaded heavier than its maximum allowable torque.
In order to compare dynamic model results, the HCMG with the same parameters, is modelled and simulated in Maxwell soft environment in system with a two-core processor and 4G RAM. The rotational speed and torque behaviour transferred to the outer rotor is shown in Figure 15. The transient simulation of 7.5 s performance of the HCMG using FEM method took about 21 h comparing with MEC method which was about 30 min. Comparing the results shows that only in fast speed and load switching and also during loss of synchronism, the results differ a bit. Also, it can be seen there is a bit ripple type oscillation in the transferred torque which is related to the MEC circuit, tangential elements number and with selecting a larger number, results will be more similar. But this action adds to the simulation time which is still faster than the FEM method. The selected reluctance matrix shows a good performance for fast time-domain analysis ( Table 2), because there is just a need for once matrix inversion (at time zero) which is the main challenge of solving the MEC model in electrical machines. So only the MMF matrix is changing with time or angular position and the inverted matrix is multiplied by the new MMF matrix at instant to find the instantaneous flux distribution.

| CONCLUSION
This article proposed the modelling of coaxial MG with Halbach PM array based on 2D equivalent magnetic circuit (MEC) modelling approach. From the obtained results, distribution of magnetic fields within the MG and the pullout torque have been evaluated and based on FEA validation, it is illustrated that the proposed analytical model can estimate the system operation with high accuracy. Compared with the FEA method, the proposed method takes less computational effort in estimating the radial and tangential components in different points of HCMGs presenting a dynamic behaviour. As a result, this approach will be a useful tool for HCMGs design especially for dynamic performance-oriented designs.

Parameter Value
Mesh resolution 1 mm