Exploitation of independent stator and rotor geometrical periodicities in electrical machines using the Schur complement

Abstract

In this paper we present a-priori model reduction technique that enables to take full advantage of the periodicity existing in the stator and rotor geometrical structures of electrical machines in order to reduce the computational time. Firstly, a change of basis is performed by applying two distinct discrete Fourier transformations on the stator and rotor periodic structures independently. Secondly, the Schur complement is introduced in the new spectral basis, to allow a parallel solving of the resulting block-diagonal matrix systems. Moreover, the using of a matrix-free Krylov method based on the conjugate gradient solver has verified an efficient solving of the equation system associated to the stator-rotor interface. Furthermore, in the peculiar case of balanced supply conditions, a model order reduction can be carried out by considering only the dominant discrete Fourier transform components. This model reduction approach is applied on a buried permanent magnet machine and has successfully shown its efficiency under balanced and unbalanced conditions.

Keywords

Finite-element analyses model reduction geometrical periodicity discrete Fourier transform electromagnetic fields rotating electrical machines buried permanent magnet machines Schur complement

1. Introduction

The main challenge in numerical modeling is to search for a compromise between a good accuracy and a reduction in memory storage and computing time requirements. To advance towards this objective, the model reduction techniques become increasingly used, and this especially in electromagnetic fields computation. These can be a-posteriori model reduction techniques like proper orthogonal decomposition [1] and perturbation zooming methods [2] or a-priori reduction techniques like the ones dealing with the magnetic symmetry or further geometrical periodicity exhibited in most rotating electrical machines [3–5].

The difficulties encountered in the a-posteriori techniques remain the determination of a well-conditioned reduced basis as in the case of the model order reduction methods [6] and the association of equitable boundary conditions as in the case of most subdomain reduction methods [7]. In the a-priori techniques, we can list likewise several limitations. In fact, the reduction technique using magnetic symmetries takes advantage of the symmetrical magnetic flux lines distribution. Let’s take the example of a permanent magnet machine with N _s stator teeth and 2p rotor permanent magnets. In case of a kind of consistent distribution of the 2p permanent magnets with respect to the N _s stator teeth, which means a greatest common divisor (GCD) of N _s and 2p greater than 1, the modeling of a single 1∕GCD (N _s, 2p) section is enough to determine effectively the full model solution by simple periodic or anti-periodic relationships [8]. Nevertheless, such magnetic symmetry can be exploited only in case of balanced supply conditions. The geometrical periodicity on the other hand, considers only the materials permeability regardless the sources distribution and the supply regime whether it is under balanced or unbalanced conditions. In the case of a buried permanent magnet machine (BPMM) for example with N _s and N _r stator and rotor teeth respectively, one can model only 1∕GCD (N _s, N _r) machine section [9]. The modeling of this single elementary section which can regenerate the full geometrical model by simple rotation transformations, is enough to determine the full model solution if one uses the group representation theory [3,4]. However, such classical exploitation of the geometrical periodicity does not allow to take advantage of the total geometrical periodicity existing separately in the stator and rotor domains.

To tackle this issue, we investigate in this paper a-priori model reduction technique that takes full advantage of the geometrical periodicities existing in the stator and rotor structures, and goes further than the classical approach by allowing to only model a single 1∕N _s elementary section of the stator with a single 1∕N _r elementary section of the rotor. When it deals with the symmetrical properties of the block circulant matrices, this reduction technique is akin to the group representation theory [10,11] and can be applied regardless of the sources distribution. In a first step, this reduction technique is implemented by means of two distinct discrete Fourier transforms (DFT) applied to the two distinct periodic structures of the stator and rotor domains. In a second step when the rotor motion is taken into account, the Schur complement is used to solve the matrix system obtained from the coupling of both stator and rotor problems [12,13]. We should denote that the continuous rotor motion is performed using the spatial Fourier interpolation method (SFIM) [14,15] which is an extension of the locked step method [16]. We should also point out that the presented approach in this paper is applied to reduce the computational time in the case of linear finite element analysis. In fact, the homogeneous magnetic property of the materials with a linear behavior retains effectively the geometrical model periodicity. Nevertheless, such geometrical periodicity is broken in the case of a non-linear analysis. In the non-linear case, some techniques have been proposed to keep the block circulant property by using iterative algorithms based on the fixed-point technique or the transmission-line method. Such techniques allow not only to apply the DFT but also to account for the non-linearity by introducing an additional fictitious source term on the righthand side of the equation system [5,17].

The proposed model reduction technique is applied to study a 9 teeth/8 poles BPMM. In fact, this particular machine does not present neither any magnetic symmetry nor any consistent geometrical periodicity between the stator and rotor structures, and thus requires a full modeling with the classical approach.

2. Classical exploitation of geometrical periodicity in the finite element (FE) modeling

2.1. Linear magnetostatic FE modeling

Let us consider a domain 𝛺 of boundary 𝛤 (𝛤 = Γ₁ ∪ Γ₂ and ${\rm\Gamma}_{1}\cap {\rm\Gamma}_{2}=\{∅\}$ ). In magnetostatics, the Maxwell’s equations that describe the problem and the associated medium relationship are given by the following: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \text{curl}∼H=J_{s},\quad \text{div}∼B=0\end{eqnarray}$ (1) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle H={\nu}(B-B_{r})\end{eqnarray}$ (2) with B the magnetic flux density, H the magnetic field, J _s the excitation source current density, B _r the residual magnetization in the permanent magnets, and 𝜈 the magnetic reluctivity property of the materials. To insure the uniqueness of the solution, Neumann and Dirichlet boundary conditions must be considered on Γ₁ and Γ₂: $\begin{eqnarray}\displaystyle H\wedge n=0\quad \text{on }{\rm\Gamma}_{1};\qquad B.n=0\quad \text{on }{\rm\Gamma}_{2} & & \displaystyle\end{eqnarray}$ (3) with n the outward unit normal vector.

To solve the previous problem, the magnetic potential A (B = curl A) is introduced. The unicity of A is imposed with the addition of a gauge condition; this can be a Coulomb gauge which is implicitly verified in our case of 2D modeling. Adding now the magnetic potential to (1) and (2), the new resulting equation to solve will be: $\begin{eqnarray}\displaystyle \text{curl}∼{\nu}(\text{curl}∼A)=\text{curl}∼{\nu}B_{r}+J_{s}. & & \displaystyle\end{eqnarray}$ (4)

2.2. Incidence of the geometrical periodicity in the FE modeling problems

Let us generalize and consider the case of an electromagnetic device made up of N periodic sections (Fig. 1). We suppose that the periodicity holds only on the permeability but not on the source produced by the currents or the residual magnetic flux density.

Fig. 1.

Full model (left), modeling of one elementary section: the crossed points represent the inner nodes, the round points represent the ones common with the previous section and the squared points represent the ones subjected to Dirichlet conditions (right).

Since the partial differential equation (4) cannot be determined analytically, solving technique by FE discretization is used. Therefore, to take into account the geometrical periodicity we have to fulfill a numerical periodicity homologous to the geometrical one. To do so, we impose that all the sections are discretized by the same mesh and have thus the same number and distribution of nodes. To achieve this, it is sufficient to discretize only the elementary periodic section, as the one presented in Fig. 1, in such a way that the full model can be deduced by a simple concatenation of the elementary model. n is supposed to be the number of the degrees of freedom in each section. It represents the total number of inner nodes excepting in the one hand, the nodes located on the boundary where a Dirichlet condition is imposed and in the other hand, the ones located on the common side with the previous section (Fig. 1).

