Analytical Modeling of Static Eccentricities in Axial Flux Permanent-Magnet Machines with Concentrated Windings

Yunkai Huang 1, Baocheng Guo 1,*, Ahmed Hemeida 2 and Peter Sergeant 2 1 Engineering Research Center for Motion Control of Ministry of Education, Southeast University, Nanjing 210096, China; huangyk@seu.edu.cn 2 Department of Electrical Energy, Systems and Automation, Ghent University, Ghent B-9000, Belgium; a.hemeida@live.com (A.H.); Peter.Sergeant@ugent.be (P.S.) * Correspondence: guobaocheng1986@gmail.com; Tel.: +86-150-0514-5661


Introduction
Axial flux permanent magnet machines (AFPMMs) have a number of distinct advantages over radial flux permanent magnet machines (RFPMMs).Axial flux permanent magnet (AFPM) can be designed to have a higher power-to-weight ratio, resulting in a need for less core material and reduced complexity to adjust the air gap [1].Moreover, because of the disc shaped rotor and stator structure, an AFPM is smaller in size than its RFPMMs counterparts.This particular structural characteristic makes it easy to address the space limitation in some applications such as electric vehicles (EV) [2], hybrid electric vehicles (HEV) [3], wind turbos [4] and flywheel storage systems [5], etc.
However, due to the smaller contact surface between the rotor and the shaft, it is more difficult to design a rotor-shaft mechanical joint with a high mechanical integrity [6], which makes AFPMMs more susceptible to imperfect assembly issues such as static eccentricities (SEs), dynamic angular misalignment [7] and static/dynamic eccentricity misalignment [8].Under angular eccentric conditions, the rotor is inclined and the air gap is asymmetric.If the rotor shaft assembly is sufficiently rigid, which means that the level of static eccentricity (SE) does not change during periodic operation, the axial force becomes unbalanced causing vibrations and unbalanced magnetic forces (UMFs).This Figure 1a shows the construction of the test machine used in this paper.The AFPM machine consists of two rotors and one segmented stator.Soft magnetic composite (SMC) is used for the stator core and axially magnetized NdFeB permanent magnets (PMs) are mounted on the surface of each rotor, which is made of #45 steel, the SMC modules are fastened by resin and holders.Three phase concentrated windings are distributed on each stator segment to avoid overlapping.This machine topology benefits from advantages including short end windings with high winding fill factor, small stator iron loss under high-speed operation and a convenient manufacturing and assembly process, plus reduced stator weight due to the absence of the stator yoke.The dimensions and specifications of the investigated machine are shown in Table 1.

Calculation of the Magnetic Field
In this paper, the Quasi-3D-method is adopted in the analytical model.The machine is divided into a certain number of layers.Thus, an axial machine could be considered as a combination of several individual linear machines, which allows separate analysis of each plane.One 2D model developed in 2D polar coordinates is shown in Figure 1d and the average radius of a particular ith layer is given by: where ns is the number of slices.
To improve the accuracy of the results, six slices are considered to calculate the healthy and fault conditions in this paper.The processes of quasi 3D method and two poles of a 2D model in Cartesian coordinates are illustrated in Figure 1b-d.
As for the magnetic flux calculation of air gap, the same approach in [20] is used to calculate the magnetic induction in the air gap caused by the magnets and the armature current.Moreover, to take the slot effect into account, the SC mapping [21] is adopted.
By applying Maxwell's equations on selected slice and by summing the air gap flux caused by the armature current and magnets, the flux distribution can be obtained.The following assumptions are made during calculations in order to reduce the complexity of the computations: (1) The magnetic material has a uniform magnetization and the relative recoil permeability μr is constant and has a value close to unity such as in NdFeB materials.(2) For the computation of armature reaction field, the magnet regions are regarded as free space.

Calculation of the Magnetic Field
In this paper, the Quasi-3D-method is adopted in the analytical model.The machine is divided into a certain number of layers.Thus, an axial machine could be considered as a combination of several individual linear machines, which allows separate analysis of each plane.One 2D model developed in 2D polar coordinates is shown in Figure 1d and the average radius of a particular ith layer is given by: where n s is the number of slices.
To improve the accuracy of the results, six slices are considered to calculate the healthy and fault conditions in this paper.The processes of quasi 3D method and two poles of a 2D model in Cartesian coordinates are illustrated in Figure 1b-d.
As for the magnetic flux calculation of air gap, the same approach in [20] is used to calculate the magnetic induction in the air gap caused by the magnets and the armature current.Moreover, to take the slot effect into account, the SC mapping [21] is adopted.
By applying Maxwell's equations on selected slice and by summing the air gap flux caused by the armature current and magnets, the flux distribution can be obtained.The following assumptions are made during calculations in order to reduce the complexity of the computations: (1) The magnetic material has a uniform magnetization and the relative recoil permeability µ r is constant and has a value close to unity such as in NdFeB materials.(2) For the computation of armature reaction field, the magnet regions are regarded as free space.
Energies 2016, 9, 892 4 of 19 (3) Magnetic saturation is absent and the rotor iron cores have infinite magnetic permeability.(4) Eddy current effects are neglected, which avoids the need for the complex eddy current field formulation.

Model of Permanent Magnets (PMs)
For NdFeB, the magnetic induction in the PMs is written as: where µ rec = µ 0 µ r is the recoil permeability.As shown in Figure 2, the magnetization vector is assumed to be along the axial direction and may be described by Fourier series containing only cosine terms: where τ p is pole pitch in circumferential direction and α p is the ratio of magnet pole arc to pole pitch, h m is axial thickness of magnets, L is the axial distance between rotor and back plates and M n is shown as follows: For a PM machine with linear demagnetization characteristic, the scalar magnetic potentials (ϕ) in both the air space and the PMs are governed by Laplace's equation when the rectangular coordinate system is adopted [21].
The components of magnetic strength are related to ϕ by: The following boundary conditions are applied to solution of Equation ( 5): By applying the previous boundary conditions, the axial and circumferential flux density (B y and B x respectively) in region I and II are solved.For region I (the air space): where:

