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

Abstract

The exact mathematical modeling of electric machines has always been an effective 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 different ring-like subdomains and solving the Maxwell equations, the magnetic field 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 different 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 verified by the finite 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.

1. 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 independently. The electric machine can work as a generator and a motor. When the electric machine serves as a generator, it can absorb power from the ICE or wheels (depending on working modes) to get the battery charged. When the electric machine serves as a motor, the power flows from the battery to the electric machine to drive the wheels. Hence, the ICE can always work at its highest efficiency to save fuel by alternating the working modes of the electric machine. This application scenario has drawn more and more attention as environmental problems become severe [17]. HEVs do offer a chance to alleviate the exhaust gas emission problem caused by fuel vehicles [18]. Moreover, compared to its counterpart, namely the mechanical gearbox with an electric machine, MGMs not only save space, but also improve efficiency and reduce noise and vibration by eliminating the physical contact of two gear sets [19,20].

Just like permanent magnet synchronous machines (PMSMs), MGMs utilize permanent magnets as the magnetic sources instead of using the electrical excitation method. Thus, the carbon brush structure can be eliminated and the durability of electric machines can be enhanced. However, the rare earth elements make the price of PMs extremely expensive [21]. To solve this problem, a consequent-pole structure can be adopted. The consequent-pole structure can not only reduce the flux linkage, but also improve the torque density [22,23,24]. Two different topologies of consequent-pole MGMs, i.e., single consequent-pole magnetic-geared machines (SCP-MGM) and dual consequent-pole magnetic-geared machines (DCP-MGM) have been proposed [25], but their mathematical modeling has not been well studied.

Although the MGMs offer many new possibilities for electric machines, their magnetic field distribution is much more complex compared to traditional electric machines with one rotor. Many scholars have focused on the magnetic field distribution calculation of MGMs [26,27,28]. Yet, to the best of author’s knowledge, no literature has studied the magnetic field distribution of consequent-pole MGMs. The introduction of soft magnetic material (SMM) to replace the PM part will make the magnetic field distribution of consequent-pole MGMs even more complicated. Research [29] has solved the magnetic field distribution of a PMSM with PMs inserted into the SMM part, but did not consider of the saturation of SMM. Additionally, the subdomain division method [30] is not suitable for MGMs, since too many subdomains increase the calculation time rapidly. Research [31] has proposed a new harmonic modeling method (HMM) to calculate the magnetic distribution of electric machines. By introducing complex Fourier series and a convolution matrix of permeability, HMM can reduce the number of subdomains to within ten. This is because the total number of these ring-like subdomains will not increase with the increase of modulator pieces and slots.

In this paper, two consequent-pole MGMs were studied using HMM. The paper is organized as follows. Section 2 discusses the configurations and operating principles of consequent-pole MGMs. Mathematical models of SCP-MGM and DCP-MGM considering iron saturation are then proposed and elaborated in Section 3. Finally, the effectiveness of proposed HMM is validated by using finite element analysis (FEA) in Section 4.

2. Configurations and Operating Principles of SCP-MGM and DCP-MGM

When the MGM (either SCP-MGM or DCP-MGM) is applied in HEV, its inner rotor can be connected to the ICE, while the outer rotor can be connected to a permanent magnet synchronous machine (PMSM), which will be further connected to the differential to drive the wheels; the battery provides energy to the windings of both the MGM and PMSM via an inverter. The whole system configuration can be seen in Figure 1. The MGM together with the PMSM can be regarded as the E-CVT in a Toyota Prius. They can cooperate with each other according to different working conditions of HEVs [32]. Briefly speaking, either the torque from the ICE or the torque on the outer rotor driven by AC current can be the prime power to drive the HEV, and they can also work together to enhance the output power. Additionally, the battery can be charged under a regenerative braking state. The concept that a PMSM is added after the CP-MGM is derived from that in E-CVT [33]. The PMSM in Figure 1 is used to regulate the performances of the outer rotor. For instance, it can be used to drag the outer 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]:

$$P_i+P_s=Q$$

(1)

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:

$$P_i{\omega }_i-Q{\omega }_o=P_s{w}_s=60f$$

(2)

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. Operating parameters of MGM under different modes of HEV.