Now, the geometrical periodicity on the material characteristics of the N sections is in agreement with the numerical discretization periodicity. Taking into account the latter, and applying a 2D FE formulation of the linear magnetostatic problem described by the equation (4), a linear system of matrix equations is obtained: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle \left(\begin{array}{@{}ccccccccccc@{}}S & & S^{\prime } & & 0 & & \cdots ∼ & & 0 & & S^{\prime \prime }\\ S^{\prime \prime } & & S & & S^{\prime } & & 0 & & \cdots ∼ & & 0\\ \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots \\ \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots \\ \cdots \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots \\ \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots ∼ & \cdots \\ S^{\prime } & & 0 & & \cdots ∼ & & 0 & & S^{\prime \prime } & & S\end{array}\right)\left(\begin{array}{@{}c@{}}A_{0}\\ A_{1}\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ A_{N-1}\\ \end{array}\right)=\left(\begin{array}{@{}c@{}}F_{0}\\ F_{1}\\ \vdots \\ \vdots \\ \vdots \\ \vdots \\ F_{N-1}\end{array}\right)\end{eqnarray}$ (5) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [S_{c}][A]=[F]\end{eqnarray}$ (6) with [S _c] a block circulant matrix, [A] the vector of nodal magnetic potential unknowns of size ([N × n] × 1) and [F] the source vector. The matrix [S] of size (n × n) account for the inner nodes contribution of each section i (0 ≤ i ≤ N −1), [S′′] accounts for the contribution of the previous section (i −1) on the section i due to the nodes located on the common boundary between the two sections and [S ^′], likewise, accounts for the contribution of the next section (i + 1). In each section i, [A _i] and [F _i] of size (n × 1) represent respectively the vector of nodal magnetic potential unknowns and the source vector.

2.3. Model reduction by exploitation of geometrical periodicity

In this section, we will spell out the incidence of the geometrical periodicity on the model reduction. In fact, since [S _c] in (6) is a block circulant matrix, it can be decomposed, by using a DFT, under the following form [18,19]: $\begin{eqnarray}\displaystyle [S_{c}]=[W][S_{{\rm\Delta}}][W]^{-1}. & & \displaystyle\end{eqnarray}$ (7) [W] is the change of basis matrix from the original spatial basis into the spectral one and [S _Δ] is block diagonal. They are given with their following explicit forms: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [W]=[U_{N}]\otimes [I_{dn}]\end{eqnarray}$ (8) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [S_{{\rm\Delta}}]=\left(\begin{array}{@{}ccccccc@{}}\begin{array}{@{}c@{}}S_{0}\\ 0\\ 0\\ \vdots \\ \vdots \\ \vdots \\ 0\end{array} & \begin{array}{@{}c@{}}0\\ S_{1}\\ 0\\ \vdots \\ \vdots \\ \vdots \\ \cdots \end{array} & \begin{array}{@{}c@{}}\cdots \\ 0\\ S_{2}\\ \vdots \\ \vdots \\ \vdots \\ \cdots \end{array} & \begin{array}{@{}c@{}}\cdots \\ \cdots \\ 0\\ \vdots \\ \vdots \\ \vdots \\ \cdots \end{array} & \begin{array}{@{}c@{}}\cdots \\ \cdots \\ \cdots \\ \vdots \\ \vdots \\ \vdots \\ \cdots \end{array} & \begin{array}{@{}c@{}}\cdots \\ \cdots \\ \cdots \\ \vdots \\ \vdots \\ \vdots \\0\end{array} & \begin{array}{@{}c@{}}0\\ 0\\ 0\\ \vdots \\ \vdots \\ 0\\ S_{(N-1)}\end{array}\end{array}\right)\end{eqnarray}$ (9) with [I _dn] the identity matrix of size (n × n), and ⊗ the Kroenecker product symbol. The DFT matrix [U _N] of size (N × N) is given by its general term [20]: $\begin{eqnarray}\displaystyle U_{N}(c,l)=\frac{1}{\sqrt{N}}e^{j\frac{2{\pi}}{N}(c-1)(l-1)};\quad 1\leq c\leq N\text{ and }1\leq l\leq N. & & \displaystyle\end{eqnarray}$ (10)

Finally, in the spectral basis the new system to solve can be deduced combining (6) and (7): $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [S_{{\rm\Delta}}][Z]=[C],\quad \text{with}\end{eqnarray}$ (11) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z]=[W]^{-1}[A]\quad \text{and}\quad [C]=[W]^{-1}[F].\end{eqnarray}$ (12) [A] and [F], which are defined in the spatial basis, are represented now in the spectral basis by their corresponding vectors [Z] and [C] respectively. In the equation (11), [S _Δ] is a block diagonal matrix; the original system (6) of size n × N is therefore transformed into a set of N independent subsystems of size n. We don’t longer need to solve for the full model solution [A] but for the vector [Z] = [(Z _i)_0≤i≤N] by solving independently the following N reduced subsystems: $\begin{eqnarray}\displaystyle [S_{i}][Z_{i}]=[C_{i}];\quad 0\leq i\leq N-1. & & \displaystyle\end{eqnarray}$ (13)

Moreover, solving these N independent subsystems can be done using a parallelization process. Accordingly, it will lead to a considerable reduction of the computational time. Finally, we can deduct from (12) the scattered full model solution on the N different sections; it is given by [A] = [W][Z].

Fig. 2.

The 9/8 BPMM: highlighting the different sources (top left), considering the permeability of materials regardless of sources (bottom left), the two elementary cells of the stator and rotor domains (right).

Now, let’s apply the previously described method exploiting the geometrical periodicity on the modeling of an electrical machine example. We take therefore a BPMM consisting of N _s = 9 stator teeth wound by a three phase winding, and 2p = 8 rotor poles (Fig. 2). From the point of view of the distribution of the magnetic field lines, even when the stator windings are fed with a balanced supply conditions, we can expect that this machine will not present any magnetic symmetry since the GCD between N _s and 2p is equal to 1 (GCD (N _s = 9, 2p = 8) = 1) [8,9]. Therefore, we can denote that, for reasons of simplifications, the stator windings can be left unloaded. On the other side considering now the materials permeability point of view regardless the excitation sources distribution, no global consistent stator-rotor geometrical periodicity is exhibited since the rotor salient structure represents N _r = 8 rotor teeth -due to the presence of 8 buried permanent magnets- which means in fact that GCD (N _s = 9, N _r = 8) = 1 [9]. Based on the forgoing, there is neither a magnetic symmetry nor a classical consistent geometrical periodicity as the one described in Fig. 1 and therefore the full machine should be considered for a magnetostatic finite element modeling. However, we can distinguish two independent geometrical periodicities in the stator and rotor structures regardless of the sources distribution (current sources and residual magnetic flux densities). (Fig. 2). Consequently, we find 9 identical periodic cells in the stator involving each a single stator tooth and 8 identical periodic cells in the rotor representing each one rotor pole. These two elementary cells are represented in Fig. 2. In the next section 3, we will show how we can totally exploit such independent periodicities existing in the most stator and rotor domains of rotating electrical machines. This will be applied to study the 9/8 BPMM considered example, later in Section 4.

3. Total exploitation of distinct stator and rotor geometrical periodicities in rotating machines

3.1. Methodology

We suppose now the case of a rotating electromagnetic device, consisting of a fixed stator and a movable rotor, and disposing – as it is in the general case – of a distinct periodicity in the stator and the rotor geometrical structures.

Fig. 3.

Full model of a rotating electromagnetic device: The marked fictive interface separates the fixed part (stator) from the rotating part (rotor) (left), modeling of one elementary section (right).

We assume to perform the rotor motion using the SFIM [14,15] which is an extension of the locked step method [16]. To do so, the presence of a fictive interface between the stator and the rotor is assigned (Fig. 3). Since the full model of the device represents N _s identical periodic sections in the stator and N _R identical periodic sections in the rotor, the geometrical periodicity, as it has been developed before, is now exploited in the stator and rotor domains separately; the elementary model is consisting of one 1∕N _s elementary section of the stator and another 1∕N _r elementary section of the rotor (Fig. 3).

In Fig. 3, we can notice that in each elementary section of the electromagnetic device we can distinguish the elementary stator domain 𝛺_si (0 ≤ i ≤ N _s − 1) whose number of degrees of freedom is n _s from the one of the rotor domain 𝛺_rj (0 ≤ i ≤ N _r − 1) whose number of degrees of freedom is n _r. The stator-rotor interface, where the motion is considered is duplicated to create dual Γ_si and Γ_rj separation interfaces which belong respectively to 𝛺_si and 𝛺_ri. The relationship between the discretization nodes of the different interfaces Γ_si and Γ_rj allows, according to the locked steps method [16] or the SFIM [15], to consider the rotor motion. Based on the foregoing, we can summarize that there are N _s periodic sections in the stator, the degree of freedom of each is n _s and N _R periodic sections in the rotor, and the degree of freedom of each is n _r. The use of the geometrical periodicity which makes it possible to construct the stator and the rotor full models of the electromagnetic device independently by concatenation of the stator and rotor elementary models leads to obtain, according to the equations (5) and (6), the following redundant overdetermined system to solve: $\begin{eqnarray}\displaystyle \left(\begin{array}{@{}cc@{}}[S_{cs}] & [0]\\ \text{}[0] & [S_{cr}]\end{array}\right)\left(\begin{array}{@{}c@{}}[A_{s}]\\ \text{}[A_{r}]\end{array}\right)=\left(\begin{array}{@{}c@{}}[F_{s}]\\ \text{}[F_{r}]\end{array}\right). & & \displaystyle\end{eqnarray}$ (14)

