Analytical Calculation of Air-Gap Magnetic Field Distribution in IPMSMs with Mixed Eccentricity Accounting for Bridge Saturation

: This paper presents an analytical model for calculating the detailed air-gap magnetic ﬁeld distribution in interior permanent magnet synchronous motors (IPMSMs) with mixed eccentricity. In order to improve the efﬁciency of model solving, a modeling strategy combining the equivalent magnetic circuit network method and the subdomain method was adopted. Speciﬁcally, the magnetic ﬁeld distribution of the rotor was modeled by using the magnetic circuit analysis method, and the magnetic ﬁeld distributions of the air-gap, slot opening and stator slot along the radial direction were modeled in different regions according to their structure ruler. Then, the inﬂuence of bridge saturation was considered. Moreover, based on the analysis of the air-gap geometric structure with mixed eccentricity, a detailed spatiotemporal analytical model of the air-gap magnetic ﬁeld was established, which provides a more accurate description of the mixed eccentricity composed of static and dynamic rotor eccentricities of different severity. The analytical models were compared with the corresponding models established by the ﬁnite element method, which proved the accuracy and validity of the models established in this paper. Finally, some key features related to radial and tangential air-gap ﬂux density were extracted, which can signiﬁcantly reﬂect the characteristics of eccentricities. The main ﬁndings reported in this paper will be of beneﬁt for developing methods for early identiﬁcation and diagnosis of eccentricity faults in IPMSMs.


Introduction
Recently, interior permanent magnet synchronous motors (IPMSMs) have been widely used in electric vehicles and other industrial fields due to their characteristics of high efficiency, high power density, carbon free and zero emission [1][2][3][4]. As a result, the requirements for safe and healthy operation with high performance are becoming higher and higher. In actual application, due to the manufacturing, installation and variation in support stiffness, rotor mixed eccentricity is a widespread phenomenon in the operation of the IPMSM, which leads to the distortion of the air-gap magnetic field of the motor and affects the performance and efficiency of the motor, and even affects the service life of the motor [5]. Therefore, it is important to clarify the variation rule and characteristics of the air-gap magnetic field in IPMSMs with mixed eccentricity, so as to further realize the early diagnosis of mixed eccentricity faults and the real-time control, which is of great significance for the healthy and efficient operation of the motor [5][6][7].
In order to build the model of air-gap magnetic field, two types of methods are usually used: finite element methods and analytical methods. The finite element method is a relatively mature method for numerical modeling. Through the refined solid models and the configuration of constraints, the air-gap magnetic field can be accurately calculated under different operation conditions. However, the real-time performance is not good, In view of the existing problems in the current analytical methods and in order to accurately calculate the air-gap magnetic field distribution in IPMSMs with mixed eccentricity, an analytical model that accounts for the influence of magnetic saturation in the bridge regions was established by combining the equivalent magnetic circuit network method and the subdomain method. This analytical model shows high accuracy, which has been verified by finite element simulation. Then, using the established model, the distribution of the air-gap magnetic field of the IPMSM with mixed eccentricity values of different severity was calculated, and features reflecting the eccentricity were extracted. The results of this study lay a foundation for developing early identification and diagnosis methods for eccentricity faults in IPMSMs.

