Analytical Modeling and Comparison of Two Consequent-Pole Magnetic-Geared Machines for Hybrid Electric Vehicles

: The exact mathematical modeling of electric machines has always been an e ﬀ ective tool for scholars to understand the working principles and structure requirements of novel machine topologies. This paper provides an analytical modeling method—the harmonic modeling method (HMM)—for two types of consequent-pole magnetic-geared machines, namely the single consequent-pole magnetic-geared machine (SCP-MGM) and the dual consequent-pole magnetic-geared machine (DCP-MGM). By dividing the whole machine domain into di ﬀ erent ring-like subdomains and solving the Maxwell equations, the magnetic ﬁeld distribution and electromagnetic parameters of the two machines can be obtained, respectively. The two machines were applied in the propulsion systems of hybrid electric vehicles (HEVs). The electromagnetic performances of two machines under di ﬀ erent operating conditions were also compared. It turns out that the DCP-MGM can reach a larger electromagnetic torque compared to that of the SCP-MGM under the same conditions. Finally, the predicted results were veriﬁed by the ﬁnite element analysis (FEA). A good agreement can be observed between HMM and FEA. Furthermore, HMM can also be applied to the mathematical modeling of other consequent-pole electric machines in further study.


Introduction
The last decade has witnessed rapid developments of magnetic gears (MGs) and electric machines that utilize the magnetic-gearing effect, which are also called magnetic-geared machines (MGMs) [1][2][3]. Ever since their invention in 2001 [4], MGs have become a research hotspot due to their high efficiency and self-protection characteristics [5][6][7].
The concept of MGMs is derived from MGs. By substituting stator windings with AC current for one rotating permanent magnet (PM) component, MGMs change one mechanical port of MGs into an electrical port. Thus, the two rotating components of the MGM and its stator windings can be regarded as a combination of a magnetic gear and an electric machine [8][9][10]. Indeed, with the introduction of another rotating component and the ability to alternate the speed ratio and torque ratio between two rotating components, MGMs have broadened the application scenarios of electric machines [11][12][13]. A good example is that MGMs can serve as the power split component (PSC) in hybrid electric vehicles (HEVs) to realize energy exchange among the internal combustion engine (ICE), wheels, and battery [14][15][16]. The ICE and electric machine can provide traction for the wheels rotor of CP-MGM to a synchronous state (the rated rotating speed) at startup state. In addition, it can deliver extra output torque to the outer rotor shaft if the output torque of the CP-MGM cannot meet the requirement. Since this paper mainly focuses on the operating modes of the CP-MGM, it is reasonable to assume that there is no power flow between the PMSM and the outer rotor shaft at the four steady states mentioned in this paper. In fact, power exchange between the PMSM and the wheels would not affect the conclusion obtained in this paper.
The working principle of MGMs is similar to that of magnetic gears. By adopting a modulator layer, the magnetic field distribution can be changed. Assuming that the pole pair number of the original magnetic field is P i , and the modulator number is Q, then a novel magnetic field will have a component that has (Q − P i ) pole pairs. Thus, the fundamental structural requirement of an MGM is [4]: where P s is the pole pair number of stator windings. Under steady working conditions, the rotating speed of two rotating rotors and the current frequency f within stator windings should then satisfy: where ω i , ω o , and ω s are the rotating speed of the inner rotor, outer rotor, and the equivalent rotating speed of stator windings. MGM cannot meet the requirement. Since this paper mainly focuses on the operating modes of the CP-MGM, it is reasonable to assume that there is no power flow between the PMSM and the outer rotor shaft at the four steady states mentioned in this paper. In fact, power exchange between the PMSM and the wheels would not affect the conclusion obtained in this paper. The working principle of MGMs is similar to that of magnetic gears. By adopting a modulator layer, the magnetic field distribution can be changed. Assuming that the pole pair number of the original magnetic field is Pi, and the modulator number is Q, then a novel magnetic field will have a component that has (Q − Pi) pole pairs. Thus, the fundamental structural requirement of an MGM is [4]: where Ps is the pole pair number of stator windings. Under steady working conditions, the rotating speed of two rotating rotors and the current frequency f within stator windings should then satisfy: where ωi, ωo, and ωs are the rotating speed of the inner rotor, outer rotor, and the equivalent rotating speed of stator windings. Since the ICE reaches its highest efficiency at the range of ~2000 r/min-3000 r/min, the rotating speed of outer rotor and the current frequency of stator winding must cooperate with the rotating speed of the inner rotor to ensure the highest efficiency of the ICE. However, if the stator windings need to provide energy for the HEV, the rotating speed of the inner rotor must be smaller than that of the outer rotor. Thus, a gearbox must come into service under hybrid mode to reduce the rotating speed of the inner rotor. Therefore, the operation modes of the proposed HEV propulsion system can be divided into four kinds, and their typical operating parameters are shown in Table 1. The rotating speed of the outer rotor is calculated according to the different driving speeds of the HEV, and the current frequency of stator winding is obtained via Equation (2). The topologies of SCP-MGM and DCP-MGM are shown in Figure 2. By substituting SMM for PMs with the same polarity, a consequent-pole structure is obtained. The name "consequent-pole" is due to SMM, and PM appears alternately on the circumferential direction. Although SMM cannot Since the ICE reaches its highest efficiency at the range of~2000 r/min-3000 r/min, the rotating speed of outer rotor and the current frequency of stator winding must cooperate with the rotating speed of the inner rotor to ensure the highest efficiency of the ICE. However, if the stator windings need to provide energy for the HEV, the rotating speed of the inner rotor must be smaller than that of the outer rotor. Thus, a gearbox must come into service under hybrid mode to reduce the rotating speed of the inner rotor. Therefore, the operation modes of the proposed HEV propulsion system can be divided into four kinds, and their typical operating parameters are shown in Table 1. The rotating speed of the outer rotor is calculated according to the different driving speeds of the HEV, and the current frequency of stator winding is obtained via Equation (2). The topologies of SCP-MGM and DCP-MGM are shown in Figure 2. By substituting SMM for PMs with the same polarity, a consequent-pole structure is obtained. The name "consequent-pole" is due to SMM, and PM appears alternately on the circumferential direction. Although SMM cannot generate a magnetic field itself, it can be easily magnetized to conduct flux lines. Hence, SMM in a consequent-pole structure can be regarded as a magnetic source to some degree. The greatest advantage of using the consequent-pole structure is saving PM material, which is the most expensive material in an electric machine. Both SCP-MGM and DCP-MGM utilize a consequent-pole structure to save PM material. The SMM part in the outer rotor of a DCP-MGM not only works as a consequent-pole structure for the PMs inserted in the outer rotor, it also modulates the magnetic field of the inner rotor. Thus, the P s -th harmonic component within the DCP-MGM is larger than that of the SCP-MGM. Additionally, the saturation of the DCP-MGM is more severe than that of the SCP-MGM. generate a magnetic field itself, it can be easily magnetized to conduct flux lines. Hence, SMM in a consequent-pole structure can be regarded as a magnetic source to some degree. The greatest advantage of using the consequent-pole structure is saving PM material, which is the most expensive material in an electric machine. Both SCP-MGM and DCP-MGM utilize a consequent-pole structure to save PM material. The SMM part in the outer rotor of a DCP-MGM not only works as a consequentpole structure for the PMs inserted in the outer rotor, it also modulates the magnetic field of the inner rotor. Thus, the Ps-th harmonic component within the DCP-MGM is larger than that of the SCP-MGM. Additionally, the saturation of the DCP-MGM is more severe than that of the SCP-MGM.