Since the problem is treated independently within the stator and the rotor subdomains, the matrix system (14) is overdetermined thanks to the duplicated nodes of the stator-rotor interface; the equations relating to these nodes are written twice: the first for their relation with the inner nodes of the stator and the second for their relation with the inner nodes of the rotor. [S _cs] and [S _cr] are two block circulant matrices linked to stator and rotor geometrical subdomains respectively. [A _s] and [A _r] which are respectively of size ([N _s × n _s] × 1) and ([N _s × n _r] × 1), represent the set of nodal magnetic potential unknowns in the stator and the rotor, and [F _s] and [F _r] are the source terms. The vectors [A _s], [A _r], [F _s] and [F _r] are given with their following explicit forms: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[A_{s}]=\left(\begin{array}{@{}c@{}}(A_{ss})_{0}\\ (A_{{\rm\Gamma}_{s}})_{0}\\ (A_{ss})_{1}\\ (A_{{\rm\Gamma}_{s}})_{1}\\ \vdots \\ \vdots \\ (A_{ss})_{(N_{s}-1)}\\ (A_{{\rm\Gamma}_{s}})_{(N_{s}-1)}\end{array}\right);\qquad [A_{r}]=\left(\begin{array}{@{}c@{}}(A_{rr})_{0}\\ (A_{{\rm\Gamma}_{r}})_{0}\\ (A_{rr})_{1}\\ (A_{{\rm\Gamma}_{r}})_{1}\\ \vdots \\ \vdots \\ (A_{rr})_{(N_{r}-1)}\\ (A_{{\rm\Gamma}_{r}})_{(N_{r}-1)}\end{array}\right)\\[54.0pt] \text{}[F_{s}]=\left(\begin{array}{@{}c@{}}(F_{ss})_{0}\\ (F_{{\rm\Gamma}_{s}})_{0}\\ (F_{ss})_{1}\\ (F_{{\rm\Gamma}_{s}})_{1}\\ \vdots \\ \vdots \\ (F_{ss})_{(N_{s}-1)}\\ (F_{{\rm\Gamma}_{s}})_{(N_{s}-1)}\end{array}\right);\qquad [F_{r}]=\left(\begin{array}{@{}c@{}}(F_{rr})_{0}\\ (F_{{\rm\Gamma}_{r}})_{0}\\ (F_{rr})_{1}\\ (F_{{\rm\Gamma}_{r}})_{1}\\ \vdots \\ \vdots \\ (F_{rr})_{(N_{r}-1)}\\ (F_{{\rm\Gamma}_{r}})_{(N_{r}-1)}\end{array}\right)\end{array} & & \displaystyle\end{eqnarray}$ (15) with $[(A_{s})_{i}]=[(A_{ss})_{i}^{T},(A_{{\rm\Gamma}_{s}})_{i}^{T}]^{T}$ represents the set of the total nodal magnetic potential unknowns in the stator domain 𝛺_si of the section i. (A _ss)_i represents the set of nodal magnetic potential unknowns in the stator domain excepting the ones on the interface Γ_Si which are represented by the vector (A _{Γ_s})_i. $[(F_{s})_{i}]=[(F_{ss})_{i}^{T},(F_{{\rm\Gamma}_{s}})_{i}^{T}]^{T}$ is the source term calculated in the stator domain 𝛺_si. (F _ss)_i is the source term calculated in the stator domain but referring to the nodes belonging to the stator domain excepting those on the interface Γ_si whose associated source term is given by (F _{Γ_s})_i. The rotor domain vectors (A _r)_j, (A _rr)_j, (A _{Γ_r})_j, (F _r)_j, (F _rr)_j and (F _{Γ_r})_j are homologous to the previously defined stator domain vectors (A _s)_i, (A _ss)_i, (A _{Γ_s})_i, (F _s)_i, (F _ss)_i and (F _{Γ_s})_i.

According to (7), the block circulant matrices [S _cs] and [S _cr] can be decomposed, using two distinct DFT, under the following forms: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[S_{cs}]=[W_{s}][S_{{\rm\Delta}s}][W_{s}]^{-1}\\ \text{}[S_{cr}]=[W_{r}][S_{{\rm\Delta}R}][W_{r}]^{-1}.\end{array} & & \displaystyle\end{eqnarray}$ (16) [S _Δs] and [S _Δr] are two block diagonal matrices and the change of basis matrices [W _s] and [W _r] are given respectively in terms of the DFT matrices [U _{N
_s}] and [U _{N
_r}] as follow: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [W_{s}]=[U_{N_{s}}]\otimes [I_{dn_{s}}]\end{eqnarray}$ (17) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [W_{r}]=[U_{N_{r}}]\otimes [I_{dn_{r}}]\end{eqnarray}$ (18) where [I _{dn
_s}] of size (n _s × n _s) and [I _{dn
_r}] of size (n _r × n _r) are the identity matrices and [U _{N
_s}] and [U _{N
_r}] are the DFT matrices of size (N _s × N _s) and (N _r × N _r) respectively.

Combining now the matrix system (14) with the equations in (16), we obtain in the spectral basis the new system to solve: $\begin{eqnarray}\displaystyle \left(\begin{array}{@{}cc@{}}[S_{{\rm\Delta}s}] & [0]\\ \text{}[0] & [S_{{\rm\Delta}r}]\end{array}\right)\left(\begin{array}{@{}c@{}}[Z_{s}]\\ \text{}[Z_{r}]\end{array}\right)=\left(\begin{array}{@{}c@{}}[C_{s}]\\ \text{}[C_{r}]\end{array}\right) & & \displaystyle\end{eqnarray}$ (19) where the vectors [Z _s], [Z _r], [C _s] and [C _r] are given by the following expressions: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z_{s}]=[W_{s}]^{-1}[A_{s}]\quad \text{and}\quad [C_{s}]=[W_{s}]^{-1}[F_{s}]\end{eqnarray}$ (20) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z_{r}]=[W_{r}]^{-1}[A_{r}]\quad \text{and}\quad [C_{r}]=[W_{r}]^{-1}[F_{r}].\end{eqnarray}$ (21)

Since [S _Δs] is a block diagonal matrix, we have therefore N _s independent sub-systems of size n _s each and since [S _Δr] is a block diagonal matrix, we have likewise N _r independent subsystems of size n _r each: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[(S_{{\rm\Delta}s})_{i}][(Z_{s})_{i}]=[(C_{s})_{i}];\quad 0\leq i\leq N_{s}-1\\ \text{}[(Z_{s})_{i}]=\left[\begin{array}{@{}c@{}}(Z_{ss})_{i}\\ (Z_{{\rm\Gamma}_{s}})_{i}\end{array}\right];\qquad [(C_{s})_{i}]=\left[\begin{array}{@{}c@{}}(C_{ss})_{i}\\ (C_{{\rm\Gamma}_{s}})_{i}\end{array}\right]\end{array} & & \displaystyle\end{eqnarray}$ (22) $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[(S_{{\rm\Delta}r})_{j}][(Z_{r})_{j}]=[(C_{r})_{j}];\quad 0\leq j\leq N_{r}-1\\ \text{}[(Z_{r})_{j}]=\left[\begin{array}{@{}c@{}}(Z_{rr})_{j}\\ (Z_{{\rm\Gamma}_{r}})_{j}\end{array}\right];\qquad [(C_{r})_{j}]=\left[\begin{array}{@{}c@{}}(C_{rr})_{j}\\ (C_{{\rm\Gamma}_{r}})_{j}\end{array}\right].\end{array} & & \displaystyle\end{eqnarray}$ (23)