Model of Armature Reaction Current
Taking into account the armature winding current effect, the fractional slot windings are regarded as thin wires.The current density J can be obtained by multiplying the value of current sheet distribution of each phase shown in Figure 3, and dividing by the thickness of slot opening tso as follows: The current sheet spatial distribution of slots in the stator in per unit (JA, JB and JC) of each phase is plotted in Figure 3 [21].Afterwards, the current density J is expressed by Fourier series to obtain the spatial distribution of the current, which could be derived as: where An and Bn are the coefficients obtained from Fourier series.As shown in Figure 4, the scalar magnetic potentials ϕ due to the current also can be described by Laplace's equation.To solve the equation, the following boundary conditions are applied: where L′ is the axial position of a typical armature winding current sheet as shown in Figure 4.

Model of Armature Reaction Current
Taking into account the armature winding current effect, the fractional slot windings are regarded as thin wires.The current density J can be obtained by multiplying the value of current sheet distribution of each phase shown in Figure 3, and dividing by the thickness of slot opening t so as follows: ))/t so (11) The current sheet spatial distribution of slots in the stator in per unit (J A , J B and J C ) of each phase is plotted in Figure 3 [21].

Model of Armature Reaction Current
Taking into account the armature winding current effect, the fractional slot windings are regarded as thin wires.The current density J can be obtained by multiplying the value of current sheet distribution of each phase shown in Figure 3, and dividing by the thickness of slot opening tso as follows: The current sheet spatial distribution of slots in the stator in per unit (JA, JB and JC) of each phase is plotted in Figure 3 [21].Afterwards, the current density J is expressed by Fourier series to obtain the spatial distribution of the current, which could be derived as: where An and Bn are the coefficients obtained from Fourier series.As shown in Figure 4, the scalar magnetic potentials ϕ due to the current also can be described by Laplace's equation.To solve the equation, the following boundary conditions are applied: where L′ is the axial position of a typical armature winding current sheet as shown in Figure 4. Afterwards, the current density J is expressed by Fourier series to obtain the spatial distribution of the current, which could be derived as: where A n and B n are the coefficients obtained from Fourier series.As shown in Figure 4, the scalar magnetic potentials ϕ due to the current also can be described by Laplace's equation.To solve the equation, the following boundary conditions are applied: (13) where L is the axial position of a typical armature winding current sheet as shown in Figure 4.The PMs are physically located in the lower region III and the armature current sheet is very close to the stator back plate.The result in region III will be: 1,3,5,... where:

Model of Stator Slotting
The effect of stator slotting can be linked to the previous flux density by defining a vector potential (λ) in each slot, the axial and circumferential components of flux density can be written in the form [18]: To this aim, the CM method is used to calculate the λ by considering the slot effect in electrical machines, but this method has a defect mentioned above.The numerical SC mapping and its Matlab SC Toolbox could draw the real slot shape [22]: where λ0 is the slotless air gap complex permeance in the T-plane, λx and λy are the circumferential and axial components of the complex relative air gap permeance in the original Z-plane.Furthermore, K, T, Z and S represent the K-plane, T-plane, Z-plane and S-plane, respectively.It should be noticed that logarithmic mapping is used for the axial flux electric machine to convert the circular geometry in the S-plane to the Z-plane.However, the AFPM machines in the quasi 3D model (Figure 1) can be regarded as linear permanent magnet machines presented by the Z-plane.Hence, the SC mapping only transforms Z-plane to T-plane.However, in order to calculate the relative air gap permeance via Hague's solution, the complex K plane is introduced in this paper.Thus, the vector potential can be deduced as: The SC transformation in this step maps one canonical domain (e.g., a rectangle, disk, bi-infinite strip, or the upper (lower) half-plane) into the interior (exterior) of the corresponding polygon.This transformation is defined as follows: The PMs are physically located in the lower region III and the armature current sheet is very close to the stator back plate.The result in region III will be: where:

Model of Stator Slotting
The effect of stator slotting can be linked to the previous flux density by defining a vector potential (λ) in each slot, the axial and circumferential components of flux density can be written in the form [18]: To this aim, the CM method is used to calculate the λ by considering the slot effect in electrical machines, but this method has a defect mentioned above.The numerical SC mapping and its Matlab SC Toolbox could draw the real slot shape [22]: where λ 0 is the slotless air gap complex permeance in the T-plane, λ x and λ y are the circumferential and axial components of the complex relative air gap permeance in the original Z-plane.Furthermore, K, T, Z and S represent the K-plane, T-plane, Z-plane and S-plane, respectively.It should be noticed that logarithmic mapping is used for the axial flux electric machine to convert the circular geometry in the S-plane to the Z-plane.However, the AFPM machines in the quasi 3D model (Figure 1) can be regarded as linear permanent magnet machines presented by the Z-plane.Hence, the SC mapping only transforms Z-plane to T-plane.However, in order to calculate the relative air gap permeance via Hague's solution, the complex K plane is introduced in this paper.Thus, the vector potential can be deduced as: The SC transformation in this step maps one canonical domain (e.g., a rectangle, disk, bi-infinite strip, or the upper (lower) half-plane) into the interior (exterior) of the corresponding polygon.This transformation is defined as follows: where A and C are unknown complex constants, n is the number of polygon corners with interior angle α t , t 1 , . . ., t k are the points in the canonical domain (in the T-plane) corresponding to the polygon corners.
The analytical solution of Equation ( 20) is very difficult for geometries with more than three vertices [22], therefore, the SC Toolbox provides a numerical solution by a library of command-line functions which could be seen in [19].In this paper, the rectangle domain (T-plane) is created by calling f = correctmap (p, alpha), p is the point and alpha is the corresponding angle as shown in Figure 5.
Energies 2016, 9, 892 7 of 18 where A and C are unknown complex constants, n is the number of polygon corners with interior angle αt, t1, …, tk are the points in the canonical domain (in the T-plane) corresponding to the polygon corners.
The analytical solution of Equation ( 20) is very difficult for geometries with more than three vertices [22], therefore, the SC Toolbox provides a numerical solution by a library of command-line functions which could be seen in [19].In this paper, the rectangle domain (T-plane) is created by calling f = correctmap (p, alpha), p is the point and alpha is the corresponding angle as shown in Figure 5.The next step aims to map the interior of annular domain in K-plane to rectangular domain in T-plane.This mapping can be calculated analytically: With the determinate rectangular domain vertices by calling vt = evalinv (f, vz), the length (Δx) and width (Δy) of the T-plane could be calculated.
Finally, Hague solution for magnetic field in an annular domain could be used, and the final equation which maps the field solution from the K-plane to Z-plane is given as Equation ( 19), where / ∂ ∂ T Z can be obtained in SC Toolbox by calling evaldiff (f, Z) and: Afterwards, the magnetic flux density at slotted air gap field could be calculated by Equation (17).The next step aims to map the interior of annular domain in K-plane to rectangular domain in T-plane.This mapping can be calculated analytically: With the determinate rectangular domain vertices by calling v t = evalinv (f, v z ), the length (∆x) and width (∆y) of the T-plane could be calculated.
Finally, Hague solution for magnetic field in an annular domain could be used, and the final equation which maps the field solution from the K-plane to Z-plane is given as Equation (19), where ∂T/∂Z can be obtained in SC Toolbox by calling evaldiff (f, Z) and: Afterwards, the magnetic flux density at slotted air gap field could be calculated by Equation (17).