Assumption and Parameter Definition
The machine structure chosen to be studied in this paper was a 24 slot 11 pole-pair SCP-MGM and DCP-MGM, as shown in Figure 2. A few assumptions must be made to simplify the mathematical modeling:

•
The geometrical shape of the machine has a radial side and a tangential side; • The magnetic field distribution is constrained in the 2D plane: the axial component is ignored;

•
The machine has infinite axial length, so the end effect is ignored; • The radial component of the permeability of SMM within a certain region is regarded as a constant; • Eddy-current effects within SMM and PMs are ignored.
Since there exists a z-direction current within the windings of the studied machines and the machine topology has a circular shape, a vector magnetic potential (VMP) Az in a polar coordinate is adopted to calculate the magnetic flux density distribution within the machines. The machine structures are then divided into several ring-like regions based on the different material interfaces, as can be seen in Figure 3, where α represents the angle of the inner PM arc, β represents the angle of the slot opening in the modulator, δ is the slot opening angle, and γ is the stator slot angle. The whole machine is divided into ten subdomains: the innermost one (region I) represents the shaft part; region II is the rotor yoke; region III is the inner consequent-pole PM; region IV is the inner air gap; region V is the modulator pieces (it should be noted that for SCP-MGMs, the gap between each two modulator pieces is air, while for DCP-MGM, bipolar PMs are inserted in that gap). Region VI is the outer air gap; region VII is the stator teeth; region VIII is the stator slots together with windings; region IX is the stator yoke; region X is the outside of the studied machines.
The angular position of the j-th PM part of the inner rotor θPM, the position of the k-th modulator piece θMod, the position of the t-th stator tooth part θtooth, and the position of the s-th stator slot part θslot can be defined, respectively, as:

Assumption and Parameter Definition
The machine structure chosen to be studied in this paper was a 24 slot 11 pole-pair SCP-MGM and DCP-MGM, as shown in Figure 2. A few assumptions must be made to simplify the mathematical modeling: • The geometrical shape of the machine has a radial side and a tangential side; • The magnetic field distribution is constrained in the 2D plane: the axial component is ignored;