3.2. Consideration of the rotor motion using SFIM

To properly solve the overall problem consisting of these different independent subsystems in the stator and the rotor domains, they must be coupled together at the stator-rotor interface. To do so, we have to express the different [(Z _{Γ_r})_j] as functions of the different [(Z _{Γ_s})_i]. In other words, where [Z _{Γ_s}] and [Z _{Γ_r}] are respectively the set of [(Z _{Γ_s})_i] and [(Z _{Γ_r})_j] we must express [Z _{Γ_r}] in terms of [Z _{Γ_s}]. This can be done by performing the rotor motion using the SFIM. It allows considering any rotation angle θ of the rotor and not only the discrete few angle steps given by θ_m = mΔθ, where $m \in \mathbb{Z}$ and Δθ is a previously fixed rotation angle step.

Unlike the locked step method which takes into account discrete rotation steps [16], the SFIM allows therefore to take into account the continuous rotor motion. Assuming that [A _{Γ_S}] = [(A _{Γ_s})_{0≤i≤N
_s−1}] and [A _{Γ_r}] = [(A _{Γ_r})_{0≤j≤N
_r−1}] are the vectors of total magnetic potential unknowns on the stator-rotor interface in the stator and rotor side respectively, the SFIM consists in expressing the DFT of both vectors [A _{Γ_S}] and [A _{Γ_r}] under the following equation [15]: $\begin{eqnarray}\displaystyle [Z_{F_{{\rm\Gamma}_{r}}}]=[D({\theta})][Z_{F_{{\rm\Gamma}_{s}}}] & & \displaystyle\end{eqnarray}$ (24) where [D (θ)] is the diagonal matrix linking the DFT [Z _{F
_{Γ_R}}] of [A _{Γ_r}] to the DFT [Z _{F
_{Γ_s}}] of [A _{Γ_s}]. It allows to add a phase lag of amplitude θ to a signal, by multiplying its kth harmonic by e ^jθk. In our case, θ represents in fact the rotor mechanical position. When m _D is the total number of the stator-rotor interface nodes, the entries of the matrix [D (θ)] of size (m _D × m _D) are given by the following expression [15]: $\begin{eqnarray}\displaystyle D(i,i)=\left\{\begin{array}{@{}l@{}}e^{j{\theta}(i-1)};\quad i∼{\epsilon}∼[1,m_{D}^{+}]\\ e^{j{\theta}(i-1-m_{D})};\quad i∼{\epsilon}∼[m_{D}^{+}+1,m_{D}]\end{array}\right.\quad \text{where }m_{D}^{+}=[(m_{D}+1)/2]. & & \displaystyle\end{eqnarray}$ (25)

Now, given that [W _F] of size (m _D × m _D) is the transformation matrix corresponding to the DFT based on the local discretization of the stator-rotor interface, the vectors [Z _{F
_{Γ_s}}] and [Z _{F
_{Γ_r}}] are given by the following equations: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z_{F_{{\rm\Gamma}_{s}}}]=[W_{F}]^{-1}[A_{{\rm\Gamma}_{s}}]\end{eqnarray}$ (26) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z_{F_{{\rm\Gamma}_{r}}}]=[W_{F}]^{-1}[A_{{\rm\Gamma}_{r}}].\end{eqnarray}$ (27)

We should denote furthermore that the vectors [A _{Γ_s}] and [A _{Γ_r}] can be expressed in function of their DFT vectors [Z _{Γ_s}] and [Z _{Γ_r}] by the following expressions respectively: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z_{{\rm\Gamma}_{s}}]=[W_{{\rm\Gamma}_{s}}]^{-1}[Z_{{\rm\Gamma}_{s}}]\end{eqnarray}$ (28) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [Z_{{\rm\Gamma}_{r}}]=[W_{{\rm\Gamma}_{r}}]^{-1}[Z_{{\rm\Gamma}_{r}}]\end{eqnarray}$ (29) where: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [W_{{\rm\Gamma}_{s}}]=[U_{N_{s}}]\otimes [I_{d(\frac{m_{D}}{N_{s}})}]\end{eqnarray}$ (30) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [W_{{\rm\Gamma}_{r}}]=[U_{N_{r}}]\otimes [I_{d(\frac{m_{D}}{N_{r}})}].\end{eqnarray}$ (31)

Combining now the equations (24), (26), (27), (28) and (29) leads to the following expression: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[Z_{{\rm\Gamma}_{r}}]=[W_{{\rm\Gamma}_{r}}]^{-1}[W_{F}][D({\theta})][W_{F}]^{-1}[W_{{\rm\Gamma}_{s}}][Z_{{\rm\Gamma}_{S}}]\\ \text{}[Z_{{\rm\Gamma}_{r}}]=[Q({\theta})][Z_{{\rm\Gamma}_{s}}].\end{array} & & \displaystyle\end{eqnarray}$ (32)

In the spectral basis, [Q (θ)] is the matrix linking the interface variables of the stator to those of the rotor at the rotor position θ and then defining the continuous rotor motion. The matrix [D (θ)] being diagonal, therefore the matrix [W _F][D (θ)][W _F]⁻¹ is circulant. However, in our case of a non-consistent periodicity in the stator and the rotor, [W _{Γ_s}] ≠ [W _{Γ_r}], and therefore the matrix [Q (θ)] = [W _{Γ_r}]⁻¹[W _F][D (θ)][W _F]⁻¹[W _{Γ_s}] is not a block diagonal matrix. In the following and for reasons of simplification we will replace [Q (θ)] by [Q].

The matrix [Q] being not a block diagonal matrix, each sub-vector [(Z _{Γ_r})_j] must be expressed as a function of all the vectors [(Z _{Γ_s})_i]. This will not have any advantage if one wants to have independent subproblems resulting from the stator-rotor coupling. For this reason, we will proceed to tackle this issue by the use of the Schur complement method.

3.3. Using of the Schur complement

The application of the Schur complement actually requires a particular form of the matrix to be solved. This indeed requires a reordering of the variables in the vectors [Z _s] and [Z _r]. After such reordering of the variables, these two vectors will be written henceforth by the following expressions: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[Z_{s}]=\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\end{array}\right)=\left(\begin{array}{@{}c@{}}(Z_{ss})_{0}\\ \vdots \\ \vdots \\ (Z_{ss})_{(N_{s}-1)}\\ (Z_{{\rm\Gamma}_{s}})_{0}\\ \vdots \\ \vdots \\ (Z_{{\rm\Gamma}_{s}})_{(N_{s}-1)}\end{array}\right)\\[54.0pt] \text{}[Z_{r}]=\left(\begin{array}{@{}c@{}}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{r}}]\end{array}\right)=\left(\begin{array}{@{}c@{}}(Z_{rr})_{0}\\ \vdots \\ \vdots \\ (Z_{rr})_{(N_{r}-1)}\\ (Z_{{\rm\Gamma}_{r}})_{0}\\ \vdots \\ \vdots \\ (Z_{{\rm\Gamma}_{r}})_{(N_{r}-1)}\end{array}\right).\end{array} & & \displaystyle\end{eqnarray}$ (33)

