*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:


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*<sup>z</sup> 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\_{\rm PM} = q\_{\rm in} + \frac{j \cdot 2\pi}{P\_i} \tag{3}
$$

$$
\theta\_{M\text{od}} = \frac{k \cdot 2\pi}{Q} + \theta\_0 - \frac{\beta}{2} \tag{4}
$$

$$
\theta\_{toch} = \mathbf{t} \cdot 2\pi/\mathbf{P} \tag{5}
$$

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

where ϕ*in*, and θ<sup>0</sup> 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,* θ<sup>0</sup> has a range of [0, 2π/*Pi*], [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; θ<sup>0</sup> 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.

**Figure 3.** Subdomain divisions of SCP-MGM and DCP-MGM.

The VMP <sup>→</sup> *<sup>A</sup>*, the magnetic flux density <sup>→</sup> *<sup>B</sup>*, the magnetic field strength <sup>→</sup> *H*, and the current density distribution <sup>→</sup> *J* in stator windings can be written in vector form as:

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

$$
\overrightarrow{B} = B\_r(r, \theta) \cdot \overrightarrow{\boldsymbol{\mu}}\_r + B\_\theta(r, \theta) \cdot \overrightarrow{\boldsymbol{\mu}}\_\theta \tag{8}
$$

$$
\overrightarrow{H} = H\_{\rm r}(r, \theta) \cdot \overrightarrow{\boldsymbol{u}}\_{\rm r} + H\_{\rm 0}(r, \theta) \cdot \overrightarrow{\boldsymbol{u}}\_{\rm 0} \tag{9}
$$

$$
\overrightarrow{J} = J\_z(r, \theta) \cdot \overrightarrow{u}\_z \tag{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) = \sum\_{n=-\infty}^{n=\infty} A\_{z,n}(r) \cdot e^{-in\theta} \tag{11}$$

$$B\_r(r, \theta) = \sum\_{n = -\infty}^{n = \infty} \mathcal{B}\_{r, n}(r) \cdot e^{-in\theta} \text{ and } B\_0(r, \theta) = \sum\_{n = -\infty}^{n = \infty} \mathcal{B}\_{0, n}(r) \cdot e^{-in\theta} \tag{12}$$

$$H\_{\mathbf{r}}(\mathbf{r},\theta) = \sum\_{n=-\infty}^{n=\infty} \hat{H}\_{\mathbf{r},\theta}(\mathbf{r}) \cdot \mathbf{e}^{-in\theta} \text{ and } H\boldsymbol{\varrho}(\mathbf{r},\theta) = \sum\_{n=-\infty}^{n=\infty} \hat{H}\boldsymbol{\varrho}\_{\mathbf{r}}(\mathbf{r}) \cdot \mathbf{e}^{-in\theta} \tag{13}$$

$$f\_z(r, \theta) = \sum\_{n = -\infty}^{n = \infty} f\_{z,n}(r) \cdot e^{-in\theta} \tag{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 Di*ff*erential Equation Solution*

The magnetic field within the machine follows quasistatic Maxwell equations:

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

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

The relationship between <sup>→</sup> *<sup>B</sup>* and <sup>→</sup> *A* can be further expressed as:

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

The radial component and tangential component matrix of magnetic flux density <sup>→</sup> *B* are then obtained in matrix form as [31]:

$$\mathbf{B}\_r = \frac{1}{r} \frac{\partial \mathbf{A}\_z}{\partial \theta} = -i \frac{1}{r} \mathbf{K} \mathbf{A}\_z \tag{18}$$

$$\mathbf{B}\_0 = -\frac{\partial \mathbf{A}\_z}{\partial r} \tag{19}$$

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

$$\mathbf{K} = \begin{bmatrix} -N & \cdots & 0 \\ \vdots & \ddots & \vdots \\ 0 & \cdots & N \end{bmatrix} \tag{20}$$

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

$$\mu(\theta) = \sum\_{n=-\infty}^{n=\infty} \hat{\mu}\_n \cdot e^{-in\theta} \tag{21}$$

Next, based on the relation between <sup>→</sup> *<sup>B</sup>* and <sup>→</sup> *H***,** as expressed below:

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

where <sup>→</sup> *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:

$$\mathbf{B}\_{\mathbf{r}} = \mu\_{\mathbf{r}, \text{cov}} \mathbf{H}\_{\mathbf{r}} + \mu\_0 \mathbf{M}\_{\mathbf{r}} \tag{23}$$

$$\mathbf{B}\_{\partial} = \mu\_{\partial, \text{cov}} \mathbf{H}\_{\partial} + \mu\_{\partial} \mathbf{M}\_{\partial} \tag{24}$$

where μ*r*,cov and μθ,cov are convolution matrices of the radial and tangential components of permeability, respectively. **M***<sup>r</sup>* and **M**<sup>θ</sup> are the radial and tangential components of magnetization intensity, respectively. **M***<sup>r</sup>* and **M**<sup>θ</sup> can all written in complex Fourier series. The convolution matrix μ*r*,cov can be defined as:

$$
\boldsymbol{\mu}\_{r, \text{cov}} = \begin{bmatrix}
\boldsymbol{\mu}\_0 & \cdots & \boldsymbol{\mu}\_{-2N} \\
\vdots & \ddots & \vdots \\
\boldsymbol{\mu}\_{2N} & \cdots & \boldsymbol{\mu}\_0 \\
\end{bmatrix} \tag{25}
$$

where μˆ *<sup>n</sup>* is the *n*-th order coefficient of the Fourier series of corresponding μ.

In Equation (24), **H**<sup>θ</sup> is continuous at the interface between two regions, but B<sup>θ</sup> 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]:

$$\boldsymbol{\mu}\_{0,\text{cov}} = \begin{bmatrix} \hat{\mu}\_0^{\text{rec}} & \cdots & \hat{\mu}\_{-2N}^{\text{rec}} \\ \vdots & \ddots & \vdots \\ \hat{\mu}\_{2N}^{\text{rec}} & \cdots & \hat{\mu}\_0^{\text{rec}} \end{bmatrix}^{-1} \tag{26}$$

From Equation (26), it can be seen that μθ,cov is acquired by replacing μˆ *<sup>n</sup>* 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.

**Figure 4.** The calculation instance of the convolution matrix with respect to relative permeability.

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

$$\mu(\theta) = \begin{cases} \mu\_0 & \theta \in [\varrho\_{out} + \frac{2(k-1)\pi}{Q}, \varrho\_{out} + \frac{2(k-1)\pi}{Q} + \beta) \\\ \mu\_{iron,k} & \theta \in [\varrho\_{out} + \frac{2(k-1)\pi}{Q} + \beta, \varrho\_{out} + \frac{2k\pi}{Q}) \end{cases} \tag{27}$$

where <sup>ϕ</sup>*out* <sup>=</sup> <sup>θ</sup><sup>0</sup> <sup>−</sup> <sup>β</sup> 2 .

When a Fourier expansion on [0, 2π] is applied to Equation (27), the expressions of μˆ *<sup>k</sup>* and μˆ*rec <sup>k</sup>* are given by:

$$\mu\_n = \begin{cases} \sum\_{k=1}^{\mathcal{Q}} \frac{\mu\_0}{2\pi i n} e^{in(\frac{2\pi}{\mathcal{Q}} + \mathcal{O}\_0)} (e^{\frac{in\mathcal{C}}{2}} - e^{-\frac{in\mathcal{C}}{2}}) + \sum\_{k=1}^{\mathcal{Q}} \frac{\mu\_{\text{inv}k}}{2\pi i n} e^{in(\frac{2\pi}{\mathcal{Q}} + \mathcal{O}\_0)} (e^{in(\frac{2\pi}{\mathcal{Q}} - \frac{\mathcal{C}}{2})} - e^{\frac{in\mathcal{C}}{2}}) n \neq 0 \\\ \underbrace{\mathcal{Q}\mu\_{00} + (\frac{2\pi}{\mathcal{Q}} - \mathcal{J}) \sum\_{k=1}^{\mathcal{Q}} \mu\_{\text{inv}k}}\_{2\pi} n = 0 \end{cases} \tag{28}$$

$$\begin{split} \mu\_{n}^{rec} &= \begin{cases} \sum\_{k=1}^{\mathbb{Q}} \frac{1}{2\pi i m \mu\_{0}} e^{j n \left(\frac{2k\pi}{\mathbb{Q}} + \theta\_{0}\right)} \left(e^{\frac{i\pi\theta}{2}} - e^{-\frac{i\pi\theta}{2}}\right) + \sum\_{k=1}^{\mathbb{Q}} \frac{1}{2\pi i m \mu\_{imk}} e^{j n \left(\frac{2k\pi}{\mathbb{Q}} + \theta\_{0}\right)} \left(e^{j n \left(\frac{2\pi}{\mathbb{Q}} - \frac{\theta}{2}\right)} - e^{\frac{i\pi\theta}{2}}\right) n \neq 0\\ \frac{\frac{\sqrt{\mu}}{\mathbb{Q}} + \left(\frac{2\pi}{\mathbb{Q}} - \theta\right)}{2\pi} \frac{1}{\frac{1}{2\pi i m \mu\_{k}}} n = 0 \end{cases} \tag{29} \\ \end{split} \tag{20}$$

The convolution matrices μ*r*,cov and μθ,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 \mathbf{A}\_z^k}{\partial r^2} + \frac{1}{r} \frac{\partial \mathbf{A}\_z^k}{\partial r} - \frac{\mathbf{V} \mathbf{A}\_z^k}{r^2} = -\mu\_{\partial, \text{av}} \mathbf{J}\_z - \frac{\mu\_0}{r} (\mathbf{M}\_\partial + i \mathbf{U} \mathbf{M}\_r) \tag{30}$$

where **V** = μθ,cov**K**μ−<sup>1</sup> *<sup>r</sup>*,cov*c***K** and **U** = μθ,cov**K**μ−<sup>1</sup> *<sup>r</sup>*,cov*c*. 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^{\nu \frac{1}{2}} + c\_2 r^{-\nu \frac{1}{2}} \tag{31}$$

where *c*1, *c*<sup>2</sup> 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***<sup>V</sup>** is defined as:

$$r^{\mathbf{V}}(i,j) = r^{\mathbf{V}(i,j)} \tag{32}$$

Hence, the complementary solution of Equation (30) **A***<sup>k</sup> z com* can be written as:

$$\left. \mathbf{A}\_z^k \right|\_{\rm com} = r^{\mathbf{V}^\frac{\tilde{l}}{2}} \mathbf{C}\_1 + r^{-\mathbf{V}^\frac{\tilde{l}}{2}} \mathbf{C}\_2 \tag{33}$$

where **<sup>C</sup>**<sup>1</sup> and **<sup>C</sup>**<sup>2</sup> are unknown coefficient vectors. Matrix *<sup>r</sup>***<sup>V</sup>** <sup>1</sup> <sup>2</sup> can be factorized as:

$$r^{\mathbf{V}^{\frac{1}{2}}} = \mathbf{P}r^{\lambda}\mathbf{P}^{-1} \tag{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:

$$\left. \mathbf{A}\_z^k \right|\_{\rm com} = \mathbf{Pr}^\lambda(\mathbf{P}^{-1} \cdot \mathbf{C}\_1) + \mathbf{Pr}^{-\lambda}(\mathbf{P}^{-1} \cdot \mathbf{C}\_2) = \mathbf{Pr}^\lambda \mathbf{D} + \mathbf{Pr}^{-\lambda} \mathbf{E} \tag{35}$$

As for the particular solution of Equation (30), **A***<sup>k</sup> z par* is given by:

$$\mathbf{A}\_z^k \Big|\_{par} = r^2 \mathbf{F} + r \mathbf{G} \tag{36}$$

where **F** = (**V** − 4**I**) −1 μθ,*cov***J***z*, **G** = μ0(**V** − **I**) −1 (**M**<sup>θ</sup> + *i***UM***r*).

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

$$\mathbf{A}\_z^k = \mathbf{A}\_z^k \big|\_{com} + \mathbf{A}\_z^k \big|\_{par} \tag{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; their coefficients can be obtained by substituting *a*, *b*, and *c* in Table 3 into the following equations:

$$\hat{\mu}\_{\text{lt}} = \frac{1}{2\pi} (\sum\_{k=1}^{c} \int\_{a}^{a+b} \mu\_0 \varepsilon^{in0} d\theta + \sum\_{k=1}^{c} \int\_{a+b}^{a+2\pi/c} \mu\_{\text{iron}}(k) \varepsilon^{in0} d\theta) \tag{38}$$

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

**Table 2.** Mathematical modeling of magnetic sources within SCP-MGMs and DCP-MGMs.

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


#### *3.3. Bondary Condition Application*

At the interfaces between two different subdomains, the radial component of <sup>→</sup> *B* and the tangential component of <sup>→</sup> *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:

$$\mathbf{B}\_r^k(\mathcal{R}\_k) = \mathbf{B}\_r^{k+1}(\mathcal{R}\_k) \tag{40}$$

$$\mathbf{H}\_{\mathcal{O}}^{k}(\mathcal{R}\_{k}) = \mathbf{H}\_{\mathcal{O}}^{k+1}(\mathcal{R}\_{k}) \tag{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:

$$\mathbf{S}\mathbf{X} = \mathbf{T} \tag{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 μ<sup>0</sup> = 1500, is assigned to a given region with SMM for the first iteration. In the *<sup>k</sup>*-th iteration, the average flux density <sup>→</sup> *B* of a specific region can be obtained by substituting μ*i*,*<sup>k</sup>* 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*+<sup>1</sup> for the (*k* + 1)-th iteration in region *i* is given as:

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

The iteration will stop only if the iteration time *ni* exceeds the maximum number of iterations *Ni* (*Ni* is set to be 50 here), or the maximum error Δ in all these regions is below the error requirement ξ, Δ is defined as:

$$\Delta = \max \left\{ \frac{\left| \overline{\mu}\_{i,k} - \overline{\mu}\_{i,\text{cal}} \right|}{\overline{\mu}\_{i,k}} \right\}\_{\prime} \text{ } i \in \{III\_{\prime} \\ V\_{\prime} \, VII\_{\prime} \, VIII \} \tag{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. 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.

**Figure 6.** Fast Fourier transform (FFT) of no-load outer-air-gap radial magnetic flux density for SCP-MGMs and DCP-MGMs.

**Table 4.** Relative permeabilities of different regions within SCP-MGMs and DCP-MGMs.


#### *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]:

$$q\rho\_k = \frac{N\_{\rm turn}L}{S\_{\rm coil}} \int\_{-\frac{\gamma}{2} + \frac{2k\pi}{P}}^{\frac{\gamma}{2} + \frac{2k\pi}{P}} \int\_{R\_7}^{R\_8} A\_z^{VIII}(r, \theta) r dr d\theta \tag{45}$$

where *Nturn* is the number of coil turns in each slot, and *Scoil* 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:

$$\boldsymbol{\upmu} = \begin{bmatrix} \psi\_A \\ \psi\_B \\ \psi\_C \end{bmatrix} = \mathbf{C}\_{\text{turn}} \begin{bmatrix} \boldsymbol{\uprho}\_1 & \boldsymbol{\uprho}\_2 & \cdots & \boldsymbol{\uprho}\_{24} \end{bmatrix} \tag{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* = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ 1 0 0 00 −1 −100 0 0 11 0 0 00 −1 −100 0 0 1 0 0 0 11 0 0 00 −1 −100 0 0 11 0 0 00 −1 −1 0 0 −1 −100 0 0 11 0 0 00 −1 −100 0 0 11 0 0 0 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (47)

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

$$\mathbf{E}\_{\rm ABC} = \frac{d\boldsymbol{\Psi}}{dt} = \frac{d\boldsymbol{\Psi}}{d\boldsymbol{\Theta}} \cdot \boldsymbol{\omega} \tag{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 *Tin* 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 *Tout* 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\_{\dot{m}} = \frac{LR\_i^2}{\mu\_0} \int\_0^{2\pi} B\_r^{IV}(\mathcal{R}\_{i\dot{\nu}}\theta) B\_0^{IV}(\mathcal{R}\_{i\dot{\nu}}\theta) d\theta \tag{49}$$

$$T\_{\rm out} = \frac{LR\_o^2}{\mu\_0} \int\_0^{2\pi} B\_r^{VI}(\mathbf{R}\_{o\prime}, \theta) B\_0^{VI}(\mathbf{R}\_{o\prime}, \theta) d\theta - \frac{LR\_i^2}{\mu\_0} \int\_0^{2\pi} B\_r^{IV}(\mathbf{R}\_{i\prime}, \theta) B\_0^{IV}(\mathbf{R}\_{i\prime}, \theta) d\theta \tag{50}$$

where *Ri* and *Ro* 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. 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 ε<sup>1</sup> and maximum difference ε<sup>2</sup> of the two methods was defined as:

$$\varepsilon\_1 = \frac{1}{N\_c} \sum\_{n=1}^{N\_c} \left| \frac{V\_{FEA,n} - V\_{HMM,n}}{V\_{FEA,n}} \right| \times 100\% \tag{51}$$

$$\varepsilon\_2 = \max \{ \left| V\_{FEA, \mathbf{u}} - V\_{HMM, \mathbf{u}} \right| \}, \mathbf{n} \in \{ 1, 2, \dots, N\_{\mathbf{c}} \} \tag{52}$$

where *VHMM,n* is the *n*-th value, calculated using the HMM, and *VFEA*,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 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 ε<sup>1</sup> and maximum difference ε<sup>2</sup> for the SCP-MGM and DCP-MGM are listed in Table 6.

**Figure 7.** On-load inner and outer air-gap radial magnetic flux density distribution of the SCP-MGM: (**a**) outer air gap; (**b**) inner air gap.

**Figure 8.** On-load inner and outer air-gap radial magnetic flux density distribution of the SCP-MGM: (**a**) outer air gap; (**b**) inner air gap.

**Figure 9.** Comparison of on-load magnetic flux density distribution of the SCP-MGM, drawn by harmonic modeling method (HMM) and finite element analysis (FEA): (**a**) HMM; (**b**) FEA.

**Figure 10.** No-load inner and outer air-gap radial magnetic flux density distribution of the DCP-MGM: (**a**) outer air gap; (**b**) inner air gap.

**Figure 11.** On-load inner and outer air-gap radial magnetic flux density distribution of the DCP-MGM: (**a**) outer air gap; (**b**) inner air gap.

**Figure 12.** Comparison of on-load magnetic flux density distribution of the DCP-MGM drawn by HMM and FEA: (**a**) HMM; (**b**) FEA.


**Table 5.** Geometrical parameters of the SCP-MGM and DCP-MGM.


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

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.

**Figure 13.** FFT of no-load outer-air-gap radial magnetic flux density for the SCP-MGM and the DCP-MGM: (**a**) SCP-MGM; (**b**) DCP-MGM.

#### *4.3. Electromagnetic Performance Analysis under Di*ff*erent 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 ω<sup>i</sup> = 1200 r/min, ω<sup>o</sup> = 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.

**Figure 14.** Analytically predicted and FEA simulated back electromotive force (EMF) of the SCP-MGM and DCP-MGM: (**a**) SCP-MGM; (**b**) DCP-MGM.

#### 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, ω<sup>i</sup> = 0 r/min, ω<sup>o</sup> = 500 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 15 and 16.

**Figure 15.** Analytically predicted and FEA simulated torque waveforms of the SCP-MGM and DCP-MGM under pure electric mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

**Figure 16.** Power distribution among the internal combustion engine (ICE), battery, and wheels under pure electric mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

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: ω<sup>i</sup> = 1200 r/min, ω<sup>o</sup> = 1015 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 17 and 18.

**Figure 17.** Analytically predicted and FEA simulated torque waveforms of SCP-MGM and DCP-MGM under pure mechanical mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

**Figure 18.** Power distribution among ICE, battery and wheels under pure mechanical mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

### 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: ω<sup>i</sup> = 1200 r/min, ω<sup>o</sup> = 2000 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 19 and 20.

**Figure 19.** Analytically predicted and FEA simulated electromagnetic torque of SCP-MGM and DCP-MGM under hybrid mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

**Figure 20.** Power distribution among ICE, battery, and wheels under hybrid mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

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: ω<sup>i</sup> = 0 r/min, ω<sup>o</sup> = 1000 r/min. The torque waveforms of the power distribution of the two machines are shown in Figures 21 and 22.

**Figure 21.** Analytically predicted and FEA simulated back EMF of SCP-MGM and DCP-MGM: (**a**) SCP-MGM; (**b**) DCP-MGM.

**Figure 22.** Power distribution among ICE, battery, and wheels under hybrid mode: (**a**) SCP-MGM; (**b**) DCP-MGM.

#### 4.3.6. Quantitative comparison between HMM and FEA

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.

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


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.

**Author Contributions:** The work presented in this paper is the output of the research projects undertaken by C.L. In specific, H.Z. and C.L. developed the topic. H.Z. carried out the calculation and simulation, analyzed the results, and wrote the paper. Z.S. gave some suggestions on the calculation process. J.Y. helped to carried out a part of the simulation.

**Funding:** This research was supported by a grant (Project No. 51677159) from the Natural Science Foundation of China (NSFC), China. Also, the work is supported by a grant (Project No. JCYJ20180307123918658) from the Science Technology and Innovation Committee of Shenzhen Municipality, Shenzhen, China. Moreover, it was also supported by a grant (Project No. CityU 21201216) from Research Grants Council of HKSAR, China.

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

#### **Appendix A**

The expressions of **H***<sup>r</sup>* and **H**<sup>θ</sup> are given by:

$$\mathbf{H}\_{\rm r} = -i\frac{1}{r}\mu\_{r,\rm cov}^{-1}\mathbf{K}\mathbf{A}\_{\rm z} - \mu\_{0}\mu\_{r,\rm cov}^{-1}\mathbf{M}\_{\rm r} \text{ and } \mathbf{H}\_{\rm 0} = -\mu\_{0,\rm cov}^{-1}\frac{\partial\mathbf{A}\_{\rm z}}{\partial r} - \mu\_{0}\mu\_{0,\rm cov}^{-1}\mathbf{M}\_{\rm 0} \tag{A1}$$

In polar coordinate, (15) can be simplified as:

$$\nabla \times \mathbf{H} = \frac{1}{r} \frac{\partial}{\partial r} (r \mathbf{H}\_{\partial}) + i \frac{\mathbf{K} \mathbf{H}\_{r}}{r} = \mathbf{J}\_{z} \tag{A2}$$

By substituting (A1) into (A2), one can obtain that:

$$\frac{\partial^2 \mathbf{A}\_z^k}{\partial r^2} + \frac{1}{r} \frac{\partial \mathbf{A}\_z^k}{\partial r} - \frac{b\_{\partial, \text{cov}} \mathbf{K} \mu\_{r, \text{cov}}^{-1} \mathbf{K} \mathbf{A}\_z^k}{r^2} = -\mu\_{\partial, \text{cov}} \mathbf{J}\_z - \frac{\mu\_0}{r} (\mathbf{M}\_\partial + i \mu\_{\partial, \text{cov}} \mathbf{K} \mu\_{r, \text{cov}}^{-1} \mathbf{M}\_r) \tag{A3}$$

The vector magnetic potential in region I~X are given in Table A1.


**Table A1.** The Vector Magnetic Potential Expressions.

Where **D**<sup>I</sup> , **D**II, ... , **D**IX and **E**<sup>I</sup> , **E**II, ... , **E**<sup>X</sup> are (N\*1) vector.

The expressions of **X**, **S**, **T** are:

**<sup>X</sup>** = (**D***<sup>I</sup>* ) *<sup>T</sup>* (**D***II*) *<sup>T</sup>* (**E***II*) *<sup>T</sup>* (**D***III*) *<sup>T</sup>* (**E***III*) *<sup>T</sup>* ··· (**D***IX*) *<sup>T</sup>* (**E***IX*) *<sup>T</sup>* (**E***X*) *<sup>T</sup> <sup>T</sup>* (A4) **S** = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ **K**1,1 **K**1,2 **K**1,3 0000 0 0 0 0 0 0 0 0 0 0 0 **K**2,1 **K**2,2 **K**2,3 0000 0 0 0 0 0 0 0 0 0 0 0 0 **K**3,2 **K**3,3 **K**3,4 **K**3,5 00 0 0 0 0 0 0 0 0 0 0 0 0 **K**4,2 **K**4,3 **K**4,4 **K**4,5 00 0 0 0 0 0 0 0 0 0 0 0 000 **K**5,4 **K**5,5 **K**5,6 **K**5,7 00 0 0 0 0 0 0 0 0 0 000 **K**6,4 **K**6,5 **K**6,6 **K**6,7 00 0 0 0 0 0 0 0 0 0 00000 **K**7,6 **K**7,7 **K**7,8 **K**7,9 000000000 00000 **K**8,6 **K**8,7 **K**8,8 **K**8,9 000000000 0000000 **K**9,8 **K**9,9 **K**9,10 **K**9,11 0000000 0000000 **K**10,8 **K**10,9 **K**10,10 **K**10,11 0000000 0000000 0 0 **K**11,10 **K**11,11 **K**11,12 **K**11,13 00000 0000000 0 0 **K**12,10 **K**12,11 **K**12,12 **K**12,13 00000 0000000 0 0 0 0 **K**13,12 **K**13,13 **K**13,14 **K**13,15 000 0000000 0 0 0 0 **K**14,12 **K**14,13 **K**14,14 **K**14,15 000 0000000 0 0 0 0 0 0 **K**15,14 **K**15,15 **K**15,16 **K**15,17 0 0000000 0 0 0 0 0 0 **K**16,14 **K**16,15 **K**16,16 **K**16,17 0 0000000 0 0 0 0 0000 **K**17,16 **K**17,17 **K**17,18 0000000 0 0 0 0 0000 **K**18,16 **K**18,17 **K**18,18 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ (A5) **<sup>T</sup>** = 0 0 (*R*2**G**1) *<sup>T</sup>* (*R*2μ*II***G**1) *<sup>T</sup>* (−*R*3**G**1) *<sup>T</sup>* (−*R*3μ*IV***G**1) *<sup>T</sup>* (*R*4**G**2) *<sup>T</sup>* (*R*4μ*IV***G**2) *<sup>T</sup>* (−*R*5**G**2) *T <sup>T</sup>* 0 0 (*R*<sup>2</sup> *<sup>T</sup>* (2*R*<sup>2</sup> *<sup>T</sup>* (−*R*<sup>2</sup> *<sup>T</sup>* (−2*R*<sup>2</sup> *<sup>T</sup>* 0 0 *<sup>T</sup>* (A6)

$$\text{where:}$$

(−μ*VIR*5**G**2)

$$\begin{bmatrix} \begin{array}{cc} \mathbf{K}\_{1,1} & \mathbf{K}\_{1,2} & \mathbf{K}\_{1,3} \end{array} \end{bmatrix} = \begin{bmatrix} \begin{array}{cc} \mathbf{I} & -\begin{pmatrix} \frac{\mathbf{R}\_1}{\|\mathbf{R}\_2\|} & -\mathbf{I} \end{pmatrix} \end{bmatrix} \end{bmatrix} \tag{A7}$$

<sup>8</sup>μ*IX***F**)

<sup>8</sup>**F**)

$$\begin{bmatrix} \begin{array}{cc} \mathbf{K}\_{2,1} & \mathbf{K}\_{2,2} & \mathbf{K}\_{2,3} \end{array} \end{bmatrix} = \begin{bmatrix} \begin{array}{cc} \boldsymbol{\mu}\_{II} & -\boldsymbol{\mu}\_{I} \begin{pmatrix} \frac{\boldsymbol{R}\_{1}}{\boldsymbol{R}\_{2}} \end{pmatrix}^{\|\mathbf{K}\|} & \boldsymbol{\mu}\_{I} \end{bmatrix} \end{bmatrix} \tag{A8}$$

<sup>7</sup>μ*open*,θ**F**)

<sup>7</sup>**F**)

$$\begin{bmatrix} \mathbf{K}\_{3,2} & \mathbf{K}\_{3,3} & \mathbf{K}\_{3,4} & \mathbf{K}\_{3,5} \end{bmatrix} = \begin{bmatrix} \mathbf{I} & \begin{pmatrix} \frac{\mathbf{R}\_1}{\mathbf{R}\_2} \end{pmatrix}^{|\mathbf{K}|} & -\mathbf{P}^{III} \begin{pmatrix} \frac{\mathbf{R}\_1}{\mathbf{R}\_2} \end{pmatrix}^{|\mathbf{K}|} & -\mathbf{P}^{III} \end{bmatrix} \tag{A9}$$

$$\begin{bmatrix} \mathbf{K}\_{4,2} & \mathbf{K}\_{4,3} & \mathbf{K}\_{4,4} & \mathbf{K}\_{4,5} \end{bmatrix} = \begin{bmatrix} \begin{array}{cc} \mu\_{\rm PM,0} & -\mu\_{\rm PM,0} \lVert \mathbf{K} \rVert \begin{pmatrix} \frac{R\_1}{R\_2} \end{pmatrix}^{\lVert \mathbf{K} \rVert} & -\mu\_{\rm II} \mathbf{P}^{\textup{III}} \lambda^{\textup{III}} \begin{pmatrix} \frac{R\_2}{R\_3} \end{pmatrix}^{\lambda^{\textup{III}}} & \mu\_{\rm II} \mathbf{P}^{\textup{III}} \lambda^{\textup{III}} \end{bmatrix} & \begin{array}{cc} \text{(A10)} \end{array}$$

$$
\begin{bmatrix}
\begin{array}{cccc}
\mathbf{K}\_{5,4} & \mathbf{K}\_{5,5} & \mathbf{K}\_{5,6} & \mathbf{K}\_{5,7}
\end{array}
\end{bmatrix} = \begin{bmatrix}
\begin{array}{cccc}
\mathbf{p}^{\mathrm{III}} & \mathbf{p}^{\mathrm{III}}
\end{array}
\begin{array}{cccc}
\mathbf{p}^{\mathrm{III}} & \begin{array}{cccc}
\mathbf{p}^{\mathrm{III}} & \begin{array}{c}
\mathbf{p}^{\mathrm{III}}
\end{array}
\end{array}
\begin{array}{cccc}
\mathbf{k}^{\mathrm{III}} & \begin{array}{c}
\mathbf{k}^{\mathrm{I}\mathbb{I}}
\end{array}
\end{array}
\begin{array}{cccc}
\mathbf{k}^{\mathrm{I}\mathbb{I}}
\end{array}
\end{bmatrix}
\end{bmatrix}
\tag{A11}
$$

$$
\begin{bmatrix}
\mathbf{K}\_{6,4} & \mathbf{K}\_{6,5} & \mathbf{K}\_{6,6} & \mathbf{K}\_{6,7}
\end{bmatrix} = \begin{bmatrix}
\begin{array}{ccccc}
\boldsymbol{\mu}\_{\rm IV} \mathbf{P}^{\rm II} \boldsymbol{\lambda}^{\rm III} & -\boldsymbol{\mu}\_{\rm IV} \mathbf{P}^{\rm III} \boldsymbol{\lambda}^{\rm III}
\end{array}
\begin{array}{ccccc}
\boldsymbol{\mu}\_{\rm PL} |\mathbf{K}| \begin{pmatrix}
\boldsymbol{R}\_{3} \\
\boldsymbol{R}\_{3}
\end{pmatrix}^{\rm IV} & \boldsymbol{\mu}\_{\rm PL} |\mathbf{K}| \begin{pmatrix}
\boldsymbol{R}\_{3} \\
\boldsymbol{R}\_{3}
\end{pmatrix}^{\rm IV} & \boldsymbol{\mu}\_{\rm PL} |\mathbf{K}|
\end{bmatrix}
\end{bmatrix} \tag{A12}
$$

$$
\begin{bmatrix}
\boldsymbol{\nu} & \boldsymbol{\nu} & \boldsymbol{\nu} & \boldsymbol{\nu} & \boldsymbol{\nu}
\end{bmatrix} \begin{bmatrix}
\boldsymbol{\nu}\_{\rm L} (|\mathbf{K}| & \boldsymbol{R}\_{\rm PL}) \boldsymbol{\lambda}^{\rm IV} & \boldsymbol{\nu}\_{\rm PL}
\end{bmatrix}
\tag{A13}
$$

 **<sup>K</sup>**7,6 **<sup>K</sup>**7,7 **<sup>K</sup>**7,8 **<sup>K</sup>**7,9 = **<sup>I</sup>** *R*<sup>3</sup> *R*4 |**K**<sup>|</sup> <sup>−</sup>**P***VR*<sup>4</sup> *R*5 λ*<sup>V</sup>* <sup>−</sup>**P***<sup>V</sup>* (A13)

$$
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{8,6} & \mathbf{K}\_{8,7} & \mathbf{K}\_{8,8} & \mathbf{K}\_{8,9} \\
\end{bmatrix}
\end{bmatrix} = \begin{bmatrix}
\begin{bmatrix}
\mu\_{\mathrm{Mod},\beta} \mathrm{|}\mathbf{K}\rangle & -\mu\_{\mathrm{Mod},\beta} \mathrm{|}\mathbf{K}\big{\big{(}}^{\frac{\mathrm{R}}{\mathrm{R}\_{4}}\big{)}\mathbf{K}\big{\big{(}}^{\frac{\mathrm{R}}{\mathrm{R}\_{5}}\big{)}} -\mu\_{\mathrm{IV}} \mathrm{P}\big{\mathrm{V}}\big{\mathrm{V}}^{\mathrm{V}}\_{\mathrm{R}\_{5}}\big{)}^{\mathrm{V}} & \mu\_{\mathrm{IV}} \mathrm{P}\big{\mathrm{V}}^{\mathrm{V}}
\end{bmatrix}
\end{bmatrix}
\tag{A14}
$$

$$
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{0,8} & \mathbf{K}\_{0,0} & \mathbf{K}\_{0,0} & \mathbf{K}\_{0,0}
\end{bmatrix}
\mathbf{n}^{\mathrm{V}} & \mathbf{n}^{\mathrm{V}}\_{1} \mathrm{(}\boldsymbol{\mathbb{R}}^{\mathrm{V}}\_{4}\big{)}^{\mathrm{V}} & \left(\boldsymbol{\mathbb{R}}\_{8}\right)^{\mathrm{[}}\mathrm{K}\big{\big{}}^{\mathrm{R}}
\end{bmatrix}
\end{bmatrix}
\tag{A15}
$$

$$
\begin{bmatrix}
\mathbf{K}\_{9,8} & \mathbf{K}\_{9,9} & \mathbf{K}\_{9,10} & \mathbf{K}\_{9,11} \\
\end{bmatrix} = \begin{bmatrix}
\mathbf{p}^V & \mathbf{p}^V \begin{pmatrix}
\frac{R\_4}{R\_5} \end{pmatrix}^V & -\begin{pmatrix}
\frac{R\_5}{R\_6} \end{pmatrix}^{|\mathbf{K}|} & -\mathbf{I} \\
& \dots & \\
\end{bmatrix} \tag{A15}
$$

$$
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{10,8} & \mathbf{K}\_{10,9} & \mathbf{K}\_{10,10} & \mathbf{K}\_{10,11} \\
\end{bmatrix}
\end{bmatrix} = \begin{bmatrix}
\begin{bmatrix}
\boldsymbol{\mu}\_{\mathrm{VI}}\mathbf{P}^{\mathrm{V}}\boldsymbol{\lambda}^{\mathrm{V}} & -\boldsymbol{\mu}\_{\mathrm{VI}}\mathbf{P}^{\mathrm{V}}\boldsymbol{\lambda}^{\mathrm{V}}
\end{bmatrix}
\begin{bmatrix}
\frac{\boldsymbol{\mathcal{R}}\_{\mathrm{S}}}{\boldsymbol{\mathcal{R}}\_{\mathrm{S}}}
\end{bmatrix}
\begin{bmatrix}
\frac{\boldsymbol{\mathcal{R}}\_{\mathrm{S}}}{\boldsymbol{\mathcal{R}}\_{\mathrm{S}}}
\end{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{1}
\end{bmatrix}
\end{bmatrix}
\begin{bmatrix}
\mathbf{K}
\end{bmatrix}
\end{bmatrix}
\tag{A16}
$$

$$
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{1}\boldsymbol{\omega}\dots\mathbf{K}\_{10} & \mathbf{K}\_{10,10} & \mathbf{K}\_{10,10}
\end{bmatrix}
\end{bmatrix} = \begin{bmatrix}
\begin{bmatrix}
\mathbf{I}\_{1}
\end{bmatrix}
\begin{bmatrix}
\mathbf{K}
\end{bmatrix}
\end{bmatrix}
\begin{bmatrix}
\mathbf{K}
\end{bmatrix}
\mathbf{J}^{\mathrm{V}\mathrm{I}}
\end{bmatrix}
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}
\end{bmatrix}
\end{bmatrix}
\tag{A17}
$$

$$\begin{bmatrix} \mathbf{K}\_{11,10} & \mathbf{K}\_{11,11} & \mathbf{K}\_{11,12} & \mathbf{K}\_{11,13} \end{bmatrix} = \begin{bmatrix} \mathbf{I} & \begin{pmatrix} \frac{R\varsigma}{R} \end{pmatrix}^{\vert \mathbf{K} \vert} & -\mathbf{P}^{VII} \begin{pmatrix} \frac{R\varsigma}{R} \end{pmatrix}^{\lambda^{\prime III}} & -\mathbf{P}^{VIII} \end{bmatrix} \tag{A17}$$

$$\begin{bmatrix} \begin{array}{cccc} \mathbf{K}\_{12,10} & \mathbf{K}\_{12,11} & \mathbf{K}\_{12,12} & \mathbf{K}\_{12,13} \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \begin{array}{cccc} \mu\_{\text{Opra},0}|\mathbf{K}| & \mu\_{\text{Opra},0}|\mathbf{K} \Big(\frac{\mathbf{R}\_{\text{S}}}{\mathbf{R}\_{\text{S}}}\Big)^{|\mathbf{K}|} & -\mu\_{VI}\mathbf{P}^{VII}\lambda^{VI}\Big(\frac{\mathbf{R}\_{\text{S}}}{\mathbf{R}\_{\text{T}}}\Big)^{\lambda^{VI}} & \mu\_{VI}\mathbf{P}^{VII}\lambda^{VII} \end{bmatrix} \end{bmatrix} \tag{A18}$$

$$
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{13,12} & \mathbf{K}\_{13,13} & \mathbf{K}\_{13,14} & \mathbf{K}\_{13,15}
\end{bmatrix}
\end{bmatrix} = 
\begin{bmatrix}
\mathbf{p}^{\text{VII}}\lambda^{\text{VIII}} & \mathbf{p}^{\text{VIII}}\left(\frac{\mathbf{R}\_{8}}{\mathbf{R}\_{7}}\right)^{\lambda^{\text{II}}} & -\mathbf{p}^{\text{VIII}}\left(\frac{\mathbf{R}\_{7}}{\mathbf{R}\_{8}}\right)^{\lambda^{\text{III}}} & -\mathbf{p}^{\text{VIII}}
\end{bmatrix}
\tag{A19}
$$

$$
\begin{bmatrix}
\begin{array}{cccc}
\mathbf{K}\_{112} & \mathbf{K}\_{113} & \mathbf{K}\_{114} & \mathbf{K}\_{115}
\end{bmatrix}
\end{bmatrix} - 
\begin{bmatrix}
\begin{bmatrix}
\mathbf{K}\_{12} & \mathbf{p}^{\text{III}}\lambda^{\text{III}} & \mathbf{p}^{\text{III}}\lambda^{\text{III}}
\end{bmatrix}
\end{bmatrix}
\end{bmatrix}
\begin{bmatrix}
\begin{array}{cccc}
\mathbf{k}^{\text{VII}} & \mathbf{k}^{\text{VIII}} \\
\mathbf{k}^{\text{II}}
\end{bmatrix}
\end{bmatrix}
\tag{A19}
$$

$$\begin{bmatrix} \mathbf{K}\_{4\cup 2} & \mathbf{K}\_{4\cup 3} & \mathbf{K}\_{4\cup 4} & \mathbf{K}\_{4\cup 5} \end{bmatrix} - \begin{bmatrix} \boldsymbol{\mu}\_{\text{Sat}\beta} \mathbf{P}^{\text{VII}} \boldsymbol{\lambda}^{\text{VII}} & -\boldsymbol{\mu}\_{\text{Sat}\beta} \mathbf{P}^{\text{III}} \boldsymbol{\lambda}^{\text{III}} \left(\frac{\mathbf{R}}{\mathbf{R}}\right)^{\mathbf{A}^{\text{IV}}} & -\boldsymbol{\mu}\_{\text{Opm}\beta} \mathbf{P}^{\text{III}} \boldsymbol{\lambda}^{\text{III}} \left(\frac{\mathbf{R}}{\mathbf{R}}\right)^{\mathbf{A}^{\text{IV}}} & \boldsymbol{\mu}\_{\text{Opm}\beta} \mathbf{P}^{\text{III}} \boldsymbol{\lambda}^{\text{III}} \end{bmatrix} \tag{A20}$$
 
$$\begin{bmatrix} \mathbf{K}\_{15,14} & \mathbf{K}\_{15,15} & \mathbf{K}\_{15,16} & \mathbf{K}\_{15,17} \end{bmatrix} = \begin{bmatrix} \mathbf{P}^{\text{VIII}} & \mathbf{P}^{\text{VIII}} \left(\frac{\mathbf{R}}{\mathbf{R}}\right)^{\mathbf{A}^{\text{III}}} & -\left(\frac{\mathbf{R}\_{\mathbf{q}}}{\mathbf{R}\mathbf{q}}\right)^{\mathbf{K} \mathbf{I}} & -\mathbf{I} \end{bmatrix} \tag{A21}$$

$$\begin{bmatrix} \begin{bmatrix} \mathbf{K}\_{16,14} & \mathbf{K}\_{16,15} & \mathbf{K}\_{16,16} & \mathbf{K}\_{16,17} \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \begin{array}{cc} \boldsymbol{\mu}\_{\text{IX}} \mathbf{P}^{\text{VII}} \boldsymbol{\lambda}^{\text{VIII}} & -\boldsymbol{\mu}\_{\text{IX}} \mathbf{P}^{\text{VIII}} \boldsymbol{\lambda}^{\text{VIII}} \left(\frac{\boldsymbol{\mu}\_{\text{7}}}{\boldsymbol{\lambda}\_{6}}\right)^{\text{IX}\text{II}} & \boldsymbol{\mu}\_{\text{Sil},0} \|\mathbf{K}\left(\frac{\boldsymbol{\mu}\_{\text{8}}}{\boldsymbol{\lambda}\_{6}}\right)^{\text{IK}} & \boldsymbol{\mu}\_{\text{Sil},0} \|\mathbf{K}\| \end{bmatrix} \end{bmatrix} \tag{A.22}$$

$$\begin{bmatrix} \mathbf{K}\_{17,16} & \mathbf{K}\_{17,17} & \mathbf{K}\_{17,18} \end{bmatrix} = \begin{bmatrix} \mathbf{I} & -\begin{pmatrix} \frac{R\_8}{R\_9} \end{pmatrix}^{\vert \mathbf{K} \rangle} & -\mathbf{I} \end{bmatrix} \tag{A23}$$

$$\begin{bmatrix} \begin{array}{cc} \mathbf{K}\_{18,16} & \mathbf{K}\_{18,17} & \mathbf{K}\_{18,18} \end{array} \end{bmatrix} = \begin{bmatrix} \begin{array}{cc} \boldsymbol{\mu}\_{\boldsymbol{X}} & -\boldsymbol{\mu}\_{\boldsymbol{X}} \begin{pmatrix} \frac{\boldsymbol{R}\_{8}}{\mathbf{K}\_{9}} \end{pmatrix}^{|\mathbf{K}|} & \boldsymbol{\mu}\_{\boldsymbol{I}\boldsymbol{X}} \end{bmatrix} \end{bmatrix} \tag{A24}$$

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