•
The machine has infinite axial length, so the end effect is ignored; • The radial component of the permeability of SMM within a certain region is regarded as a constant; • Eddy-current effects within SMM and PMs are ignored.
Since there exists a z-direction current within the windings of the studied machines and the machine topology has a circular shape, a vector magnetic potential (VMP) A z in a polar coordinate is adopted to calculate the magnetic flux density distribution within the machines. The machine structures are then divided into several ring-like regions based on the different material interfaces, as can be seen in Figure 3, where α represents the angle of the inner PM arc, β represents the angle of the slot opening in the modulator, δ is the slot opening angle, and γ is the stator slot angle. The whole machine is divided into ten subdomains: the innermost one (region I) represents the shaft part; region II is the rotor yoke; region III is the inner consequent-pole PM; region IV is the inner air gap; region V is the modulator pieces (it should be noted that for SCP-MGMs, the gap between each two modulator pieces is air, while for DCP-MGM, bipolar PMs are inserted in that gap). Region VI is the outer air gap; Energies 2019, 12, 1888 5 of 25 region VII is the stator teeth; region VIII is the stator slots together with windings; region IX is the stator yoke; region X is the outside of the studied machines.
The angular position of the j-th PM part of the inner rotor θ PM , the position of the k-th modulator piece θ Mod , the position of the t-th stator tooth part θ tooth , and the position of the s-th stator slot part θ slot can be defined, respectively, as: where ϕ in , and θ 0 the initial angular positions of the inner rotor and outer rotor, respectively. Due to the symmetrical structure of the inner rotor and outer rotor, ϕ in , θ 0 has a range of [0, 2π/P i ], [0, 2π/Q], respectively. Specifically, ϕ in is defined as zero when the lower edge of the PM in the inner rotor coincides with the positive direction of the angular axis; θ 0 is defined as zero when the center of the slot of the outer rotor (air in SCP-MGM and PM in DCP-MGM) coincides with the positive direction of the angular axis, as shown in Figure 2. The VMP A  , the magnetic flux density B  , the magnetic field strength  H , and the current density distribution  J in stator windings can be written in vector form as: To simplify the solving process, all the parameters related to magnetic field are expressed in terms of complex Fourier series. Thus, the vector amplitude above can be further expressed as: To simplify the solving process, all the parameters related to magnetic field are expressed in terms of complex Fourier series. Thus, the vector amplitude above can be further expressed as: where n represents the n-th order coefficient of the corresponding Fourier series. It should be noticed that, in numerical calculation, a reasonable harmonic order N is used to truncate the infinite Fourier series. If N is too small, the Fourier series will have a large error, if N is too large, the calculation time will be rather long.

