Analysis of an IPMSM Hybrid Magnetic Equivalent Circuit

: The most common type of electric vehicle traction motor is the interior permanent magnet synchronous motor (IPMSM). For IPMSM designs, engineers make use of the magnetic equivalent circuit method, which is a lumped constant circuit method, and the ﬁnite element method, which is a distributed constant circuit method. The magnetic equivalent circuit method is useful for simple design through fast and intuitive parameters, but it cannot derive the distribution of the magnetic ﬁeld. The ﬁnite element method can derive an accurate magnetic ﬁeld distribution, but it takes a long time and is difﬁcult to use for analysis of intuitive design parameters. In this study, the magnetic equivalent circuit method and Carter’s coefﬁcient were combined for rotor structure design and accurate identiﬁcation and analysis of circuit constants. In this paper, this design method is called the hybrid magnetic equivalent circuit method. Intuitive design parameters are derived through this hybrid magnetic equivalent circuit method. The air gap ﬂux density distribution according to rotor shape, no-load-induced voltage, and cogging torque was analyzed and compared to results of the ﬁnite element method. The proposed method was found to achieve a short solving time and acceptably accurate results. Author Contributions: Conceptualization, I.-S.S.; Data curation, I.-S.S. and B.-W.J.; Formal analysis, I.-S.S.; Investigation, B.-W.J.; Supervision, K.-C.K.; Writing—original draft, B.-W.J.; Writing—review& editing, K.-C.K.


Introduction
Recently, the demand for electric vehicles (EVs) has increased worldwide in response to eco-friendly policies and stricter emission regulations. The traction motor for most EVs is the interior permanent magnet synchronous motor (IPMSM) [1][2][3][4]. The IPMSM field does not require an excitation current because a permanent magnet is used. Thus, the rotor loss is lower than that of other motors and high-efficiency operation is possible. In addition, the output density is higher than that of other motors because magnetic torque using the arranging force of the field and armature occurs at the same time as the reluctance torque using the salient polarity of the rotor. Furthermore, a wide operating range is made possible by its field-weakening operation characteristics. Because of these characteristics, an IPMSM has multiple ratings and is suitable as a traction motor for EVs, which require maximum distance traveled on one battery charge (or range). At the same speed, high torque in the motor indicates high output; thus, torque is a parameter that determines the EV's range. On the other hand, torque ripple is a source of electrical noise and vibration and therefore is a factor that affects driver and passenger riding comfort [5][6][7][8].
IPMSMs have two operation areas. One is a constant-torque range, and the other is a constant-output operation range. The change in this range is quite diverse depending on the motor design method [9][10][11]. To design EV IPMSMs, engineers mainly use the magnetic equivalent circuit method, which is a lumped constant circuit method, and the finite element method (FEM), which is a distributed constant circuit method. In the magnetic equivalent circuit, the maximum speed, demagnetization details, and maximum torque are simply determined by obtaining the flux per pole, d-axis inductance, and q-axis inductance. After the simple design, the final design is made in consideration of torque ripple, cogging torque, and voltage containing harmonic components with respect to the spatial and temporal magnetic field changes through the FEM, the spatial harmonic method, and other methods.
The magnetic equivalent circuit method is useful for simple design through fast and intuitive parameters but cannot derive the magnetic field distribution. The FEM can derive an accurate magnetic field distribution, but it takes a long time, and it is difficult to analyze intuitive design parameters. Many studies have analyzed the magnetic field distribution in the air gap using an analytical circuit model to predict accuracy results [12][13][14]. To provide intuitive solutions in extremely short times, field quantities such as leakage flux paths and saturation, which are important in iron loss calculations, are evaluated using analytical circuit methods [15][16][17][18][19][20]. In addition, this analytical method can be applied to the analysis of many devices, such as transformers, DC/DC converter-coupled inductors, switched reluctance motors, and interior permanent magnet motors where saturation and losses are important [21].
The FEM is used to verify the theory and the expertise of the designer and must be used properly. In other words, the advantages and limitations of existing design techniques are clear [22][23][24][25].
In this study, the magnetic equivalent circuit method and Carter's coefficient are mixed for rotor structure design. Carter's coefficient helps to obtain accurate results based on defined parameters. This proposed method will be called the hybrid magnetic equivalent circuit (HMEC). It is a design method that combines the spatial harmonic method and the magnetic circuit method. It can provide accurate solution results in a shorter time than FEM, and intuitive design parameters are derived through the proposed HMEC method. Additionally, it will represent the air gap flux density distribution according to rotor shape, no-load-induced voltage, and cogging torque. In the conclusion, the proposed HMEC method is compared with the FEM method and verified that the obtained results are appropriate. Figure 1 shows three regions of the slotless flux line distribution by permanent magnet. θ p , θ a , and θ b represent, respectively, the pole arc angle, the angle between the starting point of the bridge from the center of the q-axis, and the angle between the bridge start point and end point. The slotless flux line distribution by the permanent magnet has a constant value in the region of the pole arc, and the magnetic field is close to zero in the region of the q-axis magnetic path and falls constantly in the region of a bridge. Therefore, if it is treated as a line function that falls constantly in a barrier region, the function of the air flux density by permanent magnet is derived as Equation (1) and shown in Figure 2 (the slotless air gap flux density distribution). If the Fourier series is applied to the function of Equation (1), Equations (2) and (3) are derived. The average air gap flux density is derived through the magnetic equivalent circuit.