Methods
The analytical model of air-gap magnetic field distribution in the IPMSM with mixed eccentricity can be established based on the normal model whose air-gap length is replaced by the air-gap geometric model with mixed eccentricity. In this section, the analytical model of the air-gap magnetic field in the IPMSM is established with the following assumptions: (1) The end effect, the rotor core conductivity and the eddy current effect are ignored. (2) The magnetic permeability of the rest of the iron core is infinite except for magnetic saturation in the bridge. (3) The stator slots and slot openings are approximately fanshaped. The normal air-gap magnetic field considering bridge saturation is composed of the no-loaded air-gap magnetic field of the permanent magnet, the armature reaction air-gap magnetic field and the equivalent magnetic field of the effect of bridge saturation. During the establishment of each magnetic field model, the motor is divided into four subdomains along the radial direction: rotor, air-gap, slot opening and stator slot. The air-gap flux density contains temporal and spatial variables and is described with a radial component and a tangential component. The analytical method for efficient calculation of air-gap magnetic field is formed by combining the equivalent magnetic circuit network method and the subdomain method. The flux is produced by the scalar magnetic potential on the rotor surface, and the part entering the air-gap is solved by using the equivalent magnetic circuit network method. The equivalent magnetic circuit is shown in Figure 1. The permanent magnet magnetomotive force F M , equivalent reluctance R M , extreme flux leakage Φ b and interpolar flux leakage Φ i can be expressed by the motor parameters as where B M is the residual magnetic induction intensity of the permanent magnet, h M is the thickness of the permanent magnet, µ 0 is the vacuum permeability, µ r is the relative permeability of the permanent magnet, L M is the width of the permanent magnet, L Z is the axial length of the iron core, B S is the magnetic induction intensity of bridge saturation, h b is the thickness of the end of bridge, and h i is the thickness of the interpolar magnetic in the bridge. The scalar magnetic potential of the rotor surface can be expressed as where Φ gr is the magnetic flux entering the air-gap. It can be calculated as The general solution of vector magnetic potential in air-gap is where A gp is the vector magnetic potential of the permanent magnet with no-load, and A gpk , B gpk , C gpk and D gpk are the Fourier series undetermined coefficients of the air-gap magnetic field of the permanent magnet, R S is the inner radius of the stator, k is the harmonic order of the air-gap magnetic field, α is the circumferential angle of the air-gap space, and r is the length from the center of the radial rotor to the inner boundary of the stator. Then, the radial component and tangential component of the flux density are described by Equation (9) and Equation (10), respectively by substituting Equation (8) into Equations (9) and (10) where B rpg and B tpg are the radial and tangential components of air-gap flux density of permanent magnet with no-load, respectively, and k is the harmonic frequency of air-gap flux density. The Fourier undetermined coefficients can be solved by the interface conditions of each subdomain of the IPMSM under no-load conditions. The interface conditions of each subdomain satisfy the conditions that: the vector magnetic potential and tangential flux density are both continuous on the interface between the stator slot and slot opening; the vector magnetic potential and tangential flux density are also continuous on the interface between the slot opening and air-gap; the scalar magnetic potential and air-gap magnetic flux are continuous on the interface between rotor and air-gap. According to the above interface conditions, the constraint Equations of the undetermined coefficients can be derived as where K 11~K88 are the coefficient matrixes under each boundary condition, respectively, which depend on the motor structure parameters, and Y is determined by the rotor magnetic field. A gpk , B gpk , C gpk and D gpk are the undetermined coefficients of Fourier series of air-gap flux density general solution, respectively. A sopi and B sopi are the undetermined coefficients of flux density Fourier general solution in the slot opening domain, respectively, and i is the ith stator slot. A spi is the undetermined coefficient of the Y Fourier solution in the stator slot region [19,20]. By solving the matrix Equation (13), each undetermined coefficient can be obtained, and substituting them into Equations (11) and (12), the radial and tangential components of the air-gap flux density of the permanent magnet of the IPMSM can be obtained.

Model of the Armature Reaction Magnetic Field
When analyzing the armature reaction magnetic field, there is current excitation in the stator slot winding, and the vector magnetic potential in the stator slot region satisfies the Poisson Equation: The form of general solutions of armature reaction air-gap flux density is the same as that of air-gap flux density under no-load, but the value of undetermined coefficient is different. The interface conditions of each subdomain of armature reaction under load of the IPMSM are as follows: The vector magnetic potential and tangential flux density are continuous on the interface between stator slot domain and slot opening domain; the vector magnetic potential and tangential flux density are continuous on the interface between slot opening domain and air-gap domain; the tangential air-gap flux density is 0 on the interface between air-gap domain and rotor domain. According to the interface conditions, the constraint equations of the undetermined coefficients can be derived and described as where K 12~K77 are the equation coefficient matrixes under each boundary condition, depending on the motor structure parameters, and Y 1~Y4 are determined by the armature current density under load conditions. A gak , B gak , C gak and D gak are the undetermined Fourier series coefficients of the air-gap flux density general solution, respectively. A soai and B soai are the undetermined coefficients of the flux density Fourier general solution in the stator slot opening domain, respectively, and i is the ith stator slot. A sai is the undetermined coefficient of the flux density Fourier general solution in the stator slot region [19,20]. By solving Equation (15), each undetermined coefficient can be obtained. Furthermore, the radial and tangential components of the air-gap flux density of the armature reactive of the IPMSM can be obtained.