Model of Static Eccentricity
SE is a misalignment condition of stator and rotor axis, thus, the rotor axis is deflected from that of the stator which causes a non-uniform air gap.In a healthy state (Figure 6a), the air gap length g 0 is uniform along the circumferential direction.However, with the occurrence of SE (Figure 6b,c), the rotor shaft experiences a deflection and the symmetry of the rotor deviates from that of the stator by an angle β.The air gap varies from small to large around the circumferential of the stator but the relative position does not vary with time.In other words, in the case of SE the position of the minimum air gap length is fixed in space.

Model of Static Eccentricity
SE is a misalignment condition of stator and rotor axis, thus, the rotor axis is deflected from that of the stator which causes a non-uniform air gap.In a healthy state (Figure 6a), the air gap length g0 is uniform along the circumferential direction.However, with the occurrence of SE (Figure 6b,c), the rotor shaft experiences a deflection and the symmetry of the rotor deviates from that of the stator by an angle β.The air gap varies from small to large around the circumferential of the stator but the relative position does not vary with time.In other words, in the case of SE the position of the minimum air gap length is fixed in space.The first step in modeling the eccentricity is to evaluate the air gap length in SE condition [23], the SE factor can be defined as: The air-gap length at mean radius Rmid can be written as follows: where θ is the stator position measured from a reference point (γ0) of the minimal air gap.According to Figure 6b, the air gap length in all positions can be understood as: where g' is the air gap deviation at the Rmid: Therefore, the air gap length in all positions can be written as follows: Because of the high ratio of machine diameter to length in the AFPMs, 40% SEF means that the maximum declined distance of rotor plate is 1.2 mm in this case, which is already a significant value compared to the air gap length of 3 mm.
In this paper, as quasi-2D model processed in Figure 1, the minimal air gap happens at γ0 = 0 rad and the maximum air gap at θ= π rad.Finally, the air gap length can be deduced as follows:

Bilinear Mapping
Bilinear mapping [24] is a particular mapping method which is used for the treatment of the boundaries which are circular but eccentric or intersecting.The purpose of this process is to map the The first step in modeling the eccentricity is to evaluate the air gap length in SE condition [23], the SE factor can be defined as: The air-gap length at mean radius R mid can be written as follows: where θ is the stator position measured from a reference point (γ 0 ) of the minimal air gap.According to Figure 6b, the air gap length in all positions can be understood as: where g' is the air gap deviation at the R mid : Therefore, the air gap length in all positions can be written as follows: Because of the high ratio of machine diameter to length in the AFPMs, 40% SEF means that the maximum declined distance of rotor plate is 1.2 mm in this case, which is already a significant value compared to the air gap length of 3 mm.
In this paper, as quasi-2D model processed in Figure 1, the minimal air gap happens at γ 0 = 0 rad and the maximum air gap at θ= π rad.Finally, the air gap length can be deduced as follows:

Bilinear Mapping
Bilinear mapping [24] is a particular mapping method which is used for the treatment of the boundaries which are circular but eccentric or intersecting.The purpose of this process is to map the pair of eccentric circles (E-plane) into a pair of concentric circles (S-plane).After mapping E-plane to the S-plane, the following mapping could be done by previous study.The general form of this CM can be defined as: where parameters C and K are determined by the rotor and stator radii and the magnitude of eccentricity.
Markovic used this method to calculate a non-salient motor under SE.Bilinear maps preserve the scalar magnetic potential but they cannot conserve the vector field potential (Ω v ) because flux density in concentric plane is coupled with the coordinate system.Since vector field potential calculation is complex and it is very difficult to estimate the scale factor required to map magnetic field into desired eccentric domain, in this paper, the following form is used to convert the air gap flux density of concentric AFPM motor in Z-plane to the eccentric model in E-plane [25]: where λ e is the relative permeance function due to the effect of SE and can be defined as: The flux density in the slotless air gap with eccentricity (B ecc ) and without eccentricity (B no-ecc ) can be estimated by modifying the air gap length in Equation ( 9).
The calculation processes can be seen in Figure 7.
pair of eccentric circles (E-plane) into a pair of concentric circles (S-plane).After mapping E-plane to the S-plane, the following mapping could be done by previous study.The general form of this CM can be defined as: where parameters C and K are determined by the rotor and stator radii and the magnitude of eccentricity.
Markovic used this method to calculate a non-salient motor under SE.Bilinear maps preserve the scalar magnetic potential but they cannot conserve the vector field potential (Ωv) because flux density in concentric plane is coupled with the coordinate system.Since vector field potential calculation is complex and it is very difficult to estimate the scale factor required to map magnetic field into desired eccentric domain, in this paper, the following form is used to convert the air gap flux density of concentric AFPM motor in Z-plane to the eccentric model in E-plane [25]: where λe is the relative permeance function due to the effect of SE and can be defined as: The flux density in the slotless air gap with eccentricity (Becc) and without eccentricity (Bno-ecc) can be estimated by modifying the air gap length in Equation ( 9).
The calculation processes can be seen in Figure 7.