Partical Differential Equation Solution
The magnetic field within the machine follows quasistatic Maxwell equations: The relationship between → B and → A can be further expressed as: The radial component and tangential component matrix of magnetic flux density → B are then obtained in matrix form as [31]: where K represents the harmonic order coefficient diagonal matrix that is related to N, given by: Similar to Equations (11)- (14), the relative permeability of each region can also be expressed in a complex Fourier series form: Next, based on the relation between → B and → H, as expressed below: where → M is the magnetization vector. The first item on the right is a product of two Fourier series, which can be rewritten in matrix form by using the Cauchy product theorem: where µ r,cov and µ θ,cov are convolution matrices of the radial and tangential components of permeability, respectively. M r and M θ are the radial and tangential components of magnetization intensity, respectively. M r and M θ can all written in complex Fourier series. The convolution matrix µ r,cov can be defined as: whereμ n is the n-th order coefficient of the Fourier series of corresponding µ.
In Equation (24), H θ is continuous at the interface between two regions, but B θ is discontinuous at the interface. Hence, the matrix µ θ,cov cannot be settled using Equation (25). Instead, a fast Fourier factorization is applied to calculate µ θ,cov , for the sake of keeping the rate of convergence the same for the left and right side of Equation (24) [34]. µ θ,cov can be given by [31]: From Equation (26), it can be seen that µ θ,cov is acquired by replacingμ n with the corresponding n-th order Fourier coefficient of 1/µ θ for each element, and there is a matrix inversion outside.
The region V in SCP-MGM is used to illustrated the convolution matrix with respect to relative permeability, as shown in Figure 4.
where μ n is the n-th order coefficient of the Fourier series of corresponding μ.
In Equation (24), Hθ is continuous at the interface between two regions, but Bθ is discontinuous at the interface. Hence, the matrix μθ,cov cannot be settled using Equation (25). Instead, a fast Fourier factorization is applied to calculate μθ,cov, for the sake of keeping the rate of convergence the same for the left and right side of Equation (24) [34]. μθ,cov can be given by [31]: From Equation (26), it can be seen that μθ,cov is acquired by replacing μ n with the corresponding n-th order Fourier coefficient of 1 θ μ / for each element, and there is a matrix inversion outside.
The region V in SCP-MGM is used to illustrated the convolution matrix with respect to relative permeability, as shown in Figure 4.
Air Air Air Air The relative permeability distribution on the circumferential direction can be expressed as: When a Fourier expansion on [0, 2π] is applied to Equation (27), the expressions of μ k and μ rec k are given by:  The relative permeability distribution on the circumferential direction can be expressed as: where ϕ out = θ 0 − β 2 . When a Fourier expansion on [0, 2π] is applied to Equation (27), the expressions ofμ k andμ rec k are given by:μ The convolution matrices µ r,cov and µ θ,cov can then be obtained by substituting Equations (28) and (29) into (25) and (26).
Equation (30) is a Cauchy-Euler differential equation system [35]. The general solution of a single differential equation in Equation (30) is given by: where c 1 , c 2 are unknown coefficients. Similarly, the general solution of the differential equation system in Equation (30) can be written in a matrix form, where the new element (i, j) in matrix r V is defined as: Hence, the complementary solution of Equation (30) A k z com can be written as: where C 1 and C 2 are unknown coefficient vectors. Matrix r V 1 2 can be factorized as: where λ is the eigenvalue matrix of V 1/2 , and matrix P is the eigenvector matrix of V 1/2 . Therefore, Equation (33) can be simplified as: As for the particular solution of Equation (30), A k z par is given by: The general solution of VMP in Equation (30) is the sum of the complementary solution of Equation (35), and particular solution Equation (36): The expression of VMP in each subdomain is given in Appendix A.
According to the geometrical characteristic of SCP-MGMs and DCP-MGMs, the magnetization intensity only exists in region III and region V, and the current density distribution only exists in region VIII. Their distribution waveforms and Fourier series coefficients can be seen in Table 1. Additionally, Energies 2019, 12, 1888 9 of 25 for regions I, II, IV, VI, IX, and X, there only exists one material type, so the coefficients within the convolution matrix are a constant. However, the permeability distributions are different in region III, V, VII, and VIII, as shown in Table 2; their coefficients can be obtained by substituting a, b, and c in Table 3 into the following equations: Additionally, for regions I, II, IV, VI, IX, and X, there only exists one material type, so the coefficients within the convolution matrix are a constant. However, the permeability distributions are different in region III, V, VII, and VIII, as shown in Table 2; their coefficients can be obtained by substituting a, b, and c in Table 3 into the following equations:

Sources Illustrative Waveforms Fourier Series Coefficients
Inner PM (Region III) Outer PM (Only for DCP-MGM) (Region V) Additionally, for regions I, II, IV, VI, IX, and X, there only exists one material type, so the coefficients within the convolution matrix are a constant. However, the permeability distributions are different in region III, V, VII, and VIII, as shown in Table 2; their coefficients can be obtained by substituting a, b, and c in Table 3 into the following equations:

Sources Illustrative Waveforms Fourier Series Coefficients
Inner PM (Region III) ) n 0

Stator windings (Region VIII)
Energies 2019, 12, x FOR PEER REVIEW 9 of 25 Additionally, for regions I, II, IV, VI, IX, and X, there only exists one material type, so the coefficients within the convolution matrix are a constant. However, the permeability distributions are different in region III, V, VII, and VIII, as shown in Table 2; their coefficients can be obtained by substituting a, b, and c in Table 3 into the following equations:

Sources Illustrative Waveforms Fourier Series Coefficients
Inner PM (Region III)   Additionally, for regions I, II, IV, VI, IX, and X, there only exists one material type, so the coefficients within the convolution matrix are a constant. However, the permeability distributions are different in region III, V, VII, and VIII, as shown in Table 2; their coefficients can be obtained by substituting a, b, and c in Table 3 into the following equations:

Sources Illustrative Waveforms Fourier Series Coefficients
Inner PM (Region III) where k represents the k-th subdomain of the proposed machine, and 2 ≤ k ≤ 9. Suppose the subdomain number of the proposed machine is L, and the harmonic order to be calculated is N. By applying the above boundary conditions to each interface of the machines, a system of 2*(L − 1) linear equations with 2*(L − 1) unknowns can be obtained, and each unknown is an (N × 1) vector. The system of linear equations can be further written in matrix form, as below: where the expressions of S, X, and T in this paper are given in Appendix A.
As long as the coefficient matrix S is invertible, the unknown vector X can be acquired. The numerical solution of Equation (42) can be obtained in the MATLAB software.