Model of the Magnetic Field in the Bridge Regions
The effect of the bridge saturation will affect the armature reactive magnetic field, which is essentially due to the magnetic voltage drop caused by the armature reactive magnetic field passing through the bridge of the rotor. The subdomain method combined with the equivalent magnetic circuit method is used to calculate the magnetic field in the bridge regions. The motor is divided into three subdomains: stator, air-gap and rotor. The interface conditions of each subdomain are as follows: the tangential air-gap flux density at the interface between the air-gap and stator is 0; the scalar magnetic potential at the interface between the air-gap and rotor is continuous; the air-gap flux density at the interface between the rotor domain and air-gap domain is continuous. The constraint equations of the undetermined coefficient are obtained as where K 11~K55 depend on the structural parameters of the motor, and Y depends on the flux leakage in the bridge regions. A gbk , B gbk , C gbk and D gbk are the undetermined Fourier series coefficients of the air-gap magnetic general solution, respectively [19,20]. By solving Equation (16), the undetermined coefficients can be calculated. Further, the radial component and tangential component of the air-gap flux density of the rotor magnetic field in the bridge regions of the IPMSM can be found. B rbg and B tbg can be obtained, whose Fourier expressions are the same as those of Equations (11) and (12). By superimposing the air-gap magnetic field under no-load conditions, armature reaction magnetic field and rotor magnetic field in the bridge regions, the analytical model of the normal air-gap magnetic field can be described as

Geometric Model of the Air-Gap with Mixed Eccentricity
The eccentricity caused in the process of motor manufacturing and assembly is an objective factor that cannot be ignored in motor design optimization and performance analysis. Similarly, both static eccentricity and dynamic eccentricity co-exist at the same time during motor operation. Thus, the accurate geometric model of the air-gap with mixed eccentricity is established in this paper. It is assumed that the boundary line of the rotor cross-section is almost an ideal circle, and the shape and position of the stator and rotor do not change in the axial direction. Moreover, the position of the minimum air-gap relative to the stator during static eccentricity remains unchanged [29], and the position of the minimum air-gap relative to the rotor during dynamic eccentricity remains unchanged too [30].
As shown in Figure 2, when there is mixed eccentricity in the motor, the air-gap length g(α,t) is a function of spatial angle α and time t. R s and R r are motor stator inner diameter and rotor outer diameter, respectively, δ is the eccentric angle, O S is the center of the stator circle, O R is the center of the rotor circle, and O R' is the center of the motion track of the rotor center. d SE = ||O S O R || is the static eccentricity distance, while d DE = ||O S O R ||is the dynamic eccentricity distance. We draw a straight line through point O R to intersect the circle O S and the circle O R at points A and B, respectively. The length of line segment AB is g(α,t), and the angle between the line segment O R A and the polar axis is α, namely, the angle of the air-gap space. The angle between the line segment O R O R and the pole axis is ωt, namely, the rotor rotation angle. According to the trigonometric relationship, Equation (19) can be yielded: Further, the length of O s O R can be yielded as and in ∆A O S O R' According to Equation (23), it can be obtained that It can be found in Figure 2 that the air-gap length can be described as and the ratio of the static eccentricity can be described as while the ratio of dynamic eccentricity can be described as where g 0 is the length of air-gap in normal state. The air-gap length of mixed eccentricity can be obtained by combining Equations (19)-(27) and can be described as Obviously, Equation (28) is a general expression suitable for static eccentricity, dynamic eccentricity, and mixed eccentricity. By substituting g(α,t) calculated from Equation (28) into Equations (16) and (17), the air-gap length becomes a function of α and t. The analytical model reflecting the mixed eccentric air-gap magnetic field is finally obtained.

Model Validation Method
In order to verify the accuracy of the analytical model established in this paper, the finite element model of the IPMSM was established as the benchmark model for comparison, which is shown in Figure 3. In the modeling, the assumptions are as follows: (1) The material is isotropic; (2) the permeability of the material is uniform and does not vary with temperature. According to the actual structure of the IPMSM and the division of permeance units, the flux density distribution of the air-gap is calculated by equivalent magnetic network modeling, and the schematic diagram of the basic modeling process is shown in Figure 3c. Please refer to [31] for more detail on the modeling process.
Several typical geometric structures of the permeance units and the corresponding calculation formula of the permeance are shown in Table 1. In consideration of computational accuracy and speed, two-layer air-gap permeance units are adopted and the mesh thickness of each layer is 0.35 mm.
The main parameters of the IPMSM analyzed in this paper are shown in Table 2. The air-gap magnetic field distributions with and without mixed eccentricity are calculated using established analytical models and are both compared with corresponding results calculated with the finite element models. The radial and tangential components of the air-gap magnetic field are compared separately. The accuracy of the analytical model is characterized by the root mean square percentage error (RMSPE), described as where N d is the number of samples of the radial or tangential component of air-gap magnetic field in one period, S Ai is the ith value of the radial or tangential component sequence calculated with the analytical model, and S Ei is the ith value of the radial or tangential component sequence calculated with the finite element model.   In this section, the accuracy of calculation results of air-gap magnetic field without eccentricity was examined. The spatiotemporal distributions of the radial flux density and tangential flux density of the air-gap magnetic field obtained by analytical method and FEA are shown in Figure 4, where the initial spatial angle of the temporal distribution calculation is 3.7 • in this paper. According to Equation (29), the RMSPEs of the spatial radial component and the spatial tangential component are 0.9% and 0.89%, respectively, while the RMSPEs of the temporal radial component and the temporal tangential component are 1.01% and 1.1%, respectively. This shows that the results calculated by the proposed analytical model are consistent with the results of FEA, which demonstrates the accuracy of the proposed analytical model.