Using the expressions of the vectors [Z _s] and [Z _r] in (33) by separating the variables associated within the inner domain from those located on the interface, the matrix system in the equation (19) can be now written in the following form: $\begin{eqnarray}\displaystyle \left(\begin{array}{@{}cccc@{}}\begin{array}{@{}c@{}}[S_{{\rm\Delta}ss}]\\ 0\\ \text{}[S_{{\rm\Delta}s{\rm\Gamma}_{s}}]\\ 0\end{array} & \begin{array}{@{}c@{}}0\\ \text{}[S_{{\rm\Delta}rr}]\\ 0\\ \text{}[S_{{\rm\Delta}r{\rm\Gamma}_{r}}]\end{array} & \begin{array}{@{}c@{}}[S_{{\rm\Delta}{\rm\Gamma}_{s}s}]\\ 0\\ \text{}[S_{{\rm\Delta}{\rm\Gamma}_{s}}]\\ 0\end{array} & \begin{array}{@{}c@{}}0\\ \text{}[S_{{\rm\Delta}{\rm\Gamma}_{r}r}]\\ 0\\ \text{}[S_{{\rm\Delta}{\rm\Gamma}_{r}}]\end{array}\end{array}\right)\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\\ \text{}[Z_{{\rm\Gamma}_{r}}]\end{array}\right)=\left(\begin{array}{@{}c@{}}[C_{ss}]\\ \text{}[C_{rr}]\\ \text{}[C_{{\rm\Gamma}_{s}}]\\ \text{}[C_{{\rm\Gamma}_{r}}]\end{array}\right) & & \displaystyle\end{eqnarray}$ (34) $\begin{eqnarray}\displaystyle [S_{T}]\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\\ \text{}[Z_{{\rm\Gamma}_{r}}]\end{array}\right)=\left(\begin{array}{@{}c@{}}[C_{ss}]\\ \text{}[C_{rr}]\\ \text{}[C_{s}]\\ \text{}[C_{{\rm\Gamma}_{r}}]\end{array}\right) & & \displaystyle \nonumber\end{eqnarray}$ where the bock-diagonal matrices [S _Δss], [S _{ΔΓ_s s}], [S _{ΔsΓ_s}] and [S _{ΔΓ_s}] and the vectors [C _ss] and [C _{Γ_s}] can be deduced respectively from the block diagonal matrix [S _Δs] and the vector [C _s] in (19) by means of simple separation and reordering of the variables according to (33). Same explanation is used for the determination of the matrices [S _Δrr], [S _{ΔΓ_r r}], [S _{ΔrΓ_r}] and [S _{ΔΓ_r}] and the vectors [C _rr] and [C _{Γ_r}] from the matrix [S _Δr] and the vector [C _r]. The equation (34) represents the full model problem which is characterized by an overdetermined algebraic system due to the duplication of the variable on the stator-rotor interface, as we have mentioned before. To tackle this issue in order to obtain a consistent algebraic equation system, we can introduce the equation (32); [Z _{Γ_r}] and [Z _{Γ_s}] are dependent and therefore a change of basis matrix [P], which is given in the following, must be used in the system of equation (34): $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\\ \text{}[Z_{{\rm\Gamma}_{r}}]\end{array}\right)=[P]\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\end{array}\right)\\ \left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\\ \text{}[Z_{{\rm\Gamma}_{r}}]\end{array}\right)=\left(\begin{array}{@{}ccc@{}}\begin{array}{@{}c@{}}[I_{dn_{s}}]\\ 0\\ 0\\ 0\end{array} & \begin{array}{@{}c@{}}0\\ \text{}[I_{dn_{r}}]\\ 0\\ 0\end{array} & \begin{array}{@{}c@{}}0\\ 0\\ \text{}[I_{dm}]\\ \text{}[Q]\end{array}\end{array}\right)\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\end{array}\right).\end{array} & & \displaystyle\end{eqnarray}$ (35)

Using now the equation (35) in (34), and multiplying by the transposed matrix of [P] in order to add the contribution from both sides of the stator-rotor interface, the new matrix equation system that yields will be therefore: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[P]^{T}[S_{T}][P]\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\end{array}\right)=[P]^{T}\left(\begin{array}{@{}c@{}}[C_{ss}]\\ \text{}[C_{rr}]\\ \text{}[C_{{\rm\Gamma}_{s}}]\\ \text{}[C_{{\rm\Gamma}_{r}}]\end{array}\right)\\[27.60004pt] \left(\begin{array}{@{}ccc@{}}\begin{array}{@{}c@{}}[S_{{\rm\Delta}ss}]\\[5.0pt] 0\\ \text{}[S_{{\rm\Delta}s{\rm\Gamma}_{s}}]\end{array} & \begin{array}{@{}c@{}}0\\ \text{}[S_{{\rm\Delta}rr}]\\[5.0pt] \text{}[Q]^{T}[S_{{\rm\Delta}r{\rm\Gamma}_{r}}]\end{array} & \begin{array}{@{}c@{}}[S_{{\rm\Delta}{\rm\Gamma}_{s}s}]\\[5.0pt] \text{}[S_{{\rm\Delta}{\rm\Gamma}_{r}r}][Q]\\[5.0pt] \text{}[S_{{\rm\Delta}{\rm\Gamma}_{s}}]+[Q]^{T}[S_{{\rm\Delta}{\rm\Gamma}_{r}}][Q]\end{array}\end{array}\right)\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\end{array}\right)=\left(\begin{array}{@{}c@{}}[C_{ss}]\\ \text{}[C_{rr}]\\ \text{}[C_{{\rm\Gamma}_{s}}]+[Q]^{T}[C_{{\rm\Gamma}_{r}}]\end{array}\right).\end{array} & & \displaystyle\end{eqnarray}$ (36)

This particular type of matrix system equation can be solved using the Schur complement [12,13]. This can be done by using a simple substitution method. Let’s rewrite therefore the first and the second equation in (36) by determining the expressions of [Z _ss] and [Z _rr] in terms of [Z _{Γ_s}] as the following: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[Z_{ss}]=[S_{{\rm\Delta}ss}]^{-1}([C_{ss}]-[S_{{\rm\Delta}{\rm\Gamma}_{s}s}][Z_{{\rm\Gamma}_{s}}])\\ \text{}[Z_{rr}]=[S_{{\rm\Delta}rr}]^{-1}([C_{rr}]-[S_{{\rm\Delta}{\rm\Gamma}_{r}r}][Q][Z_{{\rm\Gamma}_{s}}]).\end{array} & & \displaystyle\end{eqnarray}$ (37)

Now replacing [Z _ss] and [Z _rr] by their values given by the equations in (37), in the third equation of the full system (36), we will obtain therefore the following system to solve whose size is equal to the number of nodes on the interface, and which is very small in comparison with the initial full problem.: $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle [S_{sr}][Z_{{\rm\Gamma}_{s}}]=[C_{sr}]\nonumber\end{eqnarray}$ $\begin{eqnarray}\displaystyle \text{}[S_{sr}]=[S_{{\rm\Delta}{\rm\Gamma}_{s}}]+[Q]^{T}[S_{{\rm\Delta}{\rm\Gamma}_{r}}][Q]-[S_{{\rm\Delta}s{\rm\Gamma}_{s}}][S_{{\rm\Delta}ss}]^{-1}[S_{{\rm\Delta}{\rm\Gamma}_{s}s}]-[Q]^{T}[S_{{\rm\Delta}r{\rm\Gamma}_{r}}][S_{{\rm\Delta}rr}]^{-1}[S_{{\rm\Delta}{\rm\Gamma}_{r}r}][Q] & & \displaystyle\end{eqnarray}$ (38) $\begin{eqnarray}\displaystyle [C_{sr}]=[C_{{\rm\Gamma}_{s}}]+[Q]^{T}[C_{{\rm\Gamma}_{r}}]+[S_{{\rm\Delta}s{\rm\Gamma}_{s}}][S_{{\rm\Delta}ss}]^{-1}[C_{ss}]+[Q]^{T}[S_{{\rm\Delta}r{\rm\Gamma}_{r}}][S_{{\rm\Delta}rr}]^{-1}[C_{rr}]. & & \displaystyle \nonumber\end{eqnarray}$

At the end, the Schur complement will lead us to the following system that can be solved in the opposite direction now with a back substitution method: $\begin{eqnarray}\displaystyle \left(\begin{array}{@{}ccc@{}}\begin{array}{@{}c@{}}[S_{{\rm\Delta}ss}]\\ 0\\ 0\end{array} & \begin{array}{@{}c@{}}0\\ \text{}[S_{{\rm\Delta}rr}]\\ 0\end{array} & \begin{array}{@{}c@{}}[S_{{\rm\Delta}{\rm\Gamma}_{s}s}]\\ \text{}[S_{{\rm\Delta}{\rm\Gamma}_{r}r}][Q]\\ \text{}[S_{sr}]\end{array}\end{array}\right)\left(\begin{array}{@{}c@{}}[Z_{ss}]\\ \text{}[Z_{rr}]\\ \text{}[Z_{{\rm\Gamma}_{s}}]\end{array}\right)=\left(\begin{array}{@{}c@{}}[C_{ss}]\\ \text{}[C_{rr}]\\ \text{}[C_{sr}]\end{array}\right). & & \displaystyle\end{eqnarray}$ (39)