Start Set initial value
Calculate the slotless magnetic flux density under healthy condition via Eq.8 and Eq.9 Air

Simulation Verification
To validate the analytical model, an AFPMM with nine slots and six poles has been studied.The basic parameters are shown in Table 1.The air gap flux distribution is one of the most important characteristics in electrical machines because any change of this feature would affect the performance of motor.In this section, as the first step of investigation, the air gap flux density under healthy condition is presented.The results of the analytical model are compared with the results calculated by FE model.As the second step, the flux density under SE at the same conditions is also compared

Simulation Verification
To validate the analytical model, an AFPMM with nine slots and six poles has been studied.The basic parameters are shown in Table 1.The air gap flux distribution is one of the most important characteristics in electrical machines because any change of this feature would affect the performance of motor.In this section, as the first step of investigation, the air gap flux density under healthy condition is presented.The results of the analytical model are compared with the results calculated by FE model.As the second step, the flux density under SE at the same conditions is also compared with the FE results.Moreover, the effect of eccentricity on air gap flux of AFPM will be discussed in the final step.

Air Gap Flux under Healthy Condition
The axial and circumferential components of complex relative permeance in the middle of air gap are shown in Figure 8.It shows that, under healthy condition, the complex relative air gap permeance is periodic for each slot pitch.
Energies 2016, 9, 892 10 of 18 with the FE results.Moreover, the effect of eccentricity on air gap flux of AFPM will be discussed in the final step.

Air Gap Flux under Healthy Condition
The axial and circumferential components of complex relative permeance in the middle of air gap are shown in Figure 8.It shows that, under healthy condition, the complex relative air gap permeance is periodic for each slot pitch.With Equation ( 17), the axial and circumferential components of the flux density in the slotted air gap are given by: The axial and circumferential flux density at the center of air gap under no-load condition are shown in Figure 9.With Equation ( 17), the axial and circumferential components of the flux density in the slotted air gap are given by: The axial and circumferential flux density at the center of air gap under no-load condition are shown in Figure 9. Figure 10 shows the flux density at the rated current (10 A), 20 krpm.The figures show that the analytical model is in a very good agreement with the FE model.with the FE results.Moreover, the effect of eccentricity on air gap flux of AFPM will be discussed in the final step.

Air Gap Flux under Healthy Condition
The axial and circumferential components of complex relative permeance in the middle of air gap are shown in Figure 8.It shows that, under healthy condition, the complex relative air gap permeance is periodic for each slot pitch.With Equation ( 17), the axial and circumferential components of the flux density in the slotted air gap are given by: The axial and circumferential flux density at the center of air gap under no-load condition are shown in Figure 9.  with the FE results.Moreover, the effect of eccentricity on air gap flux of AFPM will be discussed in the final step.

Air Gap Flux under Healthy Condition
The axial and circumferential components of complex relative permeance in the middle of air gap are shown in Figure 8.It shows that, under healthy condition, the complex relative air gap permeance is periodic for each slot pitch.With Equation ( 17), the axial and circumferential components of the flux density in the slotted air gap are given by: The axial and circumferential flux density at the center of air gap under no-load condition are shown in Figure 9.

Air Gap Flux Density under SE Condition
In order to verify the SE model, the rotor of AFPM is assumed with 40% eccentricity.It is shown in previous section that the slotless flux density caused by PMs under SE and healthy condition could be solved analytically, thus, the relative permeance function of SE can be determined.The FE model is represented by two rotor disc to ensure the real flux distribution.The 3D FE model and its grid could be seen in Figure 11.The rotor axis is set as shown in Figure 6a.

Air Gap Flux Density under SE Condition
In order to verify the SE model, the rotor of AFPM is assumed with 40% eccentricity.It is shown in previous section that the slotless flux density caused by PMs under SE and healthy condition could be solved analytically, thus, the relative permeance function of SE can be determined.The FE model is represented by two rotor disc to ensure the real flux distribution.The 3D FE model and its grid could be seen in Figure 11.The rotor axis is set as shown in Figure 6a.The axial and circumferential components of one side rotor disk with 40% eccentricity under noload and rated load condition conditions are shown in Figures 12 and 13, respectively.It can be seen that the axial and circumferential components of the flux density between the FE model and analytical method do not agree well due to the interaction between the circumferential and axial flux density, but the error still remains in a reasonable range.Overall, the proposed method in this paper is an effective way to investigate the eccentric influence of AFPM.The axial and circumferential components of one side rotor disk with 40% eccentricity under no-load and rated load condition conditions are shown in Figures 12 and 13, respectively.It can be seen that the axial and circumferential components of the flux density between the FE model and analytical method do not agree well due to the interaction between the circumferential and axial flux density, but the error still remains in a reasonable range.Overall, the proposed method in this paper is an effective way to investigate the eccentric influence of AFPM.

Air Gap Flux Density under SE Condition
In order to verify the SE model, the rotor of AFPM is assumed with 40% eccentricity.It is shown in previous section that the slotless flux density caused by PMs under SE and healthy condition could be solved analytically, thus, the relative permeance function of SE can be determined.The FE model is represented by two rotor disc to ensure the real flux distribution.The 3D FE model and its grid could be seen in Figure 11.The rotor axis is set as shown in Figure 6a.The axial and circumferential components of one side rotor disk with 40% eccentricity under noload and rated load condition conditions are shown in Figures 12 and 13, respectively.It can be seen that the axial and circumferential components of the flux density between the FE model and analytical method do not agree well due to the interaction between the circumferential and axial flux density, but the error still remains in a reasonable range.Overall, the proposed method in this paper is an effective way to investigate the eccentric influence of AFPM.