Calculation of the Slotless Air Gap Flux Density Distribution by Permanent Magnet
where B g , n, α, B gn and N p are, respectively, the average flux density of the air gap, harmonic order, angle of rotor position, Fourier coefficient of flux density, and number of poles.
where g B , n , α , gn B and p N are, respectively, the average flux density of the air gap, harmonic order, angle of rotor position, Fourier coefficient of flux density, and number of poles.   Figure 3 shows the magnetic circuit of a slotless permanent magnet. The reluctance of the rotor core is quite small compared with that of the air gap, barrier, and permanent magnet, so it is neglected. Because the reluctance of the stator core is connected in series with reluctance of the air gap, as shown in Equation (4), it is treated as r K . As shown in Equation (5), the barrier leakage of the permanent magnet is treated as l K . By applying the flux distributive law, the flux equation of the air gap is derived as Equation (6). If Equation (6) is divided by the area of the air gap, the average of the air gap flux density of Equation (7) is derived.  Figure 3 shows the magnetic circuit of a slotless permanent magnet. The reluctance of the rotor core is quite small compared with that of the air gap, barrier, and permanent magnet, so it is neglected. Because the reluctance of the stator core is connected in series with reluctance of the air gap, as shown in Equation (4), it is treated as K r . As shown in Equation (5), the barrier leakage of the permanent magnet is treated as K l . By applying the flux distributive law, the flux equation of the air gap is derived as Equation (6). If Equation (6) is divided by the area of the air gap, the average of the air gap flux density of Equation (7) is derived.
Energies 2021, 14, 5011 4 of 17 where R m , R bl , R b , R g , R sy , φ g and φ r are, respectively, the reluctance of the permanent magnet, reluctance of a barrier, reluctance of a bridge, reluctance of the air gap, reluctance of the stator yoke, flux of the air gap, and residual flux.

Calculation of the Air Gap Flux Density Distribution by Permanent Magnet with Slot
The structure of a motor with slots, as shown in Figure 4, is common. Figure 4 shows the flux line distribution by permanent magnets with slots. When there is a slot, the flux can be divided into two regions. One is the flux region that goes directly from the pole to the teeth, and the other is the flux region that goes from the pole to the teeth through the slot opening. The magnetic field is constant in the teeth region, but as can be seen from the flux in the slot opening region, the length of the air gap increases toward the center of the slot opening and reaches its maximum length at the center. The distribution of the slotopening region is derived by applying the air gap length function of Carter's coefficient. In addition, the total average air gap flux density, as shown in Equation (8), is derived through a magnetic equivalent circuit, as shown in Figure 5.  are, respectively, the reluctance of the rotor teeth, reluctance of the rotor yoke, permeance of the air gap corresponding to slot k, total permeance, area of the permanent magnet, area of air gap corresponding to