The Eccentric Case
In order to show that the analytical model in this paper can accurately and effectively reflect the phenomenon of eccentricity, the spatiotemporal radial and tangential distributions of the air-gap flux density under the conditions of ε S = 50%, and ε D = 50% were calculated. Figure 5 shows the spatiotemporal distributions of radial and tangential components of the air-gap flux density. Equation (29) was also used to analyze the model error.
The calculation results show that the RMSPEs of the spatial radial and tangential components are 1.02% and 0.97%, respectively. Additionally, RMSPEs of the temporal radial and tangential components are 1.11% and 1.12%, respectively. Therefore, the analytical model of the air-gap magnetic field distribution with mixed eccentricity is consistent with the FEA model, which verifies the accuracy and effectiveness of the proposed modelling method. In addition, from the spatial diagram of the radial components, the peak values of the magnetic flux density of the air-gap are mainly affected by the magnetic bridge of the rotor and are close to 2T. Especially in the case of mixed eccentricity, the tightest of the air-gap can reach a saturation state. Once saturated, the maximum magnetic density reaches 2.07T, which is consistent with the results of finite element analysis.

Influence of the Mixed Eccentricity on the Air-Gap Magnetic Field
In order to investigate the influence of mixed eccentricity on the spatiotemporal distribution of the air-gap magnetic field, Figure 6 shows the three-dimensional spatiotemporal distribution of air-gap flux density under the conditions of ε S = 50% and ε D = 50%. Figure 6a shows the three-dimensional spatiotemporal distribution of radial air-gap flux density, and Figure 6b shows the three-dimensional spatiotemporal distribution of tangential air-gap flux density. The blue curves along the x axis (times) are the curves of the variation in radial and tangential flux density with time at a given space angle, while the blue curves along the y axis (space angle) are the curves of the variation of radial and tangential components of air-gap flux density with space angle at a given moment. It can be seen from the figures that, under different space angles, the temporal distributions of the radial and tangential flux density are completely different, while the spatial distributions of the radial and tangential flux density at different times are relatively consistent. Therefore, in the following parts of the paper, only the spatial distribution of the air-gap magnetic field is analyzed. According to the above calculation, it can be found that mixed eccentricity leads to significant changes in the distribution of the air-gap magnetic field. In order to charac-terize this change, some features representing the radial and tangential components of the magnetic field were extracted, and their variation between different degrees of mixed eccentricity was analyzed.

Analysis of the Radial Component of the Air-Gap Flux Density
The radial un-symmetry (USYr) and the radial total harmonic distortion rate (THDr) are designed for the radial component of the air-gap flux density, which can be described as where B r is the spatial distribution of radial component of air-gap flux density, and B rmax and B rmin are the maximum and minimum root mean square values of all periods of B r . B rn is the root mean square value without eccentricity. B r1 is the amplitude of fundamental component of the radial air-gap flux density, and B r (k) is the amplitude of kth harmonic component of radial air-gap flux density.
In order to study the influence of different degrees of static and dynamic eccentricities on radial air-gap flux density, the results of two situations are shown in Table 3; ε S = 20% and tε S = 50% for the first situation, while ε D = 50% and ε D = 20% for the second situation. From the table, we can find that the spatial distribution of the air-gap radial flux density is distorted by mixed eccentricity. Compared with the case without mixed eccentricity, the values of USYr and THDr increase when the mixed eccentricity occurs. In the situation of ε S = 50% and ε D = 20%, the USYr increases by 0.19%, the amplitude of spatial harmonics increases by 0.09% averagely, and the THDr also increases by 0.02%. It can be seen that the influence of static eccentricity of the rotor on air-gap radial flux density is greater than that of dynamic eccentricity. Table 3. The spatial harmonic amplitudes, USYr and THDr with different mixed eccentricities. In order to further highlight the influence of mixed eccentricity on the features of radial air-gap flux density, the degrees and steps of static and dynamic eccentricities were set to be the same, so the degree of mixed eccentricity is the same as the degree of static eccentricity or dynamic eccentricity for the analysis of the variation characteristic of Br with a mixed-eccentric fault. The air-gap flux density was calculated when the static eccentricity angle was 2 • , and the result is shown in Figure 7. The values of USYr and THDr of the radial air-gap flux density also increase with the increase in the degree of mixed eccentricity. When the degree of the mixed eccentricity is low (ε m < 40%), the values of USYr are low and change smoothly, while when the ε m > 40%, the values of USYr increase significantly.