Air Gap Flux Density under SE Condition
In order to verify the SE model, the rotor of AFPM is assumed with 40% eccentricity.It is shown in previous section that the slotless flux density caused by PMs under SE and healthy condition could be solved analytically, thus, the relative permeance function of SE can be determined.The FE model is represented by two rotor disc to ensure the real flux distribution.The 3D FE model and its grid could be seen in Figure 11.The rotor axis is set as shown in Figure 6a.The axial and circumferential components of one side rotor disk with 40% eccentricity under noload and rated load condition conditions are shown in Figures 12 and 13, respectively.It can be seen that the axial and circumferential components of the flux density between the FE model and analytical method do not agree well due to the interaction between the circumferential and axial flux density, but the error still remains in a reasonable range.Overall, the proposed method in this paper is an effective way to investigate the eccentric influence of AFPM.

Results and Discussion
In this section, some important machine characteristics, e.g., back-EMF and torque, will be presented and compared.Since these characteristics under different SE conditions only differ in magnitude, the 40% eccentricity is chosen to carry out the analysis.
As shown in Figure 14, the flux passes from one rotor disc and enters the other rotor disc.Therefore, the coils are assumed to be symmetrical with respect to the stator teeth axis.Hence, both rotors contribute with the same amount of magnetic flux (ϕ Lk ,ϕ Rk ).Therefore, the flux linkage of the kth tooth coil with N c turns is simply given by: Energies 2016, 9, 892 12 of 18

Results and Discussion
In this section, some important machine characteristics, e.g., back-EMF and torque, will be presented and compared.Since these characteristics under different SE conditions only differ in magnitude, the 40% eccentricity is chosen to carry out the analysis.
As shown in Figure 14, the flux passes from one rotor disc and enters the other rotor disc.Therefore, the coils are assumed to be symmetrical with respect to the stator teeth axis.Hence, both rotors contribute with the same amount of magnetic flux (ϕLk,ϕRk).Therefore, the flux linkage of the kth tooth coil with Nc turns is simply given by:

Machine Main Characteristics under No Load Condition
The no-load air gap total flux density under healthy condition compared with 40% eccentricity condition, the analytical results are shown in Figure 15.Comparing with the healthy conditions, it can be seen that the flux density near minimum air gap has increased but the flux density near maximum air gap has decreased.This is owing to the partly increased yet partly decreased air gap length around circumferential direction.
The back EMF and cogging torque are also important for PM machines.For this purpose, the rotor movement should be modeled to calculate these parameters.The back EMF waveform of the AFPM can be calculated from the axial component of the no-load flux density distribution (By) with the knowledge of the armature winding distribution.The voltage induced in a single coil in each layers of quasi-3D method can be calculated as: where Ψc,i is the flux linkage of the coil in ith layer and is equal to the integral of the air gap flux density distribution across one coil pitch, as given by:

Machine Main Characteristics under No Load Condition
The no-load air gap total flux density under healthy condition compared with 40% eccentricity condition, the analytical results are shown in Figure 15.

Results and Discussion
In this section, some important machine characteristics, e.g., back-EMF and torque, will be presented and compared.Since these characteristics under different SE conditions only differ in magnitude, the 40% eccentricity is chosen to carry out the analysis.
As shown in Figure 14, the flux passes from one rotor disc and enters the other rotor disc.Therefore, the coils are assumed to be symmetrical with respect to the stator teeth axis.Hence, both rotors contribute with the same amount of magnetic flux (ϕLk,ϕRk).Therefore, the flux linkage of the kth tooth coil with Nc turns is simply given by:

Machine Main Characteristics under No Load Condition
The no-load air gap total flux density under healthy condition compared with 40% eccentricity condition, the analytical results are shown in Figure 15.Comparing with the healthy conditions, it can be seen that the flux density near minimum air gap has increased but the flux density near maximum air gap has decreased.This is owing to the partly increased yet partly decreased air gap length around circumferential direction.
The back EMF and cogging torque are also important for PM machines.For this purpose, the rotor movement should be modeled to calculate these parameters.The back EMF waveform of the AFPM can be calculated from the axial component of the no-load flux density distribution (By) with the knowledge of the armature winding distribution.The voltage induced in a single coil in each layers of quasi-3D method can be calculated as: where Ψc,i is the flux linkage of the coil in ith layer and is equal to the integral of the air gap flux density distribution across one coil pitch, as given by: Comparing with the healthy conditions, it can be seen that the flux density near minimum air gap has increased but the flux density near maximum air gap has decreased.This is owing to the partly increased yet partly decreased air gap length around circumferential direction.
The back EMF and cogging torque are also important for PM machines.For this purpose, the rotor movement should be modeled to calculate these parameters.The back EMF waveform of the AFPM can be calculated from the axial component of the no-load flux density distribution (B y ) with the knowledge of the armature winding distribution.The voltage induced in a single coil in each layers of quasi-3D method can be calculated as: where Ψ c,i is the flux linkage of the coil in ith layer and is equal to the integral of the air gap flux density distribution across one coil pitch, as given by: where θ 0 is the coil starting side angle from origin and θ c is the angle of coil pitch, and B y,i is the axial component of flux density in ith layer.The back-EMF per phase is calculated in all coils with the phase winding connected in parallel.For phase A, according to the coils arrangement shown in Figure 3, the back-EMF in ith layer E A,i is given by: And based on the quasi-3-D consumption, the back EMF for the whole machine is calculated as: The back-EMF of the prototype AFPM with and without SE condition obtained from both FE and SC methods are compared in Figure 16.
It can be seen that in healthy condition, the induced voltage of three coils have the same shape and amplitude.However, monitoring back EMF is not an effective approach to detect eccentric fault in double rotor AFPMs since the back EMF is not affected by eccentricity sufficiently in SE condition.Eccentricity increases the air gap, but on the other hand, the air gap length of corresponding position of the other rotor is decreased.These two effects seem to compensate each other, given the same results.
Energies 2016, 9, 892 13 of 18 where θ0 is the coil starting side angle from origin and θc is the angle of coil pitch, and By,i is the axial component of flux density in ith layer.
The back-EMF per phase is calculated in all coils with the phase winding connected in parallel.For phase A, according to the coils arrangement shown in Figure 3, the back-EMF in ith layer EA,i is given by: , And based on the quasi-3-D consumption, the back EMF for the whole machine is calculated as: The back-EMF of the prototype AFPM with and without SE condition obtained from both FE and SC methods are compared in Figure 16.
It can be seen that in healthy condition, the induced voltage of three coils have the same shape and amplitude.However, monitoring back EMF is not an effective approach to detect eccentric fault in double rotor AFPMs since the back EMF is not affected by eccentricity sufficiently in SE condition.Eccentricity increases the air gap, but on the other hand, the air gap length of corresponding position of the other rotor is decreased.These two effects seem to compensate each other, given the same results.The back EMFs and its fast Fourier transform (FFT) under SE condition varies from 0% to 65% is shown in Figures 17 and 18 respectively.When the SE factor is increasing, the amplitude of back EMFs is slightly increasing.Moreover, the SE has a strong impact on the first order of FFT.