Calculation of the Air Gap Flux Density Distribution by Permanent Magnet with Slot
The structure of a motor with slots, as shown in Figure 4, is common. Figure 4 shows the flux line distribution by permanent magnets with slots. When there is a slot, the flux can be divided into two regions. One is the flux region that goes directly from the pole to the teeth, and the other is the flux region that goes from the pole to the teeth through the slot opening. The magnetic field is constant in the teeth region, but as can be seen from the flux in the slot opening region, the length of the air gap increases toward the center of the slot opening and reaches its maximum length at the center. The distribution of the slot-opening region is derived by applying the air gap length function of Carter's coefficient. In addition, the total average air gap flux density, as shown in Equation (8), is derived through a magnetic equivalent circuit, as shown in Figure 5. where R rt , R ry , P gk , P total , A m , A gk , A gtotal , φ gtotal and B gtotal are, respectively, the reluctance of the rotor teeth, reluctance of the rotor yoke, permeance of the air gap corresponding to slot k, total permeance, area of the permanent magnet, area of air gap corresponding to slot k, area of the total air gap, total flux of the air gap, and total average air gap flux density.   Figure 6 shows the air gap flux line distribution per slot when there are slots. The average air gap flux density-derived Equation (12) is the total average of both the slotopening and the teeth region. A transformation of the equation is necessary to obtain an accurate distribution. First, the air gap flux density in the teeth region, as in Equation (14), can be obtained by multiplying Equations (12) and (13), which is the flux distributive law of the reluctance difference of the teeth and the slot-opening regions. Finally, the total air gap flux density in the teeth region is derived as Equation (15).
are a correction factor for teeth about the air gap flux density, permeance of the air gap corresponding to the stator teeth of slot   Figure 6 shows the air gap flux line distribution per slot when there are slots. The average air gap flux density-derived Equation (12) is the total average of both the slotopening and the teeth region. A transformation of the equation is necessary to obtain an accurate distribution. First, the air gap flux density in the teeth region, as in Equation (14), can be obtained by multiplying Equations (12) and (13), which is the flux distributive law of the reluctance difference of the teeth and the slot-opening regions. Finally, the total air gap flux density in the teeth region is derived as Equation (15) are a correction factor for teeth about the air gap flux density, permeance of the air gap corresponding to the stator teeth of slot  Figure 6 shows the air gap flux line distribution per slot when there are slots. The average air gap flux density-derived Equation (12) is the total average of both the slotopening and the teeth region. A transformation of the equation is necessary to obtain an accurate distribution. First, the air gap flux density in the teeth region, as in Equation (14), can be obtained by multiplying Equations (12) and (13), which is the flux distributive law of the reluctance difference of the teeth and the slot-opening regions. Finally, the total air gap flux density in the teeth region is derived as Equation (15). where C tk , P gtk , 2P gso , A gtk , A gttotal , φ gttotal and B gttotal are a correction factor for teeth about the air gap flux density, permeance of the air gap corresponding to the stator teeth of slot k, permeance of the air gap corresponding to the slot opening of slot k, area of the stator teeth, total area of the stator teeth, total flux of the stator teeth, and total flux density of the stator teeth.
Energies 2021, 14, x FOR PEER REVIEW 7 of 20 teeth, total area of the stator teeth, total flux of the stator teeth, and total flux density of the stator teeth. The function of the air gap length of the slot-opening region is expressed as Equation (16). In order to derive the function of the slot factor, the magnetic circuit equation obtained by functionalizing the air gap length is divided into the magnetic circuit equation before functionalization. If these equations are rearranged, the function equation of the slot factor is derived as Equation (17). However, this slot factor function does not reflect the saturation of flux density of the shoe part, as shown in Figure 7. Therefore, by applying the slotopening range of Hanselman, the slot opening is changed to 0.7, which is larger than the actual 0.5, and the function is changed to Equation (18) [14]. The final slot factor function is derived as Equations (19) and (20). Applying this expression to the Fourier series yields Equations (21)- (23). Figure 6b shows the slot factor derived through the HMEC and FEM.
( )   The function of the air gap length of the slot-opening region is expressed as Equation (16). In order to derive the function of the slot factor, the magnetic circuit equation obtained by functionalizing the air gap length is divided into the magnetic circuit equation before functionalization. If these equations are rearranged, the function equation of the slot factor is derived as Equation (17). However, this slot factor function does not reflect the saturation of flux density of the shoe part, as shown in Figure 7. Therefore, by applying the slot-opening range of Hanselman, the slot opening is changed to 0.7, which is larger than the actual 0.5, and the function is changed to Equation (18) [14]. The final slot factor function is derived as Equations (19) and (20). Applying this expression to the Fourier series yields Equations (21)- (23). Figure 6b shows the slot factor derived through the HMEC and FEM. where G sl (θ), G t , θ sp , θ so , θ tr , r si , G min and N s are a function of slot factor, stator teeth pitch, slot pitch, slot-opening pitch, stator teeth pitch considering saturation of the shoe, inner radius of the stator, minimum value of the slot factor, and the number of slots.
where ( ) sl G θ , t G , sp θ , so θ , tr θ , si r , min G and s N are a function of slot factor, stator teeth pitch, slot pitch, slot-opening pitch, stator teeth pitch considering saturation of the shoe, inner radius of the stator, minimum value of the slot factor, and the number of slots.
Energies 2021, 14, x FOR PEER REVIEW 9 of 20  Figure 9 shows the flux linkage derived through the HEMC and the FEM. Equation (25) represents the flux linkage of one coil. In the case of this model, the number of slots per phase per pole is three. Therefore, the flux linkage in one phase can be obtained as in Equation (26). By differentiating Equation (26), the no-load-induced voltage can be derived as Equation (27). Figure 10 shows a comparison of the no-load-induced voltage derived from the HMEC and the FEM. The cogging torque (Equation (29)) is derived by differentiating the air gap energy of Equation (28). Figure 11 shows a comparison of the cogging torque waveform derived from the HMEC and FEM.  Figure 9 shows the flux linkage derived through the HEMC and the FEM. Equation (25) represents the flux linkage of one coil. In the case of this model, the number of slots per phase per pole is three. Therefore, the flux linkage in one phase can be obtained as in Equation (26). By differentiating Equation (26), the no-load-induced voltage can be derived as Equation (27). Figure 10 shows a comparison of the no-load-induced voltage derived from the HMEC and the FEM. The cogging torque (Equation (29)) is derived by  Figure 11 shows a comparison of the cogging torque waveform derived from the HMEC and FEM.