The determination of [Z _{Γ_s}] by solving the system (38) in a first step allows us to calculate [Z _{Γ_r}] using the equation (32) and to determine [Z _ss] and [Z _rr] in a second step and this by solving respectively the two following matrix equations: $\begin{eqnarray}\displaystyle \begin{array}{@{}l@{}}[S_{{\rm\Delta}ss}][Z_{ss}]=[C_{ss}]-[S_{{\rm\Delta}{\rm\Gamma}_{s}s}][Z_{{\rm\Gamma}_{s}}]\\ \text{}[S_{{\rm\Delta}rr}][Z_{rr}]=[C_{rr}]-[S_{{\rm\Delta}{\rm\Gamma}_{r}r}][Z_{{\rm\Gamma}_{r}}].\end{array} & & \displaystyle\end{eqnarray}$ (40)

3.4. Efficient solver using the conjugate gradient method and an implicit LU decomposition

We should denote that, to solve (38), we can avoid the implicit calculation of the matrix [S _sr]. This can be done using a matrix-free Krylov method based on the conjugate gradient solver [21]. It is an iterative method that does not deal with the direct solving of the system (38) but with an iteration process of matrix-vector multiplications. At each iteration, the matrix [S _sr] is multiplied by the kth iteration vector denoted [y _k].

Since [S _Δss] and [S _Δrr] are block-diagonal matrices with N _s and N _r matrix blocks respectively; the system (40) will lead to (N _s + N _r) independent subsystems that can be solved effectively in a parallel computation process. Moreover, this parallel computation can be used also to solve for [S _Δss]⁻¹[C _ss] and [S _Δrr]⁻¹[C _rr] within the calculation of [C _sr] and to solve for [S _Δss]⁻¹[S _{ΔΓ_s s}][y _k] and [S _Δrr]⁻¹[S _{ΔΓ_r r}][Q][y _k] within each conjugate gradient solver for the solution [Z _{Γ_s}] of the system (38). Henceforth, the solving of the matrix system (39) using a back substitution method combined with a parallelization process can lead to an effective acceleration in the computational time which will be verified in the following example.

To prevent a repetitive direct solving technique, we recommend the use of an implicit LU decomposition. The LU decomposition uses a backward and forward solving i.e. triangular solving and has verified therefore a more efficient calculation time. If we consider that (S _Δss)_i and (S _Δrr)_j are the ith (0 ≤ i ≤ N _s − 1) and the jth (0 ≤ j ≤ N _r −1) matrix blocks of the block diagonal matrices [S _Δss] and [S _Δrr] respectively, we can write: (S _Δss)_i = (L _ss)_i × (U _ss)_i and (S _Δrr)_j = (L _rr)_j × (U _rr)_j. These LU decomposition matrices are calculated one time and are stored; this avoids a significant waste of the computational time in the repetitive direct solving technique that will be requested for example in the conjugate gradient iterations in (38). All these steps to solve the system (39) can be summarized in the following diagram in Fig. 4. We should denote that the boxes with dashed borders represent the ith and the jth systems of the N _s and N _r systems respectively, that could be solved in a parallel process. Furthermore, and for reasons of simplification, the conjugate gradient solver box is represented only by the matrix- vector multiplication of the matrix [S _sr] with the iterative vector [y _k] at the kth iteration.

Fig. 4.

Diagram showing the different solving steps and highlighting the parallel computational processes presented by the boxes with dashed borders.

4. Application example

The application example is a 9/8 BPMM (Fig. 5). It consists of 9 stator teeth and 8 rotor permanent magnets. As we have mentioned in Section 2, there is neither an obvious magnetic symmetry nor a whole consistent geometrical periodicity in this studied machine. As a matter of fact, the full machine should be modeled in the classical approach. However, the stator and the rotor represent two different geometrical periodic structures. In fact, from an elementary 1/9 section of the stator and an elementary 1/8 section of the rotor we can generate the full model of the machine. Therefore, to exploit the total geometrical periodicity as we have noticed in Section 2, the elementary periodic cell in the stator turns out to be only one stator tooth and that of the rotor turns out to be only one permanent magnet (Fig. 5).

Fig. 5.

The BPMM 9/8 (left), Mesh of the elementary modeled sections: elementary 1/9 stator section and elementary 1/8 rotor section (right).

We should denote that the 2D spatial mesh of the elementary cell is made of n _s = 738 nodes in the stator section and n _r = 512 nodes in the rotor section, and that for reasons of simplifications, as we have mentioned in Section 2, the stator phases are left unloaded. From the mesh of one elementary cell we have reconstructed a full FE model (reference) that we have compared to the reduced model. The exploitation of the geometrical periodicity makes it possible to switch effectively from the large FE system into 9 independent subsystems in the stator, 8 independent subsystems in the rotor and one system representing the coupling of the stator-rotor nodes on the connecting interface.

4.1. Balanced regime

The first application is the case of a balanced permanent magnet source distribution in the rotor while the stator winding is kept unloaded. In fact, we have verified that the reduced model solution matches largely the one of the reference full model. The post processing global quantity compared between both models is the magnetic flux flowing through a tooth coil (Fig. 6). The local quantities are the flux lines distribution (Fig. 7), the magnetic flux density (Fig. 8) and the magnetic potential calculated at the air gap level in function of the angular position (Fig. 9).

Fig. 6.

Magnetic flux through a tooth coil calculated with the reference full model and the reduced model.

Fig. 7.

Flux lines distribution calculated with the reference full model (left) and the reduced model (right).

Fig. 8.

Flux density distribution calculated with the reference full model (left) and the reduced model (right).

Fig. 9.

Magnetic potential in the air gap in function of the angular position calculated with the reference full model and the reduced model.

The flux lines, the flux density distributions and the magnetic potential are calculated at a given rotor position while the magnetic flux is calculated for a full rotor mechanical revolution. The conjugate gradient solver converges after 72 iterations with a stop criterium defined by a maximum residual error of 10⁻³. The speed up resulting from the transformation of the large system into the spectral basis system is 6 when using a sequential computation and 11 if we use a parallelization process. The errors between the full and the reduced models are 0.15% in the calculation of the magnetic flux and 0.3% for the air gap magnetic potential. This error is due to the use of a matrix-free krylov iteration method which is characterized by a maximum number of iterations linked to a stopping criterium that insures a minimum residual error.

4.2. Order reduction in balanced regimes considering only Z _s4 and Z _r4

Each component of the stator and the rotor DFT components actually represents a given spectral content of the solution distribution in the stator and the rotor domain respectively (see Appendix B). According to (55) the harmonic spectrum of the stator and the rotor solution distribution can be projected to the corresponding DFT components as presented in Tables 1 and 2. The conjugate relationship between the different DFT components is explained in Appendix A. In our case, the distribution of the source in the rotor has a particular shape: it consists of 4 pairs of permanent magnets leading to a source distribution in such a way that the 4th harmonic is dominant. This can be verified by looking at the spectral representation of the magnetic potential distribution in the airgap (Fig. 10).

Table 1
Harmonic content of the stator DFT components

Z _s0 [l × N _s]⁺ = [0, 9, 18, 27, …]

Z _s1 = Z _s8 [l × N _s − 1, l × N _s + 1]⁺ = [1, 8, 10, …]

Z _s2 = Z _s7 [l × N _s − 2, l × N _s + 2]⁺ = [2, 7, 11, …]

Z _s3 = Z _s6 [l × N _s − 3, l × N _s + 3]⁺ = [3, 6, 12, …]

Z _s4 = Z _s5 [l × N _s − 4, l × N _s + 4]⁺ = [4, 5, 13, …]

Table 2

Harmonic content of the rotor DFT components

Z _r0	[l × N _r]⁺ = [0, 8, …]
Z _r1 = Z _r7	[l × N _r − 1, l × N _r + 1]⁺ = [1, 7, 9, …]
Z _r2 = Z _r6	[l × N _r − 2, l × N _r + 2]⁺ = [2, 6, 10, …]
Z _r3 = Z _r5	[l × N _r − 3, l × N _r + 3]⁺ = [3, 5, 11, …]
Z _r4	[l × N _r − 4, l × N _r + 4]⁺ = [4, 12, 18, …]

Fig. 10.

Spectral representation of the magnetic potential distribution at the air gap level with respect to the associated stator and rotor DFT components.