Saturation Consideration of Soft-Magnetic Material
For the SMM in the consequent-pole part, modulator, and stator teeth, the flux line is concentrated, thus, the saturation of the SMM must be considered. The nonlinear B-µ curve of 50JN1300 is given in Figure 4. In HMM, the relative permeability µ is obtained by an iterative method, as shown by the dot lines in Figure 5, where the number "1, 2, 3, 4" means the iteration number. The detailed iterative process is shown in Figure 6. First, an initial value, namely µ 0 = 1500, is assigned to a given region with SMM for the first iteration. In the k-th iteration, the average flux density → B of a specific region can be obtained by substituting µ i,k and solving the matrix Equation (42). A new average relative permeability µ i,cal in region i can then be acquired. Where i belongs to {III, V, VII, VIII}. The average relative permeability µ i,k+1 for the (k + 1)-th iteration in region i is given as: The iteration will stop only if the iteration time n i exceeds the maximum number of iterations N i (N i is set to be 50 here), or the maximum error ∆ in all these regions is below the error requirement ξ, ∆ is defined as: where ξ is set to be 0.05 in this paper. The saturations of SMM in SCP-MGMs and DCP-MGMs are calculated respectively using this method. The relative permeabilities of each region under on-load conditions are shown in Table 4. It can be observed that the saturation of regions V and VII is more severe in DCP-MGMs, since there are PMs inserted in the modulator, thus, there are more flux lines in the modulator and stator teeth.

Electromagnetic Parameters Calculation
Once the VMP is solved in Equation (42), the related electromagnetic parameters of the two machines can be calculated. The electromagnetic torque of the two rotors of SCP-MGMs and DCP-MGMs can be calculated by using the Maxwell stress tensor.
The flux linkage of each coil side can be given by [36]:

Electromagnetic Parameters Calculation
Once the VMP is solved in Equation (42), the related electromagnetic parameters of the two machines can be calculated. The electromagnetic torque of the two rotors of SCP-MGMs and DCP-MGMs can be calculated by using the Maxwell stress tensor.
The flux linkage of each coil side can be given by [36]:

Electromagnetic Parameters Calculation
Once the VMP is solved in Equation (42), the related electromagnetic parameters of the two machines can be calculated. The electromagnetic torque of the two rotors of SCP-MGMs and DCP-MGMs can be calculated by using the Maxwell stress tensor. The flux linkage of each coil side can be given by [36]: where N turn is the number of coil turns in each slot, and S coil is the cross section area of a single slot. When all the coils in each phase are in series, the three-phase flux linkage can be written as: where C turn is coil-connecting matrix of the proposed machine, the coil number is given in Figure 2, and C turn can be expressed as: The back electromotive force (EMF) is computed by the derivative of ψ with respect to time: where ω is the rotating speed of the magnetic field in the outer air gap. The electromagnetic torque is calculated by using a Maxwell stress tensor. Thus, the electromagnetic torque of the inner rotor T in equals the calculus of the Maxwell stress tensor of the inner air gap along the circumferential direction, and the electromagnetic torque of the outer rotor T out equals the algebraic sum of the calculus of the Maxwell stress tensor of both the inner air gap and the outer air gap along the circumferential direction. They can be expressed as: where R i and R o are the middle radius of the inner air gap and outer air gap, respectively.

Simulation Environment and Machine Parameters
To make quantitative comparison between the SCP-MGMs and DCP-MGMs, all the geometrical parameters of these two machines should be set as the same, and other parameters, such as the slot filling factor and root mean square value of the winding current should be also set as the same, as shown in Table 5. The analytical prediction of the HMM was carried out using MATLAB, and the FEA model was constructed and run in JMAG software. The FEA model had 35,597 elements and 25,368 nodes; the element size near the air gap was set as 1 mm, to maintain calculation accuracy. It took 4.6 s for the computer to obtain the magnetic field distribution for one step. The computer system configuration was as follows: Processor: Intel Core i7-4790 CPU @ 3.60 GHz; Installed Memory (RAM):  28.0 GB-System type: 64-bit Windows Operating System. Additionally, the mean error percentage ε 1 and maximum difference ε 2 of the two methods was defined as: where V HMM,n is the n-th value, calculated using the HMM, and V FEA,n is n-th the value obtained by FEA. N c is the total number of calculated points.

Comparison between HMM and FEA
Based on the calculation and simulation, a good agreement between the HMM and FEA can be observed in Figures 7-12 for SCP-MGMs and DCP-MGMs. It was also found that the difference between HMM and FEA at a no-load condition was less compared to that at a load condition. This is due to the truncation of the infinite Fourier series; the high-order components take up a greater proportion at a no-load condition, leading to a larger error for HMM. Additionally, the numerical values calculated via MATLAB were constrained by the computational accuracy. Values exceeding the computational accuracy were ignored, which led to errors of the magnetic field calculation. Generally, if the harmonic order and computational accuracy is improved, the error will decrease. However, the computational time increases rapidly with the increase of harmonic order and computational accuracy. Thus, there is a tradeoff between calculation accuracy and calculation time. In this paper, the computational accuracy of MATLAB was set as 32 bits.
The mean error percentage ε 1 and maximum difference ε 2 for the SCP-MGM and DCP-MGM are listed in Table 6. where VHMM,n is the n-th value, calculated using the HMM, and VFEA,n is n-th the value obtained by FEA. Nc is the total number of calculated points.