Calculation Process of Flux Linkage, No-Load-Induced Voltage, and Cogging Torque
where N coil , L stk , λ coil and λ ph are the turns of one coil, stack length, flux linkage of one coil, and flux linkage of one phase.         Figure 12 shows a magnetic equivalent circuit of a V-type rotor. As shown in Figure 12, unlike the bar shape, the V shape has an inner bridge and barrier part on the central axis of the pole arc. In addition, by adjusting the angle of the magnet, the area of the permanent magnet can be increased compared with the bar shape. Figure 13 shows a comparison of flux density distribution of the air gap with a bar-type and a V-type rotor with slots. Because the area of the air gap and the permanent magnet is the same, as shown in Figure 13, the difference of the two rotors is the difference in leakage flux according to the presence or absence of the inner bridge and barrier. Therefore, if the area of the permanent magnet is not increased in the V shape, the flux per pole drops as shown in Equations (30) and (31) because of the inner leakage flux. In addition, as shown in Equation (32), as the leakage of the bridge and barrier change, the value decreases and the cogging torque changes, as shown in Figure 14. Figure 14 shows a comparison of cogging torque for a bar-type and a V-type rotor with slots. Table 1 shows comparison of magnetic field distribution parameters of bar-type and V-type rotor when there is a slot.

Calculation of Air Gap Flux Density Distribution by Permanent Magnet V-Type Rotor
where R il , R ol , R ib and R ob are reluctance of the inner and outer barrier, and inner and outer bridge.
where il R , ol R , ib R and ob R are reluctance of the inner and outer barrier, and inner and outer bridge.    (a) (b)     Figure 15 shows the air gap flux density distribution of a slotless double-layer bar-type rotor. In the double-layer bar-type rotor, the air gap region corresponding to the permanent magnet should be divided into two parts.