Back EMF of Phase A [V]
Angular Position (Mech deg.)The back EMFs and its fast Fourier transform (FFT) under SE condition varies from 0% to 65% is shown in Figures 17 and 18 respectively.When the SE factor is increasing, the amplitude of back EMFs is slightly increasing.Moreover, the SE has a strong impact on the first order of FFT.
Energies 2016, 9, 892 13 of 18 where θ0 is the coil starting side angle from origin and θc is the angle of coil pitch, and By,i is the axial component of flux density in ith layer.
The back-EMF per phase is calculated in all coils with the phase winding connected in parallel.For phase A, according to the coils arrangement shown in Figure 3, the back-EMF in ith layer EA,i is given by: , And based on the quasi-3-D consumption, the back EMF for the whole machine is calculated as: The back-EMF of the prototype AFPM with and without SE condition obtained from both FE and SC methods are compared in Figure 16.
It can be seen that in healthy condition, the induced voltage of three coils have the same shape and amplitude.However, monitoring back EMF is not an effective approach to detect eccentric fault in double rotor AFPMs since the back EMF is not affected by eccentricity sufficiently in SE condition.Eccentricity increases the air gap, but on the other hand, the air gap length of corresponding position of the other rotor is decreased.These two effects seem to compensate each other, given the same results.The back EMFs and its fast Fourier transform (FFT) under SE condition varies from 0% to 65% is shown in Figures 17 and 18 respectively.When the SE factor is increasing, the amplitude of back EMFs is slightly increasing.Moreover, the SE has a strong impact on the first order of FFT.

Back EMF of Phase A [V]
Angular Position (Mech deg.)In the case of SE, the following differential flux linkages act on the phase loops [26], in this paper, phase A is disconnected the parallel path (shown in Figure 16) and: where ψ lI and ψ lII are magnetic flux linkages of two disconnected parallel paths of phase A.
By time deriving, the no-load loop EMF follows: Thus, monitoring loop EMF could be a practical approach to detect defects of AFPMs. Figure 19 shows the distribution of phase A loop EMF during one complete rotor revolution at rated speed.It can be seen that the loop EMF fluctuates at each position.
The cogging torque is another important parameter for PM machines which occurs in slotted AFPM machine due to the interaction between the PM magnetic field and the varying air gap permeance.
By assuming negligible saturation, the cogging torque is independent of armature currents and is created only by the PM magnetic field.The period of the cogging torque waveform is calculated as follows: Cogging torque period = 2π/LCM (Qs, 2p) where LCM (Qs, 2p) is the least common multiple of the stator slot number and the rotor pole number.The cogging torque is calculated by the changing rate of total air gap co-energy including the region of PMs [27]: It should be noted that the AFPM machine in this paper has two rotor disks, which means that both rotors contribute to the cogging torque.Thus, the cogging torque of each disk should be superimposed.The cogging torque waveforms obtained under healthy and SE condition are compared in Figure 20.Comparing with the results obtained from FE model, we can see that the results obtained by analytical model do not agree well, this is mainly due to the subtraction in Equation ( 40), but the trend of analytical model is still in accordance with that of FE model.
In order to validate the results, 3D FEM calculation is also depicted, which shows a good agreement with analytical simulations.From the simulation results and comparisons, the cogging torque of each disk has the same shape and amplitude, but under SE condition, the cogging torque is increased.In the case of SE, the following differential flux linkages act on the phase loops [26], in this paper, phase A is disconnected the parallel path (shown in Figure 16) and: where ψ l I and ψ l I I are magnetic flux linkages of two disconnected parallel paths of phase A. By time deriving, the no-load loop EMF follows: Thus, monitoring loop EMF could be a practical approach to detect defects of AFPMs. Figure 19 shows the distribution of phase A loop EMF during one complete rotor revolution at rated speed.It can be seen that the loop EMF fluctuates at each position.
The cogging torque is another important parameter for PM machines which occurs in slotted AFPM machine due to the interaction between the PM magnetic field and the varying air gap permeance.
By assuming negligible saturation, the cogging torque is independent of armature currents and is created only by the PM magnetic field.The period of the cogging torque waveform is calculated as follows: Cogging torque period = 2π/LCM (Q s , 2p) where LCM (Q s , 2p) is the least common multiple of the stator slot number and the rotor pole number.The cogging torque is calculated by the changing rate of total air gap co-energy including the region of PMs [27]: It should be noted that the AFPM machine in this paper has two rotor disks, which means that both rotors contribute to the cogging torque.Thus, the cogging torque of each disk should be superimposed.The cogging torque waveforms obtained under healthy and SE condition are compared in Figure 20.Comparing with the results obtained from FE model, we can see that the results obtained by analytical model do not agree well, this is mainly due to the subtraction in Equation ( 40), but the trend of analytical model is still in accordance with that of FE model.
In order to validate the results, 3D FEM calculation is also depicted, which shows a good agreement with analytical simulations.From the simulation results and comparisons, the cogging torque of each disk has the same shape and amplitude, but under SE condition, the cogging torque is increased.

Machine Main Characteristics under on Load Condition
The output torque is comprised of EM and cogging torque.Neglecting saturation, currents in each phase are updated according to rotor position.The electromagnetic torque for each layer could be calculated by integrating Maxwell stress tensor near the center of the air gap: The electromagnetic torque of machine is obtained as: Also, the average electromagnetic torque can be simply calculated from where kw is the winding factor,φf is the magnetic flux excited by the permanent magnet per pole and I is the stator current.From pervious study, the magnetic flux of each phase can be compensated by