Operation Modes	Rotating Speed of Inner Rotor ωo	Rotating Speed of Outer Rotor ωi	Current Frequency f
Pure electric mode	0 r/min	500 r/min	108.3 Hz
Pure mechanical mode	1200 r/min	1015 r/min	0 Hz (DC)
Hybrid mode	1200 r/min	2000 r/min	213.3 Hz
Regenerative braking mode	0 r/min	1000 r/min	216.6 Hz

. 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.

3. Mathematical Modeling of SCP-MGMs and DCP-MGMs

3.1. 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; 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:

$${\theta }_{PM}={\phi }_{in}+\frac{j\cdot 2\pi }{P_i}$$

(3)

$${\theta }_{Mod}=\frac{k\cdot 2\pi }{Q}+{\theta }_0-\frac{\beta }{2}$$

(4)

$${\theta }_{tooth}=t\cdot 2\pi /P$$

(5)

$${\theta }_{slot}=s\cdot 2\pi /P$$

(6)

where φ_in, and θ₀ 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, θ₀ 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; θ₀ 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 $\overrightarrow{A}$, the magnetic flux density $\overrightarrow{B}$, the magnetic field strength $\overrightarrow{H}$, and the current density distribution $\overrightarrow{J}$ in stator windings can be written in vector form as:

$$\overrightarrow{A}=A_z(r,\theta )\cdot {\overrightarrow{u}}_z$$

(7)

$$\overrightarrow{B}=B_r(r,\theta )\cdot {\overrightarrow{u}}_r+B_{\theta }(r,\theta )\cdot {\overrightarrow{u}}_{\theta }$$

(8)

$$\overrightarrow{H}=H_r(r,\theta )\cdot {\overrightarrow{u}}_r+H_{\theta }(r,\theta )\cdot {\overrightarrow{u}}_{\theta }$$

(9)

$$\overrightarrow{J}=J_z(r,\theta )\cdot {\overrightarrow{u}}_z$$

(10)

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:

$$A_z(r,\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{A}}_{z,n}(r)}\cdot e^{-in\theta }$$

(11)

$$B_r(r,\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{B}}_{r,n}(r)}\cdot e^{-in\theta }\mathrm{and}{B}_{\theta }(r,\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{B}}_{\theta ,n}(r)}\cdot e^{-in\theta }$$

(12)

$$H_r(r,\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{H}}_{r,n}(r)}\cdot e^{-in\theta }\mathrm{and}{H}_{\theta }(r,\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{H}}_{\theta ,n}(r)}\cdot e^{-in\theta }$$

(13)

$$J_z(r,\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{J}}_{z,n}(r)}\cdot e^{-in\theta }$$

(14)

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.

3.2. Partical Differential Equation Solution

The magnetic field within the machine follows quasistatic Maxwell equations:

$$\nabla \times \overrightarrow{H}=\overrightarrow{J}$$

(15)

$$\nabla \cdot \overrightarrow{B}=0$$

(16)

The relationship between $\overrightarrow{B}$ and $\overrightarrow{A}$ can be further expressed as:

$$\overrightarrow{B}=\nabla \times \overrightarrow{A}$$

(17)

The radial component and tangential component matrix of magnetic flux density $\overrightarrow{B}$ are then obtained in matrix form as [31]:

$$B_r=\frac{1}{r}\frac{\partial A_z}{\partial \theta }=-i\frac{1}{r}K{A}_z$$

(18)

$$B_{\theta }=-\frac{\partial A_z}{\partial r}$$

(19)

where K represents the harmonic order coefficient diagonal matrix that is related to N, given by:

$$K=\left[\begin{array}{ccc}-N& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & N\end{array}\right]$$

(20)

Similar to Equations (11)–(14), the relative permeability of each region can also be expressed in a complex Fourier series form:

$$\mu (\theta )={\displaystyle \sum _{n=-\infty }^{n=\infty }{\widehat{\mu }}_n}\cdot e^{-in\theta }$$

(21)

Next, based on the relation between $\overrightarrow{B}$ and $\overrightarrow{H}$, as expressed below:

$$\overrightarrow{B}=\mu \overrightarrow{H}+{\mu }_0\overrightarrow{M}$$

(22)