Analysis of the Tangential Air-Gap Flux Density Component
We define the tangential un-symmetry (USYt) of the tangential air-gap flux density component as where B t is the spatial distribution of tangential component of air-gap flux density, and B tmax and B tmin are the maximum and minimum root mean square values of all periods of B t . B tn is the root mean square value without eccentricity. When mixed eccentricity occurs, the tangential air-gap flux density curve will be distorted and become a non-stationary curve. Thus, spectrum analysis based on Fourier transform cannot accurately describe the frequency and amplitude characteristics of tangential air-gap flux density under different spatial angles. In this paper, the wavelet transform method was used to decompose the original tangential air-gap flux density, and the number of decomposition layers was 7. Among them, the low-frequency components a6 and a7 were closer to the envelope of the original signal. As the eccentricity increases, the distortion becomes serious, that is, the asymmetry tends to be serious. The high-frequency components d1, d2, and d3 are obvious periodic signals, and as the degree of eccentricity increases, the signal energy distributions converge toward the center of the samples, as shown in Figure 8.  Therefore, the feature of energy proportion is defined to describe the degree of aggregation of signal energy, which is expressed as where i = 1,2 are the coefficients of the first layer and the second layer of the high-frequency components of the signal, respectively, signal samples are equally divided into 10 intervals, and E Mi is the signal energy of the fifth and sixth intervals in the middle. E T is the total energy of signal samples. The signal energy is calculated by where N d is the calculated sample length and s is the calculated sample values. The samples of low-frequency component a7 were used to calculate the USYt, and the samples of high-frequency components d1, d2, and d3 were used to calculate the E oi . The results are shown in Figure 9. In Figure 9a, USYt increases significantly with the increase in degree of mixed eccentricity, indicating that mixed eccentricity has a serious impact on the low-frequency component of tangential air-gas flux density. In Figure 9b, E oi also increases correspondingly. The E oi of the d1 component is most sensitive to the severity of mixed eccentricity.

Influence of Bridge Saturation on Air-Gap Magnetic Field with Mixed Eccentricity
In order to analyze the influence of bridge saturation on the air-gap magnetic field distribution with mixed eccentricity, the radial spatial component of air-gap flux density was calculated with the analytical model without considering saturation. Specifically, the calculation was performed by ignoring Equation (16) while solving Equations (17) and (18). Then, the radial un-symmetry USYnr and the radial total harmonic distortion rate THDnr were further calculated for the case without considering bridge saturation and were then compared with the results calculated in Section 3.2.1 in which the bridge saturation was considered, as shown in Figure 10. It can be seen from Figure 10 that the values of USYnr are lower than those of USYr when the degree of mixed eccentricities are low. In other words, when the analytical model does not consider the bridge saturation, the calculation results show that the eccentricity has little effect on the magnetic field. So, it is difficult to identify the characteristics of the early stage of eccentricity fault. When the degree of mixed eccentricity is high, the increase in USYnr is greater than that of USYr. That is, it seems that the impact on the magnetic field is greater than in reality when the model without considering the bridge saturation, which will cause misjudgment of the degree of eccentricity. In addition, we found that the bridge saturation has little effect on THDnr.

Conclusions
Mixed eccentricity is a common fault of IPMSMs, which will distort the air-gap magnetic field and affect the motor performance. In this paper, an analytical model of the air-gap magnetic field distribution in IPMSMs with mixed eccentricity considering the effect of bridge saturation is established, and the accuracy is verified by comparison with finite element analysis. On this basis, features of the air-gap magnetic field distribution that can effectively reflect the mixed eccentricity fault are extracted. The calculation results show that features change significantly with the variation in the degree of mixed eccentricity, which can accurately reflect the evolution process and characteristics of the mixed eccentricity. In the future, we will further study the relationship between the features extracted in this paper and the monitored quantities such as electrical parameters and mechanical parameters and achieve early detection and diagnosis of mixed eccentricity.