Machine Main Characteristics under on Load Condition
The output torque is comprised of EM and cogging torque.Neglecting saturation, currents in each phase are updated according to rotor position.The electromagnetic torque for each layer could be calculated by integrating Maxwell stress tensor near the center of the air gap: The electromagnetic torque of machine is obtained as: Also, the average electromagnetic torque can be simply calculated from where kw is the winding factor,φf is the magnetic flux excited by the permanent magnet per pole and I is the stator current.From pervious study, the magnetic flux of each phase can be compensated by

Machine Main Characteristics under on Load Condition
The output torque is comprised of EM and cogging torque.Neglecting saturation, currents in each phase are updated according to rotor position.The electromagnetic torque for each layer could be calculated by integrating Maxwell stress tensor near the center of the air gap: The electromagnetic torque of machine is obtained as: Also, the average electromagnetic torque can be simply calculated from where k w is the winding factor, φ f is the magnetic flux excited by the permanent magnet per pole and I is the stator current.From pervious study, the magnetic flux of each phase can be compensated by the right and left magnetic flux, thus, SE does not have significant effects on the EM torque.This also can be proved in Equation ( 46).As shown in Figure 21, the results of the proposed method are in good agreement with the simulation results of 3D-FE model.
Energies 2016, 9, 892 16 of 18 the right and left magnetic flux, thus, SE does not have significant effects on the EM torque.This also can be proved in Equation ( 46).As shown in Figure 21, the results of the proposed method are in good agreement with the simulation results of 3D-FE model.

Experimental Results
The double rotor AFPM in Table 1 is the prototype machine investigated in this paper, as shown in Figure 22a,b.The experimental set-up and devices are presented in Figure 22c.The AFPM motor is driven by a frequency converter shown in Figure 22c to test its motor performance, and the no-load parameters.The back EMF is measured by a driving motor via a belt and the drive speed is 300 rpm, the measured back EMF is shown in Figure 23.The analytical model provides good correspondence with the experimental results shown in Figure 16.Table 2 shows the back EMF coefficient obtained by analytical model, FE model and experiment, respectively.However, there is a slight error between the analytical results and the measured value which is mainly caused by the assembly process.
In terms of the computational time, the analytical model requires only 206 s to obtain the cogging torque and electromagnetic torque with one mechanical cycle for SE condition.The 3D FE model, on the other hand, requires almost 27 h to compute the performance in one electrical cycle.Although the analytical method shows a slight error compared to the experimental result, it is still acceptable and thus, can be regarded as a meaningful approach which could save time and achieve a satisfying result.Moreover, the analytical model could also be used for further calculation, e.g., optimization.

Experimental Results
The double rotor AFPM in Table 1 is the prototype machine investigated in this paper, as shown in Figure 22a,b.The experimental set-up and devices are presented in Figure 22c.The AFPM motor is driven by a frequency converter shown in Figure 22c to test its motor performance, and the no-load parameters.
Energies 2016, 9, 892 16 of 18 the right and left magnetic flux, thus, SE does not have significant effects on the EM torque.This also can be proved in Equation ( 46).As shown in Figure 21, the results of the proposed method are in good agreement with the simulation results of 3D-FE model.

Experimental Results
The double rotor AFPM in Table 1 is the prototype machine investigated in this paper, as shown in Figure 22a,b.The experimental set-up and devices are presented in Figure 22c.The AFPM motor is driven by a frequency converter shown in Figure 22c to test its motor performance, and the no-load parameters.The back EMF is measured by a driving motor via a belt and the drive speed is 300 rpm, the measured back EMF is shown in Figure 23.The analytical model provides good correspondence with the experimental results shown in Figure 16.Table 2 shows the back EMF coefficient obtained by analytical model, FE model and experiment, respectively.However, there is a slight error between the analytical results and the measured value which is mainly caused by the assembly process.
In terms of the computational time, the analytical model requires only 206 s to obtain the cogging torque and electromagnetic torque with one mechanical cycle for SE condition.The 3D FE model, on the other hand, requires almost 27 h to compute the performance in one electrical cycle.Although the analytical method shows a slight error compared to the experimental result, it is still acceptable and thus, can be regarded as a meaningful approach which could save time and achieve a satisfying result.Moreover, the analytical model could also be used for further calculation, e.g., optimization.The back EMF is measured by a driving motor via a belt and the drive speed is 300 rpm, the measured back EMF is shown in Figure 23.The analytical model provides good correspondence with the experimental results shown in Figure 16.Table 2 shows the back EMF coefficient obtained by analytical model, FE model and experiment, respectively.However, there is a slight error between the analytical results and the measured value which is mainly caused by the assembly process.
In terms of the computational time, the analytical model requires only 206 s to obtain the cogging torque and electromagnetic torque with one mechanical cycle for SE condition.The 3D FE model, on the other hand, requires almost 27 h to compute the performance in one electrical cycle.Although the analytical method shows a slight error compared to the experimental result, it is still acceptable and thus, can be regarded as a meaningful approach which could save time and achieve a satisfying result.Moreover, the analytical model could also be used for further calculation, e.g., optimization.

Conclusions
In this paper an approach for modeling double rotor AFPM with SE has been presented.To take into account the intrinsic 3D nature of the AFPM machine, a quasi-3D method is used.The magnetic field, cogging torque and output torque of AFPM machine under healthy and 40% eccentricity conditions are quickly calculated with great accuracy compared to the 3D FE model and experimental results.In addition, it can be found that the eccentricity plays no effects on the back EMF and the electromagnetic torque.However, the cogging torque increases with SE.In addition, this paper makes a remarkable contribution to the computational time savings while maintaining as much the high accuracy the FE model can provide.Finally, an experimental validation has been carried out to verify the accuracy of the simulation results obtained by the proposed method.