Comparison between HMM and FEA
Based on the calculation and simulation, a good agreement between the HMM and FEA can be observed in Figures 7-12 for SCP-MGMs and DCP-MGMs. It was also found that the difference between HMM and FEA at a no-load condition was less compared to that at a load condition. This is due to the truncation of the infinite Fourier series; the high-order components take up a greater proportion at a no-load condition, leading to a larger error for HMM. Additionally, the numerical values calculated via MATLAB were constrained by the computational accuracy. Values exceeding the computational accuracy were ignored, which led to errors of the magnetic field calculation. Generally, if the harmonic order and computational accuracy is improved, the error will decrease. However, the computational time increases rapidly with the increase of harmonic order and computational accuracy. Thus, there is a tradeoff between calculation accuracy and calculation time. In this paper, the computational accuracy of MATLAB was set as 32 bits.      The mean error percentage ε1 and maximum difference ε2 for the SCP-MGM and DCP-MGM are listed in Table 6.      Figure 13. It can be seen that the second harmonic component was much higher for the DCP-MGM, since the inserted PM on the outer rotor produced a second magnetic field after the modulation of the consequent-pole iron part of the inner rotor. Hence, the consequent-pole structure of the inner rotor and outer rotor of the DCP-MGM could modulate the magnetic field generated by its counterpart. The electromagnetic torque of the DCP-MGM was expected to be larger than that of the SCP-MGM under the same working conditions. Additionally, the difference of each frequency component between the HMM and FEA was very small. Specifically, the error percentages of the second harmonic component for the SCP-MGM and the DCP-MGM using HMM and FEA were 4.14% and 2.15%, respectively. A fast Fourier transform (FFT) was executed on the no-load radial component of the magnetic flux density of the outer air gap for both the SCP-MGM and DCP-MGM, as shown in Figure 13. It can be seen that the second harmonic component was much higher for the DCP-MGM, since the inserted PM on the outer rotor produced a second magnetic field after the modulation of the consequent-pole iron part of the inner rotor. Hence, the consequent-pole structure of the inner rotor and outer rotor of the DCP-MGM could modulate the magnetic field generated by its counterpart. The electromagnetic torque of the DCP-MGM was expected to be larger than that of the SCP-MGM under the same working conditions. Additionally, the difference of each frequency component between the HMM and FEA was very small. Specifically, the error percentages of the second harmonic component for the SCP-MGM and the DCP-MGM using HMM and FEA were 4.14% and 2.15%, respectively.

Electromagnetic Performance Analysis under Different Working Conditions
The electromagnetic performances of the two MGMs are predicted by HMM and simulated in FEA software under different operating conditions, where the rotating speed of two rotors and the current frequency of stator windings are set to cooperate with the practical driving conditions of the HEC, as shown in Table 1. Since the ICE is connected to the inner rotor to provide the power; the outer rotor is connected to the differential, which is further connected to the wheels; the battery is connected to the stator windings, and there can be two-way power flow between the battery and the stator windings. Thus, the power transmission relation among inner rotor, outer rotor and stator windings is equal to the power transmission relation among ICE, wheels and battery. Assume the anti-clock direction is positive, the torque is positive if it's on the anti-clock direction, otherwise it's negative. The power of a component is defined as a positive one when this component is inputting energy; when a component is outputting energy, its power is defined as a negative one.

Back EMF under No-Load Condition
The amplitude of back EMF was determined by the rotating speed of both the inner rotor and the outer rotor. Figure 14 shows the no-load back EMF of the SCP-MGM and the DCP-MGM under the same operating conditions, namely ωi = 1200 r/min, ωo = 1500 r/min. It can be seen that there was an error between the predicted back EMF and FEA result; because the back EMF was calculated as the derivative of flux linkage with respect to time, a small error of flux linkage will be amplified on the back EMF. Additionally, the back EMF of the DCP-MGM was larger than that of the SCP-MGM, due to the flux enhancing effect of the PMs inserted into the modulator. The maximum errors for the SCP-MGM and DCP-MGM were 165 V and 520 V, respectively. The average errors of the SCP-MGM and DCP-MGM were 3.7% and 6.9%, respectively.