In Tables 1 and 2, $l \in \mathbb{N}$ , []⁺ represents the positive value and Z is the complex conjugate of Z. According to Tables 1 and 2, this prevailing 4th harmonic corresponds effectively to the stator DFT components Z _s4 and the rotor DFT component Z _r4. Now solving only the equation sub-system giving Z _s4, instead of the all 9 sub-systems in the stator, and the equation sub-system giving Z _r4, instead of the all 8 sub-systems in the rotor, can be seen as an effective model order reduction with a speed up of 11 even with a sequential solver, avoiding therefore the use of a complex parallelization process.

The results are in agreement not only with respect to the global quantity of the magnetic flux (Fig. 11) where the error does not exceed 1.8% but also regarding the local ones (Figs 12, 13, 14 and 15), but now with a slight low precision where the error considering the calculation of the magnetic potential in the air gap which is about 11% from the stator side (Fig. 14) and 14% from the rotor side (Fig. 15).

Fig. 11.

Magnetic flux through a tooth coil calculated with the reference full model and the reduced model with only Z _s4 and Z _r4.

Fig. 12.

Flux lines distribution calculated with the reference full model (left) and the reduced model with only Z _s4 and Z _r4 (right).

Fig. 13.

Flux density distribution calculated with the reference full model (left) and the reduced model with only Z _s4 and Z _r4 (right).

Fig. 14.

Magnetic potential in the air gap in function of the angular position calculated with the reference full model and the reduced model with only Z _s4.

Fig. 15.

Magnetic potential in the air gap in function of the angular position calculated with the reference full model and the reduced model with only Z _r4.

4.3. Unbalanced regime with a demagnetization of a rotor permanent magnet

In this section we will study the case of an unbalanced regime characterized by a demagnetization defect of a permanent magnet in the rotor (Fig. 16). In this case, we can no longer speak about a particular distribution of the solution in the machine and thus of the dominance of a single spatial harmonic as it was the case in the balanced regime. Therefore, all the DFT components in the stator and the rotor domain have to be calculated. The conjugate gradient solver converges after 97 iterations with the same stop criterium, as in balanced conditions, defined by a maximum residual error of 10⁻³. The number of iterations until convergence is not the same in the case of balanced and unbalanced conditions since the convergence depends on the source vector. Now the speed-up is slightly smaller and it is about 5 in the case of a sequential computation of the problem and about 10.2 when we use the parallel computation. The solution calculated by the reduced model is compared to the one of a classical full model in the case of the demagnetization defect: they are in a good agreement as we can see in Figs 17, 18, 19, and 20.

Fig. 16.

Demagnetization defect of a permanent magnet in the rotor.

Fig. 17.

Magnetic flux through a tooth coil calculated with the reference full model and the reduced model.

Fig. 18.

Flux lines distribution calculated with the reference full model (left) and the reduced model (right).

Fig. 19.

Flux density distribution calculated with the reference full model (left) and the reduced model (right).

Fig. 20.

Magnetic potential in the air gap in function of the angular position calculated with the reference full model and the reduced model.

The error does not exceed the 0.15% with respect to the calculation of the magnetic flux in a full mechanical revolution (Fig. 17) and the 0.3% with respect to the calculation of the magnetic potential in the airgap (Fig. 20).

The harmonic spectrum of the magnetic potential distribution in the airgap in Fig. 21 shows that this distribution is not governed by a single spatial harmonic but by different harmonics of considerable values. The solution in the stator and the rotor as well as the magnetic potential in the airgap are calculated separately in function of each stator and rotor DFT components and are represented in Figs 22 and 23.

Fig. 21.

Spectral representation of the magnetic potential distribution at the air gap level with respect to the associated stator and rotor DFT components.

Fig. 22.

Flux lines distribution and air gap magnetic potential calculated in the stator by means of the different modes Z _s0, Z _s1, Z _s2, Z _s3 and Z _s4.

Fig. 23.

Flux lines distribution and air gap magnetic potential calculated in the rotor by means of the different modes Z _r0, Z _r1, Z _r2, Z _r3 and Z _r4.

5. Conclusion

In this paper we have investigated a-priori model reduction approach based on the full exploitation of the spatial periodicity exhibited in the stator and rotor geometrical structures of the most rotating electrical machines. In fact, the use of the Schur complement has enhanced such exploitation of the geometrical periodicity by leading to solve in parallel the block diagonal matrix systems in the stator and rotor domains. Moreover, an efficient solving of the stator-rotor coupling system on the interface is carried out by means of the matrix-free Krylov method based on the conjugate gradient solver. A 9/8 BPMM is studied in balanced and unbalanced regimes; it has a particular geometrical structures where no magnetic symmetry neither a classical consistent geometrical periodicity can be used to reduce the FE modeling. However, using the presented model reduction approach, we have successfully managed to solve the matrix system with a speed up of about 6 in case of a sequential computation and 11 when we use a parallel process.

Furthermore, in the particular case of balanced supply conditions, we have verified a model order reduction of the FE model due to the DFT properties; only one DFT component that represent the prevailing harmonic of the solution distribution can be only solved for, which lead to an effective acceleration in the calculation time of about 11 within even a sequential computation of the matrix system.

Footnotes

Conjugate relationship between DFT components

Assuming the case of a vector X of dimension N given by the following form: (41) $\begin{eqnarray}\displaystyle [X]=[\begin{array}{@{}cccccc@{}}x_{0} & x_{1} & \cdots ∼ & x_{m} & \cdots ∼ & x_{N-1}\end{array}]^{T}. & & \displaystyle\end{eqnarray}$

The DFT vector [Z] of [X] will be given by the following expression using the inverse DFT matrix [U]⁻¹: (42) $\begin{eqnarray}\displaystyle [Z]=[U]^{-1}[X]=[\begin{array}{@{}cccccc@{}}Z_{0} & Z_{1} & \cdots ∼ & Z_{m} & \cdots ∼ & Z_{N-1}\end{array}]^{T} & & \displaystyle\end{eqnarray}$ where the entries of [U]⁻¹ are given by: (43) $\begin{eqnarray}\displaystyle U^{-1}(l,c)=\frac{1}{\sqrt{N}}e^{-j\times \frac{2{\pi}}{N}\times (l-1)\times (c-1)}\quad 1\leq l\leq N\text{ and }1\leq c\leq N. & & \displaystyle\end{eqnarray}$

From the equations (42) and (43), the general term Z _m of the vector [Z] will be given therefore by the following expression: (44) $\begin{eqnarray}\displaystyle Z_{m}=\mathop{\sum }_{n=0}^{N-1}x_{n}e^{-j\times \frac{2{\pi}}{N}\times m\times n}. & & \displaystyle\end{eqnarray}$

Let us first suppose the case where N, the number of samples is even, then the vector [Z] can be given by its general representation form: (45) $\begin{eqnarray}\displaystyle [Z]=\left[\begin{array}{@{}c@{}}Z_{0}\\ Z_{\frac{N}{2}-{\alpha}}\\[8.0pt] Z_{\frac{N}{2}}\\ Z_{\frac{N}{2}+{\alpha}}\end{array}\right]\quad \text{where }{\alpha}\in \left[1;2\ldots ;\frac{N}{2}-1\right]. & & \displaystyle\end{eqnarray}$

Based on the equation (44), the expressions of both general terms $Z_{\frac{N}{2}-{\alpha}}$ and $Z_{\frac{N}{2}+{\alpha}}$ will be given by: (46) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle Z_{\frac{N}{2}-{\alpha}}=\mathop{\sum }_{n=0}^{N-1}x_{n}e^{-j\times {\pi}\times n}\times e^{j\times \frac{2{\pi}}{N}\times {\alpha}\times n}\end{eqnarray}$ (47) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle Z_{\frac{N}{2}+{\alpha}}=\mathop{\sum }_{n=0}^{N-1}x_{n}e^{-j\times {\pi}\times n}\times e^{-j\times \frac{2{\pi}}{N}\times {\alpha}\times n}\end{eqnarray}$ with n being an integer which takes the values 1,2, …N −1, the term e ^−j×π×n is then a real number that can take the values 1 or −1. We can then deduce, from (46) and (47), that $Z_{\frac{N}{2}-{\alpha}}$ and $Z_{\frac{N}{2}+{\alpha}}$ are conjugate: (48) $\begin{eqnarray}\displaystyle Z_{\frac{N}{2}+{\alpha}}=\overline{Z_{\frac{N}{2}-{\alpha}}};\quad {\alpha}\in \left[1;2\ldots ;\frac{N}{2}-1\right]. & & \displaystyle\end{eqnarray}$