Conclusions
In this paper an approach for modeling double rotor AFPM with SE has been presented.To take into account the intrinsic 3D nature of the AFPM machine, a quasi-3D method is used.The magnetic field, cogging torque and output torque of AFPM machine under healthy and 40% eccentricity conditions are quickly calculated with great accuracy compared to the 3D FE model and experimental results.In addition, it can be found that the eccentricity plays no effects on the back EMF and the electromagnetic torque.However, the cogging torque increases with SE.In addition, this paper makes a remarkable contribution to the computational time savings while maintaining as much the high accuracy the FE model can provide.Finally, an experimental validation has been carried out to verify the accuracy of the simulation results obtained by the proposed method.

Figure 1 .
Figure 1.Construction of the investigated axial flux permanent magnet machine (AFPMM) (a) and the principle of quasi-3D method (b-d).

Figure 1 .
Figure 1.Construction of the investigated axial flux permanent magnet machine (AFPMM) (a) and the principle of quasi-3D method (b-d).

Figure 2 .
Figure 2. Representation of field regions divided by magnets.

Figure 3 .
Figure 3. Current sheet spatial distribution in space of each phase in per unit for AFPM with six poles and nine stator slots.

Figure 2 .
Figure 2. Representation of field regions divided by magnets.

Figure 2 .
Figure 2. Representation of field regions divided by magnets.

Figure 3 .
Figure 3. Current sheet spatial distribution in space of each phase in per unit for AFPM with six poles and nine stator slots.

Figure 3 .
Figure 3. Current sheet spatial distribution in space of each phase in per unit for AFPM with six poles and nine stator slots.

Figure 4 .
Figure 4. Representation of field regions caused by armature reaction current sheet.

Figure 4 .
Figure 4. Representation of field regions caused by armature reaction current sheet.

Figure 7 .
Figure 7. Calculation processes of magnet flux density.

Figure 7 .
Figure 7. Calculation processes of magnet flux density.

Figure 8 .
Figure 8. Complex relative permeance without eccentricity in the middle of air gap.(a) Axial component; and (b) Circumferential component.

Figure 9 .Figure 10 .
Figure 9. Axial and circumferential components of permanent magnetic (PM) flux density at the center of air gap under no-load condition.

Figure 8 .
Figure 8. Complex relative permeance without eccentricity in the middle of air gap.(a) Axial component; and (b) Circumferential component.

Figure 8 .
Figure 8. Complex relative permeance without eccentricity in the middle of air gap.(a) Axial component; and (b) Circumferential component.

Figure 9 .Figure 10 .
Figure 9. Axial and circumferential components of permanent magnetic (PM) flux density at the center of air gap under no-load condition.

Figure 9 .
Figure 9. Axial and circumferential components of permanent magnetic (PM) flux density at the center of air gap under no-load condition.

Figure 8 .
Figure 8. Complex relative permeance without eccentricity in the middle of air gap.(a) Axial component; and (b) Circumferential component.

Figure 9 .Figure 10 .
Figure 9. Axial and circumferential components of permanent magnetic (PM) flux density at the center of air gap under no-load condition.

Figure 10 .
Figure 10.Axial and circumferential components of PM flux density at the center of air gap under rated load.Figure 10.Axial and circumferential components of PM flux density at the center of air gap under rated load.

Figure 11 .
Figure 11.The 3D FE model and its meshing.

Figure 12 .Figure 13 .
Figure 12.Axial and circumferential components of flux density at the air gap (l = L − 1 mm) under no-load condition.

Figure 11 .
Figure 11.The 3D FE model and its meshing.

Figure 11 .
Figure 11.The 3D FE model and its meshing.

Figure 12 .Figure 13 .
Figure 12.Axial and circumferential components of flux density at the air gap (l = L − 1 mm) under no-load condition.

Figure 12 .
Figure 12.Axial and circumferential components of flux density at the air gap (l = L − 1 mm) under no-load condition.

Figure 11 .
Figure 11.The 3D FE model and its meshing.

Figure 12 .Figure 13 .
Figure 12.Axial and circumferential components of flux density at the air gap (l = L − 1 mm) under no-load condition.

Figure 13 .
Figure 13.Axial and circumferential components of flux density at the air gap (l = L − 1 mm) under rated-load condition.Figure 13.Axial and circumferential components of flux density at the air gap (l = L − 1 mm) under rated-load condition.

Figure 14 .
Figure 14.Right and left magnetic flux in the tooth.

Figure 15 .
Figure 15.Total flux density comparison at the air gap under no-load condition.

Figure 14 .
Figure 14.Right and left magnetic flux in the tooth.

Figure 14 .
Figure 14.Right and left magnetic flux in the tooth.

Figure 15 .
Figure 15.Total flux density comparison at the air gap under no-load condition.

Figure 15 .
Figure 15.Total flux density comparison at the air gap under no-load condition.

Figure 18 .
Figure 18.The fast Fourier transform (FFT) of back EMF under SE condition varies from 0% to 65%.

Figure 18 .
Figure 18.The fast Fourier transform (FFT) of back EMF under SE condition varies from 0% to 65%.

Figure 19 .
Figure 19.Connections among the phase to form the phase winding (up) and the waveforms of the loop EMFs of phase A during one complete rotor revolution at rated speed (down).

Figure 19 . 18 Figure 19 .
Figure 19.Connections among the phase to form the phase winding (up) and the waveforms of the loop EMFs of phase A during one complete rotor revolution at rated speed (down).

Figure 22 .
Figure 22.Experimental test rig.(a) The permanent magnet rotor; (b) The prototype machine; and (c) The test rig.

Figure 22 .
Figure 22.Experimental test rig.(a) The permanent magnet rotor; (b) The prototype machine; and (c) The test rig.

Figure 22 .
Figure 22.Experimental test rig.(a) The permanent magnet rotor; (b) The prototype machine; and (c) The test rig.

Figure 23 .
Figure 23.Experimental Back EMF waveform of AFPMMs (a) and the comparison with analytical model (b).

Figure 23 .
Figure 23.Experimental Back EMF waveform of AFPMMs (a) and the comparison with analytical model (b).

Table 2 .
Comparison of Back EMF with Experiment Data.

Table 2 .
Comparison of Back EMF with Experiment Data.