Electromagnetic Performance Analysis under Different Working Conditions
The electromagnetic performances of the two MGMs are predicted by HMM and simulated in FEA software under different operating conditions, where the rotating speed of two rotors and the current frequency of stator windings are set to cooperate with the practical driving conditions of the HEC, as shown in Table 1. Since the ICE is connected to the inner rotor to provide the power; the outer rotor is connected to the differential, which is further connected to the wheels; the battery is connected to the stator windings, and there can be two-way power flow between the battery and the stator windings. Thus, the power transmission relation among inner rotor, outer rotor and stator windings is equal to the power transmission relation among ICE, wheels and battery. Assume the anti-clock direction is positive, the torque is positive if it's on the anti-clock direction, otherwise it's negative. The power of a component is defined as a positive one when this component is inputting energy; when a component is outputting energy, its power is defined as a negative one.

Back EMF under No-Load Condition
The amplitude of back EMF was determined by the rotating speed of both the inner rotor and the outer rotor. Figure 14 shows the no-load back EMF of the SCP-MGM and the DCP-MGM under the same operating conditions, namely ω i = 1200 r/min, ω o = 1500 r/min. It can be seen that there was an error between the predicted back EMF and FEA result; because the back EMF was calculated as the derivative of flux linkage with respect to time, a small error of flux linkage will be amplified on the back EMF. Additionally, the back EMF of the DCP-MGM was larger than that of the SCP-MGM, due to the flux enhancing effect of the PMs inserted into the modulator. The maximum errors for the SCP-MGM and DCP-MGM were 165 V and 520 V, respectively. The average errors of the SCP-MGM and DCP-MGM were 3.7% and 6.9%, respectively. Under this mode, the HEV is at a low driving speed, the inner rotor is locked, and the ICE does not come into service; stator windings only work to provide the power that the HEV needs. Thus, ωi = 0 r/min, ωo = 500 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 15 and 16.  Under this mode, the HEV is at a low driving speed, the inner rotor is locked, and the ICE does not come into service; stator windings only work to provide the power that the HEV needs. Thus, ω i = 0 r/min, ω o = 500 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 15 and 16. Under this mode, the HEV is at a low driving speed, the inner rotor is locked, and the ICE does not come into service; stator windings only work to provide the power that the HEV needs. Thus, ωi = 0 r/min, ωo = 500 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 15 and 16.  Under this mode, the HEV is running at a medium speed, so the battery does not output power anymore and the ICE comes into use. However, to maintain the magnetic field, the stator windings are electrified with DC current [14]. The rotating speeds of two rotors are: ω i = 1200 r/min, ω o = 1015 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 17 and 18. Under this mode, the HEV is running at a medium speed, so the battery does not output power anymore and the ICE comes into use. However, to maintain the magnetic field, the stator windings are electrified with DC current [14]. The rotating speeds of two rotors are: ωi = 1200 r/min, ωo = 1015 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 17 and 18. When the HEV needs to further accelerate, the ICE alone is not enough to provide the power that the HEV needs. The SCP-MGM and DCP-MGM can then switch into hybrid mode. In this mode, both ICE and battery provide the energy for the wheels. The rotating speeds of two rotors are: ωi = Under this mode, the HEV is running at a medium speed, so the battery does not output power anymore and the ICE comes into use. However, to maintain the magnetic field, the stator windings are electrified with DC current [14]. The rotating speeds of two rotors are: ωi = 1200 r/min, ωo = 1015 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 17 and 18. When the HEV needs to further accelerate, the ICE alone is not enough to provide the power that the HEV needs. The SCP-MGM and DCP-MGM can then switch into hybrid mode. In this mode, both ICE and battery provide the energy for the wheels. The rotating speeds of two rotors are: ωi = 1200 r/min, ωo = 2000 r/min. The torque waveforms of the power distribution of the two machines When the HEV needs to further accelerate, the ICE alone is not enough to provide the power that the HEV needs. The SCP-MGM and DCP-MGM can then switch into hybrid mode. In this mode, both ICE and battery provide the energy for the wheels. The rotating speeds of two rotors are: ω i = 1200 r/min, ω o = 2000 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 19 and 20. When the HEV deaccelerates, the ICE stops working. The magnetic field generated via the stator windings provides a resistance for the wheels, and thus the power flows from the wheels to the battery. At this time, the power flows from the stator windings to the battery. Thus, the battery can be charged under this mode. The rotating speeds of two rotors are: ωi = 0 r/min, ωo = 1000 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 21 and 22. When the HEV deaccelerates, the ICE stops working. The magnetic field generated via the stator windings provides a resistance for the wheels, and thus the power flows from the wheels to the battery. At this time, the power flows from the stator windings to the battery. Thus, the battery can be charged under this mode. The rotating speeds of two rotors are: ω i = 0 r/min, ω o = 1000 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 21 and 22.