Calculation of Flux Density Distribution of the Air Gap by Permanent Magnet Double-Layer Type Rotor
Equation (33) shows the magnetic flux density distribution by permanent magnet. Applying this expression to the Fourier series yields Equation (34). Although the flux per pole drops when the same permanent magnet is used, the harmonic order is reduced because the air gap flux density can be made closer to sinusoidal. The double-layer rotor can be classified into three types. As shown in Figure 15, one is a double-layer bar (DB) structure, the second is a double-layer V (DV) structure, and the last is a delta (mixed bar and V-type) structure. The number of permanent magnets used, area of the air gap, and area of the permanent magnet are the same. Therefore, the amount of flux per pole in the two regions varies depending on the presence of the inner barrier and bridge. By applying these conditions, the flux density of the air gap for each two-layer shape in Table 2 was calculated. Figure 16 shows the magnetic equivalent circuit of a permanent magnet when there is no slot for each two-layer rotor. Table 2 shows equations of parameters of air gap flux density distribution according to double-layer rotor type. Because the reluctance of the rotor and stator core were ignored and there was an effective area difference of the air gap between the HMEC and FEM, the value calculated by the HMEC had an error ranging from 2.54% to 3.53% compared with that of the FEM, as shown in Table 3, which shows a comparison of the magnetic field distribution parameters of double layer rotors when there is a slot. In addition, the flux density of the air gap in the DB without the inner bridge and barrier was the highest, and that of the DV was the lowest. The delta shape had a medium flux per pole, but there was no barrier or bridge in the two-layer part, so the inner magnetic flux was large, and the outer magnetic flux was rather small compared with those of the double-layered V-shape. The comparison of flux density distribution when there is a slot for each two-layer shape derived through the HMEC and FEM is represented in Figure 17. As shown in Table 4 between the HMEC and FEM, the value calculated by the HMEC had an error ranging from 2.54% to 3.53% compared with that of the FEM, as shown in Table 3, which shows a comparison of the magnetic field distribution parameters of double layer rotors when there is a slot. In addition, the flux density of the air gap in the DB without the inner bridge and barrier was the highest, and that of the DV was the lowest. The delta shape had a medium flux per pole, but there was no barrier or bridge in the two-layer part, so the inner magnetic flux was large, and the outer magnetic flux was rather small compared with those of the double-layered V-shape. The comparison of flux density distribution when there is a slot for each two-layer shape derived through the HMEC and FEM is represented in Figure 17. As shown in Table 4  between the HMEC and FEM, the value calculated by the HMEC had an error ranging from 2.54% to 3.53% compared with that of the FEM, as shown in Table 3, which shows a comparison of the magnetic field distribution parameters of double layer rotors when there is a slot. In addition, the flux density of the air gap in the DB without the inner bridge and barrier was the highest, and that of the DV was the lowest. The delta shape had a medium flux per pole, but there was no barrier or bridge in the two-layer part, so the inner magnetic flux was large, and the outer magnetic flux was rather small compared with those of the double-layered V-shape. The comparison of flux density distribution when there is a slot for each two-layer shape derived through the HMEC and FEM is represented in Figure 17. As shown in Table 4

Conclusions
In this paper, a hybrid magnetic equivalent circuit method combined with Carter's coefficient is proposed. It is a design method that combines the spatial harmonic method and the magnetic circuit method. The magnetic field distribution according to the rotor type was derived through the hybrid magnetic equivalent circuit method, and the cogging torque and no-load counter electromotive force were derived through the magnetic field distribution map. Through the hybrid magnetic equivalent circuit method, the existing design method was enhanced with reduced analysis time and intuitive design parameters. However, an error occurred because of the difference in the effective air gap area, as a result of ignoring the iron core resistance of the rotor and stator and the effect of bridge magnetic saturation. However, because this error was less than 5% in error range compared with the FEM, the proposed hybrid magnetic equivalent circuit method is considered valid.

Conclusions
In this paper, a hybrid magnetic equivalent circuit method combined with Carter's coefficient is proposed. It is a design method that combines the spatial harmonic method and the magnetic circuit method. The magnetic field distribution according to the rotor type was derived through the hybrid magnetic equivalent circuit method, and the cogging torque and no-load counter electromotive force were derived through the magnetic field distribution map. Through the hybrid magnetic equivalent circuit method, the existing design method was enhanced with reduced analysis time and intuitive design parameters. However, an error occurred because of the difference in the effective air gap area, as a result of ignoring the iron core resistance of the rotor and stator and the effect of bridge magnetic saturation. However, because this error was less than 5% in error range compared with the FEM, the proposed hybrid magnetic equivalent circuit method is considered valid.

Conflicts of Interest:
The authors declare no conflict of interest.