In the other case where the sampling number N is odd, the vector [Z] can be given its general representation form: (49) $\begin{eqnarray}\displaystyle [Z]=\left[\begin{array}{@{}c@{}}Z_{0}\\ Z_{\frac{N+1}{2}-{\alpha}}\\ Z_{\frac{N-1}{2}+{\alpha}}\end{array}\right]\quad \text{where }{\alpha}\in \left[1;2\ldots ;\frac{N}{2}-1\right]. & & \displaystyle\end{eqnarray}$

In the same way, based on the equation (44) the expressions of both general terms $Z_{\frac{N+1}{2}-{\alpha}}$ and $Z_{\frac{N-1}{2}+{\alpha}}$ will be given by: (50) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle Z_{\frac{N+1}{2}-{\alpha}}=\mathop{\sum }_{n=0}^{N-1}x_{n}e^{-j\times {\pi}\times n}\times e^{-j\times \frac{{\pi}}{N}\times (1-2{\alpha})\times n}\end{eqnarray}$ (51) $\begin{eqnarray}\displaystyle \text{} & \text{} & \displaystyle Z_{\frac{N-1}{2}+{\alpha}}=\mathop{\sum }_{n=0}^{N-1}x_{n}e^{-j\times {\pi}\times n}\times e^{j\times \frac{{\pi}}{N}\times (1-2{\alpha})\times n}.\end{eqnarray}$

We can then deduce from (50) and (51) that $Z_{\frac{N+1}{2}-{\alpha}}$ and $Z_{\frac{N-1}{2}+{\alpha}}$ are conjugate: (52) $\begin{eqnarray}\displaystyle Z_{\frac{N-1}{2}+{\alpha}}=\overline{Z_{\frac{N+1}{2}-{\alpha}}};\quad {\alpha}\in \left[1;2\ldots ;\frac{N-1}{2}\right]. & & \displaystyle\end{eqnarray}$

Spectral content of DFT components

Let’s consider the case of a signal x (t) representing a sinusoidal function of frequency kf corresponding to the kth harmonic of a reference signal of frequency f. When this given kth harmonic signal whose mathematical representation is x (t) = A _ksin(k2πft), is discretizing through a sampling frequency F _s = N ∗ f, the sampling vector consisting of the N samples will be given by the following form: (53) $\begin{eqnarray}\displaystyle [X]=\left[\begin{array}{@{}c@{}}0\\ A_{k}\sin (k2{\pi}/N)\\ \vdots \\ A_{k}\sin (k2{\pi}/N\times m)\\ \vdots \\ A_{k}\sin -(k2{\pi}/N\times (N-1))\end{array}\right]. & & \displaystyle\end{eqnarray}$

Performing now in a following step, the DFT of the vector [X] which leads to the vector [Z] = [U ⁻¹][X] whose general term, using equation (44), is given by the following expression: (54) $\begin{eqnarray}\displaystyle Z_{m}=\mathop{\sum }_{n=0}^{N-1}\frac{A_{k}}{2j\sqrt{N}}[e^{\frac{j\times 2{\pi}\times (k-m)}{N}\times n}-e^{-\frac{j\times 2{\pi}\times (k+m)}{N}\times n}]. & & \displaystyle\end{eqnarray}$

We can show that this series in (54) will be not equal to zero if there exists an integer l such that: (55) $\begin{eqnarray}\displaystyle k\pm m=N\times l. & & \displaystyle\end{eqnarray}$

In the excepted cases presented by the relation (55), $e^{\frac{j\times ∼2{\pi}\times (k-m)}{N}\times n}$ or $e^{-\frac{j\times ∼2{\pi}\times (k+m)}{N}\times n}$ will be equal to 1 regardless the value of n; Z _m will be not equal to zero and it will include therefore the spectral information on the harmonic k.

References

Chatterjee

, An introduction to the proper orthogonal decomposition, Current Science 78(7) (2000), 808–817.

Sabariego

R.V.

and Dular

, A Perturbation approach for the modeling of eddy current nondestructive testing problems with differential probes, IEEE Transactions on Magnetics 43(4) (2007), 1289–1292.

Bossavit

, The exploitation of geometrical symmetry in 3-D Eddy-currents computation, IEEE Transactions on Magnetics 21(6) (1985), 2307–2309.

Lobry

and Broche

, Exploitation of the geometrical symmetry in the boundary element method with the group representation theory, IEEE Transactions on Magnetics 30(1) (1994), 118–123.

Lobry

Trecat

and Broche

, Symmetry and TLM method in nonlinear magnetostatics, IEEE Transactions on Magnetics 32(3) (1996), 702–705.

Noor

A.K.

and Peters

J.M.

, Reduced basis technique for nonlinear analysis of structures, AIAA Journal 18(1) (1980), 455–462.

Al Eit

Dular

Bouillault

Marchand

and Krebs

, Perturbation finite element method for efficient copper losses calculation in switched reluctance machines, IEEE Transactions on Magnetics 53(6) (2017), 1–4.

Chandrupatla

T.R.

Belegundu

A.D.

Ramesh

and Ray

, Introduction to Finite Elements in Engineering, Prentice Hall, Upper Saddle River, NJ, 2002.

Larsen

Gensior

and Hein

, Torque ripple reduction based on current control for a flux switching permanent magnet machine, in: International Symposium on Power Electronics Power Electronics, Electrical Drives, Automation and Motion, Sorrento, 2012, pp. 712–717.

10.

Hamermesh

, Group Theory and its Application to Physical Problems, Addison-Wesley Publishing Company, 1962.

11.

Jansens

and Boon

, Theory of finite groups, Applications in Physics (1967).

12.

Haynsworth

E.V.

, On the Schur complement, Basel Mathematical Notes, #BNB 20, June 1968.

13.

Zhang

, The Schur Complement and its Applications, Springer, 2005.

14.

De Gersem

and Weiland

, Harmonic weighting functions at the sliding interface of a finite-element machine model incorporating angular displacement, IEEE Transactions on Magnetics 40(2) (2004), 545–548.

15.

Montier

Clenet

Henneron

and Goursaud

, Rotation movement based on the spatial Fourier interpolation method (SFIM), IEEE Transactions on Magnetics 53(6) (2017), 1–4.

16.

Preston

T.W.

Reece

A.B.J.

and Sangha

P.S.

, Induction motor analysis by time-stepping techniques, IEEE Transactions on Magnetics 24(1) (1988), 471–474.

17.

Chiampi

Chiarabaglio

and Repetto

, A Jiles-Atherton and fixed-point combined technique for time periodic magnetic field problems with hysteresis, IEEE Transactions on Magnetics 31(6) (1995), 4306–4311.

18.

Davis

P.J.

, Circulant Matrices, Chelsea Publishing N. Y., 1994.

19.

Olson

Shaw

Shi

Pierre

and Parker

, Circulant matrices and their application to vibration analysis, Applied Mechanics Reviews, American Society of Mechanical Engineers 66(4) (2014).

20.

Dickinson

and Steiglitz

, Eigenvectors and functions of the discrete Fourier transform, IEEE Transactions on Acoustics, Speech, and Signal Processing 30(1) (1982), 25–31.

21.

Hestenes

M.R.

and Stiefel

E.L.

, Methods of conjugate gradients for solving linear systems, Journal of Research of the National Bureau of Standards, Section B 49 (1952), 409–436.

Z _s0	[l × N _s]⁺ = [0, 9, 18, 27, …]
Z _s1 = Z _s8	[l × N _s − 1, l × N _s + 1]⁺ = [1, 8, 10, …]
Z _s2 = Z _s7	[l × N _s − 2, l × N _s + 2]⁺ = [2, 7, 11, …]
Z _s3 = Z _s6	[l × N _s − 3, l × N _s + 3]⁺ = [3, 6, 12, …]
Z _s4 = Z _s5	[l × N _s − 4, l × N _s + 4]⁺ = [4, 5, 13, …]