Regenerative Braking Mode (Mode 4)
When the HEV deaccelerates, the ICE stops working. The magnetic field generated via the stator windings provides a resistance for the wheels, and thus the power flows from the wheels to the battery. At this time, the power flows from the stator windings to the battery. Thus, the battery can be charged under this mode. The rotating speeds of two rotors are: ωi = 0 r/min, ωo = 1000 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 21 and 22.  Table 7 shows the mean error percentage and maximum difference of the torque waveforms between FEA and HMM of the SCP-MGM and DCP-MGM under different operating modes. It can be observed that the mean error percentage was below 8%, which is acceptable for the electromagnetic torque prediction of SCP-MGMs and DCP-MGMs. The error of torque prediction was mainly caused by the error of the air-gap magnetic flux density prediction on both the radial and tangential direction, since the torque was calculated using the calculus of the product of the radial component and tangential component of air-gap magnetic flux density.
From the above comparisons between the SCP-MGM and DCP-MGM under various working conditions, it can be seen that the torque ratio between outer rotor and inner rotor was about 1.18, which is equal to the pole-pair ratio of outer rotor to inner rotor. The ripple rate of the outer rotor was larger than that of the inner rotor, which can be alleviated by adopting a skewed stator when these machines are applied in practical applications. Additionally, the electromagnetic torque of the DCP-MGM was about 1.88 times of that of SCP-MGM. However, the torque per unit weight of the PM of the DCP-MGM was only about 0.72 of that of SCP-MGM. Considering the fact that the consequent-pole structure already saves half the PM material that would otherwise have been used, the total usage of PM material for DCP-MGM is reasonable. Thus, the DCP-MGM is more suitable to be used in the propulsion systems of HEVs.   Table 7 shows the mean error percentage and maximum difference of the torque waveforms between FEA and HMM of the SCP-MGM and DCP-MGM under different operating modes. It can be observed that the mean error percentage was below 8%, which is acceptable for the electromagnetic torque prediction of SCP-MGMs and DCP-MGMs. The error of torque prediction was mainly caused by the error of the air-gap magnetic flux density prediction on both the radial and tangential direction, since the torque was calculated using the calculus of the product of the radial component and tangential component of air-gap magnetic flux density. From the above comparisons between the SCP-MGM and DCP-MGM under various working conditions, it can be seen that the torque ratio between outer rotor and inner rotor was about 1.18, which is equal to the pole-pair ratio of outer rotor to inner rotor. The ripple rate of the outer rotor was larger than that of the inner rotor, which can be alleviated by adopting a skewed stator when these machines are applied in practical applications. Additionally, the electromagnetic torque of the DCP-MGM was about 1.88 times of that of SCP-MGM. However, the torque per unit weight of the PM of the DCP-MGM was only about 0.72 of that of SCP-MGM. Considering the fact that the consequent-pole structure already saves half the PM material that would otherwise have been used, the total usage of PM material for DCP-MGM is reasonable. Thus, the DCP-MGM is more suitable to be used in the propulsion systems of HEVs.

Discussion of HMM
The greatest advantage of HMM is that the modeling process is simpler compared to the conventional subdomain method [30]. Since it unitizes Fourier series of relative permeability, the boundary conditions become simpler and the number of unknowns becomes less. For instance, the consequent-pole structure leads to a very complex general solution if the subdomain method is used [29]. Additionally, magnetic saturation can be taken into account to an extent in HMM, and a magnetic flux density distribution figure can be obtained. However, the HMM also has some drawbacks. First, HMM can only provide a simplified magnetic saturation model, since the magnetic saturation varies from point to point in some SMM parts. The relative permeability in the radial direction is also regarded as a constant, but, in reality, it could change. To the authors' knowledge, no existing analytical method can really reflect the magnetic saturation at every point within an electric machine. Secondly, the use of convolution sometimes leads to an extremely large or small value, which further leads to a large error in the final result due to the accuracy limit of the numerical calculation.

Conclusions
In this paper, an analytical modeling method for consequent-pole MGMs was proposed and elaborated. Two machine topologies, namely the SCP-MGM and DCP-MGM, were analyzed and compared quantitatively using HMM and FEA. A good agreement was achieved for these two methods. Furthermore, these two machines were embedded into the propulsion system of HEVs under different operating conditions. By inserting extra PMs in the modulator, the electromagnetic torque of the DCP-MGM increased greatly compared to its counterpart. Therefore, the DCP-MGM has the potential to be applied in the propulsion systems of HEVs. Additionally, the proposed HMM can also be applied in the mathematical modeling of other consequent-pole electric machines.  The expressions of X, S, T are: K 9,8 K 9,9 K 9,10 K 9,11 = P V P V R 4 |K| −I (A15) K 10,8 K 10,9 K 10,10 K 10,11 = µ VI P V λ V −µ VI P V λ V R 4