where $\overrightarrow{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:

$$B_r={\mathsf{\mu }}_{r,\mathrm{cov}}{H}_r+{\mu }_0{M}_r$$

(23)

$$B_{\theta }={\mathsf{\mu }}_{\theta ,\mathrm{cov}}{H}_{\theta }+{\mu }_0{M}_{\theta }$$

(24)

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:

$${\mathsf{\mu }}_{r,\mathrm{cov}}=\left[\begin{array}{ccc}{\widehat{\mu }}_0& \cdots & {\widehat{\mu }}_{-2N}\\ ⋮& \ddots & ⋮\\ {\widehat{\mu }}_{2N}& \cdots & {\widehat{\mu }}_0\end{array}\right]$$

(25)

where ${\widehat{\mu }}_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]:

$${\mathsf{\mu }}_{\theta ,\mathrm{cov}}={\left[\begin{array}{ccc}{\widehat{\mu }}_0^{rec}& \cdots & {\widehat{\mu }}_{-2N}^{rec}\\ ⋮& \ddots & ⋮\\ {\widehat{\mu }}_{2N}^{rec}& \cdots & {\widehat{\mu }}_0^{rec}\end{array}\right]}^{-1}$$

(26)

From Equation (26), it can be seen that μ_θ,cov is acquired by replacing ${\widehat{\mu }}_n$ with the corresponding n-th order Fourier coefficient of $1/{\mu }_{\theta }$ 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.

The relative permeability distribution on the circumferential direction can be expressed as:

$$\mu (\theta )=\{\begin{array}{ll}{\mu }_0\hfill & \theta \in [{\phi }_{out}+\frac{2(k-1)\pi }{Q},{\phi }_{out}+\frac{2(k-1)\pi }{Q}+\beta )\hfill \\ {\mu }_{iron,k}\hfill & \theta \in [{\phi }_{out}+\frac{2(k-1)\pi }{Q}+\beta ,{\phi }_{out}+\frac{2k\pi }{Q})\hfill \end{array}$$

(27)

where ${\phi }_{out}={\theta }_0-\frac{\beta }{2}$.

When a Fourier expansion on [0, 2π] is applied to Equation (27), the expressions of ${\widehat{\mu }}_k$ and ${\widehat{\mu }}_k^{rec}$ are given by:

$${\widehat{\mu }}_n=\{\begin{array}{l}{\displaystyle \sum _{k=1}^Q\frac{{\mu }_0}{2\pi in}{e}^{in(\frac{2k\pi }{Q}+{\theta }_0)}(e^{\frac{in\beta }{2}}-e^{-\frac{in\beta }{2}})}+{\displaystyle \sum _{k=1}^Q\frac{{\mu }_{iron,k}}{2\pi in}{e}^{in(\frac{2k\pi }{Q}+{\theta }_0)}(e^{in(\frac{2\pi }{Q}-\frac{\beta }{2})}-e^{\frac{in\beta }{2}})}n\ne 0\hfill \\ \frac{Q\beta {\mu }_0+(\frac{2\pi }{Q}-\beta ){\displaystyle \sum _{k=1}^Q{\mu }_{iron,k}}}{2\pi }n=0\hfill \end{array}$$

(28)

$${\widehat{\mu }}_n^{rec}=\{\begin{array}{l}{\displaystyle \sum _{k=1}^Q\frac{1}{2\pi in{\mu }_0}{e}^{in(\frac{2k\pi }{Q}+{\theta }_0)}(e^{\frac{in\beta }{2}}-e^{-\frac{in\beta }{2}})}+{\displaystyle \sum _{k=1}^Q\frac{1}{2\pi in{\mu }_{iron,k}}{e}^{in(\frac{2k\pi }{Q}+{\theta }_0)}(e^{in(\frac{2\pi }{Q}-\frac{\beta }{2})}-e^{\frac{in\beta }{2}})}n\ne 0\hfill \\ \frac{\frac{Q\beta }{{\mu }_0}+(\frac{2\pi }{Q}-\beta ){\displaystyle \sum _{k=1}^Q\frac{1}{{\mu }_{iron,k}}}}{2\pi }n=0\hfill \end{array}$$

(29)

The convolution matrices ${\mathsf{\mu }}_{r,\mathrm{cov}}$ and ${\mathsf{\mu }}_{\theta ,\mathrm{cov}}$ can then be obtained by substituting Equations (28) and (29) into (25) and (26).

Combining Equations (18)–(26) together, the VMP satisfies:

$$\frac{{\partial }^2{A}_z^k}{\partial r^2}+\frac{1}{r}\frac{\partial {{\rm A}}_z^k}{\partial r}-\frac{V{{\rm A}}_z^k}{r^2}=-{\mathsf{\mu }}_{\theta ,cov}{J}_z-\frac{{\mu }_0}{r}(M_{\theta }+iU{M}_r)$$

(30)

where $V={\mathsf{\mu }}_{\theta ,\mathrm{cov}}K{\mathsf{\mu }}_{r,\mathrm{cov}c}^{-1}K$ and $U={\mathsf{\mu }}_{\theta ,\mathrm{cov}}K{\mathsf{\mu }}_{r,\mathrm{cov}c}^{-1}$. The derivation process of Equation (30) is given in Appendix A.

Equation (30) is a Cauchy–Euler differential equation system [35]. The general solution of a single differential equation in Equation (30) is given by:

$$y=c_1{r}^{v^{\frac{1}{2}}}+c_2{r}^{-v^{\frac{1}{2}}}$$

(31)

where c₁, c₂ 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:

$$r^V(i,j)=r^{V(i,j)}$$

(32)

Hence, the complementary solution of Equation (30) ${{{\rm A}}_z^k|}_{com}$ can be written as:

$${{{\rm A}}_z^k|}_{com}=r^{V^{\frac{1}{2}}}{C}_1+r^{-V^{\frac{1}{2}}}{C}_2$$

(33)

where C₁ and C₂ are unknown coefficient vectors. Matrix $r^{V^{\frac{1}{2}}}$ can be factorized as:

$$r^{V^{\frac{1}{2}}}=P{r}^{\lambda }{P}^{-1}$$

(34)

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:

$${{{\rm A}}_z^k|}_{com}=P{r}^{\lambda }(P^{-1}\cdot C_1)+P{r}^{-\lambda }(P^{-1}\cdot C_2)=P{r}^{\lambda }D+P{r}^{-\lambda }E$$

(35)

As for the particular solution of Equation (30), ${{{\rm A}}_z^k|}_{par}$ is given by:

$${{{\rm A}}_z^k|}_{par}=r^{2}F+rG$$

(36)

where $F={(V-4I)}^{-1}{\mu }_{\theta ,cov}{J}_z$, $G={\mu }_0{(V-I)}^{-1}(M_{\theta }+iU{M}_r)$.

The general solution of VMP in Equation (30) is the sum of the complementary solution of Equation (35), and particular solution Equation (36):

$${{\rm A}}_z^k={{{\rm A}}_z^k|}_{com}+{{{\rm A}}_z^k|}_{par}$$

(37)

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, 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. Mathematical modeling of magnetic sources within SCP-MGMs and DCP-MGMs.

Sources	Illustrative Waveforms	Fourier Series Coefficients
Inner PM (Region III)		$M_{ri}(n)=\{\begin{array}{ll}\frac{\alpha P_i{B}_{r1}}{2\pi {\mu }_0}\hfill & \mathrm{n}=0\hfill \\ \frac{B_{r1}}{2n\pi {\mu }_{0}i}{e}^{in{P}_i{\phi }_i}\cdot \hfill & \\ (e^{in{P}_i\alpha }-1)& \mathrm{n}\ne 0\hfill \end{array}$
Outer PM (Only for DCP-MGM) (Region V)		$M_{ro}(n)=\{\begin{array}{ll}\frac{\beta Q{B}_{r2}}{2\pi {\mu }_0}\hfill & \mathrm{n}=0\hfill \\ \frac{B_{r1}}{2n\pi {\mu }_{0}i}{e}^{inQ{\theta }_0}\cdot \hfill & \\ (e^{\frac{inQ\beta }{2}}-e^{\frac{-inQ\beta }{2}})& \mathrm{n}\ne 0\end{array}$
Stator windings (Region VIII)		$\begin{array}{l}J(n)={\displaystyle \sum _{k=1}^P\frac{J(k)}{2n\pi i}{e}^{in\frac{2k\pi }{P_s}}}\cdot (e^{\frac{in\gamma }{2}}-e^{\frac{-in\gamma }{2}})\end{array}$

; their coefficients can be obtained by substituting a, b, and c in

Table 3. Mathematical modeling of the permeability distribution of different regions within SCP-MGMs and DCP-MGMs.

Regions	Illustrative Waveforms	Coefficients
III		a = φin, b = α, c = Pi
V		a = φin, b = β, c = Q
VII		a = δ/2, b = 2π/Ps − δ, c = Ps
VIII		a = γ/2, b = 2π/Ps − γ, c = Ps

into the following equations:

$${\widehat{\mu }}_n=\frac{1}{2\pi }({\displaystyle \sum _{k=1}^c{\displaystyle {\int }_a^{a+b}{\mu }_0{e}^{in\theta }d\theta }}+{\displaystyle \sum _{k=1}^c{\displaystyle {\int }_{a+b}^{a+2\pi /c}{\mu }_{iron}(k)e^{in\theta }d\theta }})$$

(38)

$${\widehat{\mu }}_n^{rec}=\frac{1}{2\pi }({\displaystyle \sum _{k=1}^c{\displaystyle {\int }_a^{a+b}\frac{e^{in\theta }}{{\mu }_0}d\theta }}+{\displaystyle \sum _{k=1}^c{\displaystyle {\int }_{a+b}^{a+2\pi /c}\frac{e^{in\theta }}{{\mu }_{iron,k}}d\theta }})$$

(39)

3.3. Bondary Condition Application

At the interfaces between two different subdomains, the radial component of $\overrightarrow{B}$ and the tangential component of $\overrightarrow{H}$ should be continuous across the boundary. Due to their ring-like shapes, each subdomain only interfaces with two subdomains at most, and region I and X have only one interface. Hence, all the boundary conditions can be written as:

$$B_r^k(R_k)=B_r^{k+1}(R_k)$$

(40)

$$H_{\theta }^k(R_k)=H_{\theta }^{k+1}(R_k)$$

(41)

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:

$$SX=T$$

(42)

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.

3.4. 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 μ₀ = 1500, is assigned to a given region with SMM for the first iteration. In the k-th iteration, the average flux density $\overrightarrow{B}$ of a specific region can be obtained by substituting ${\overline{\mu }}_{i,k}$ and solving the matrix Equation (42). A new average relative permeability ${\overline{\mu }}_{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:

$${\mu }_{i,k+1}=\frac{{\overline{\mu }}_{i,k}+{\overline{\mu }}_{i,kcal}}{2}$$

(43)

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:

$$\mathsf{\Delta }=\max \left\{\frac{\left|{\overline{\mu }}_{i,k}-{\overline{\mu }}_{i,cal}\right|}{{\overline{\mu }}_{i,k}}\right\},i\in \left\{III,V,VII,VIII\right\}$$

(44)

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. Relative permeabilities of different regions within SCP-MGMs and DCP-MGMs.

Regions	SCP-MGM	DCP-MGM
III	27.46	27.92
V	1309.15	1257.49
VII	888.05	838.35
VIII	1013.2	917.24

. 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.

3.5. 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]:

$${\phi }_k=\frac{N_{turn}L}{S_{coil}}{\displaystyle {\int }_{-\frac{\gamma }{2}+\frac{2k\pi }{P}}^{\frac{\gamma }{2}+\frac{2k\pi }{P}}{\displaystyle {\int }_{R_7}^{R_8}{A}_z^{VIII}(r,\theta )rdrd\theta }}$$

(45)

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:

$$\mathsf{\psi }=\left[\begin{array}{c}{\psi }_A\\ {\psi }_B\\ {\psi }_C\end{array}\right]=C_{turn}\left[\begin{array}{cccc}{\phi }_1& {\phi }_2& \cdots & {\phi }_{24}\end{array}\right]$$

(46)

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:

$$C_{turn}=\left[\begin{array}{cccccccccccccccccccccccc}1& 0& 0& 0& 0& -1& -1& 0& 0& 0& 0& 1& 1& 0& 0& 0& 0& -1& -1& 0& 0& 0& 0& 1\\ 0& 0& 0& 1& 1& 0& 0& 0& 0& -1& -1& 0& 0& 0& 0& 1& 1& 0& 0& 0& 0& -1& -1& 0\\ 0& -1& -1& 0& 0& 0& 0& 1& 1& 0& 0& 0& 0& -1& -1& 0& 0& 0& 0& 1& 1& 0& 0& 0\end{array}\right]$$

(47)

The back electromotive force (EMF) is computed by the derivative of ψ with respect to time:

$$E_{ABC}=\frac{d\mathsf{\psi }}{dt}=\frac{d\mathsf{\psi }}{d\theta }\cdot \omega$$

(48)

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:

$$T_{in}=\frac{L{R}_i^2}{{\mu }_0}{\displaystyle {\int }_0^{2\pi }{B}_r^{IV}(R_i,\theta )}{B}_{\theta }^{IV}(R_i,\theta )d\theta$$

(49)

$$T_{out}=\frac{L{R}_o^2}{{\mu }_0}{\displaystyle {\int }_0^{2\pi }{B}_r^{VI}(R_o,\theta )}{B}_{\theta }^{VI}(R_o,\theta )d\theta -\frac{L{R}_i^2}{{\mu }_0}{\displaystyle {\int }_0^{2\pi }{B}_r^{IV}(R_i,\theta )}{B}_{\theta }^{IV}(R_i,\theta )d\theta$$

(50)

where R_i and R_o are the middle radius of the inner air gap and outer air gap, respectively.

4. Validation and Comparison

4.1. 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. Geometrical parameters of the SCP-MGM and DCP-MGM.

Parameters	Symbols	Values	Units
Number of inner PM pole pairs	Pi	11	-
Number of modulator pieces	Q	13	-
Number of stator slots	P	24	-
Number of stator winding pole pairs	Ps	2	-
Axial length	L	110	mm
Outer radius of shaft	R1	25	mm
Inner radius of inner PM	R2	55.5	mm
Outer radius of inner rotor	R3	63.2	mm
Inner radius of outer rotor	R4	63.8	mm
Outer radius of outer rotor	R5	74.4	mm
Inner radius of stator	R6	75	mm
Radius of stator slot bottom	R7	78	mm
Outer radius of stator slot	R8	109	mm
Outer radius of stator	R9	120	mm
Angle of inner PM arc	α	0.286	rad
Angle of Modulator piece	β	0.242	rad
Angle of slot opening	δ	0.032	rad
Angle of stator slot	γ	0.168	rad
Slot filling factor	Fa	60%	-
Current density in stator windings	ID	5	A/mm2

. 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 ε₁ and maximum difference ε₂ of the two methods was defined as:

$${\epsilon }_1=\frac{1}{N_c}{\displaystyle \sum _{n=1}^{N_c}\left|\frac{V_{FEA,n}-V_{HMM,n}}{V_{FEA,n}}\right|}\times 100\%$$

(51)

$${\epsilon }_2=\max \left\{\left|V_{FEA,n}-V_{HMM,n}\right|\right\},\mathrm{n}\in \{1,2,\dots {,\mathrm{N}}_c\}$$

(52)

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.

4.2. Comparison between HMM and FEA

Based on the calculation and simulation, a good agreement between the HMM and FEA can be observed in Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 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 ε₁ and maximum difference ε₂ for the SCP-MGM and DCP-MGM are listed in

Table 6. Mean error percentage and maximum difference of the magnetic flux density distributions between the FEA and HMM of the SCP-MGM and DCP-MGM.

State	SCP-MGM				DCP-MGM
	Inner Rotor		Outer Rotor		Inner Rotor		Outer Rotor
	ε1	ε2	ε1	ε2	ε1	ε2	ε1	ε2
No-load	9.3%	0.55 T	15.9%	0.11 T	11.4%	1.06 T	14.4%	0.70 T
On-load	14.2%	1.62 T	14.3%	0.83 T	20.8%	3.11 T	16.7%	1.36 T

.

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.

4.3. 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.

4.3.1. 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.

4.3.2. Pure Electric Mode (Mode 1)

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 Figure 15 and Figure 16.

4.3.3. Pure Mechanical Mode (Mode 2)

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 Figure 17 and Figure 18.

4.3.4. Hybrid Mode (Mode 3)

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 Figure 19 and Figure 20.

4.3.5. 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 Figure 21 and Figure 22.

4.3.6. Quantitative comparison between HMM and FEA

Table 7. 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.

Mode	SCP-MGM				DCP-MGM
	Inner Rotor		Outer Rotor		Inner Rotor		Outer Rotor
	ε1	ε2	ε1	ε2	ε1	ε2	ε1	ε2
1	3.59%	12.1 Nm	2.64%	14.8 Nm	3.55%	22.5 Nm	4.46%	30.3 Nm
2	3.00%	7.85 Nm	7.39%	37.2 Nm	2.57%	14.9 Nm	5.88%	48.6 Nm
3	2.75%	8.9 Nm	3.61%	15.2 Nm	2.27%	15.3 Nm	3.27%	30.3 Nm
4	5.73%	18.1 Nm	5.06%	30.0 Nm	3.50%	28.4 Nm	4.56%	43.0 Nm

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.

4.3.7. 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.

5. 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.