Model Order Reduction for Rotating Electrical Machines

TL;DR

This paper develops an adaptive model order reduction method using proper orthogonal decomposition for simulating non-symmetric rotating electrical machines efficiently, with error certification and demonstrated effectiveness.

Contribution

It introduces an adaptive reduction technique specifically tailored for non-symmetric machines, addressing challenges not covered by existing symmetric models.

Findings

Effective reduction of simulation costs for non-symmetric machines

Error estimator certifies solution accuracy

Numerical examples validate the method's efficiency

Abstract

The simulation of electric rotating machines is both computationally expensive and memory intensive. To overcome these costs, model order reduction techniques can be applied. The focus of this contribution is especially on machines that contain non-symmetric components. These are usually introduced during the mass production process and are modeled by small perturbations in the geometry (e.g., eccentricity) or the material parameters. While model order reduction for symmetric machines is clear and does not need special treatment, the non-symmetric setting adds additional challenges. An adaptive strategy based on proper orthogonal decomposition is developed to overcome these difficulties. Equipped with an a posteriori error estimator the obtained solution is certified. Numerical examples are presented to demonstrate the effectiveness of the proposed method.

Figures

Tables

Table 1. Table 1: Performance summary for the different settings and index set 𝒦 𝒦 \mathcal{K} .

Setting	FEM (sec.)	ROM (sec.)	Basis $({\Omega }_s,{\Omega }_r)$	Iter.	Speedup	Overhead
\svhline sym	$328.56$	$7.04$	$(12,6)$	$1$	$46.67$	$37$%
rot	$337.17$	$12.76$	$(12,6)$	$2$	$26.42$	$29$%
stat	$341.28$	$46.98$	$(13,16)$	$6$	$7.26$	$41$%
rot_stat	$328.30$	$49.85$	$(13,16)$	$6$	$6.58$	$47$%

Table 2. Table 2: Performance summary for the different settings and index set ℳ ℳ \mathcal{M} .

Setting	FEM (sec.)	ROM (sec.)	Basis $({\Omega }_s,{\Omega }_r)$	Iter.	Speedup	Overhead
\svhline sym	$328.56$	$14.36$	$(13,6)$	$2$	$22.88$	$38$%
rot	$337.17$	$20.60$	$(13,6)$	$3$	$16.36$	$34$%
stat	$341.28$	$32.74$	$(13,16)$	$4$	$10.42$	$44$%
rot_stat	$328.30$	$42.17$	$(13,16)$	$6$	$7.78$	$37$%

Equations

$$\vec{\nabla }\times \left(\nu \vec{\nabla }\times \vec{A}(\vartheta )\right)={\vec{J}}_{\mathrm{src}}(\vartheta )-\nabla \times \left(\nu {\mathbf{B}}_{\mathrm{rem}}\right),\text{mbox}on{\Omega }_{\mathrm{s}}$$

$$-\nabla \cdot (\nu \nabla A_z(\vartheta ))=j_{\mathrm{src}}(\vartheta )-j_{\mathrm{pm}},\text{mbox}on\Omega (\vartheta ),$$

$$A_z^{\mathrm{N}}=\sum _{i=1}^{\mathrm{N}}{a}_i{w}_i(x,y)=\sum _{i=1}^{\mathrm{N}}{a}_i\frac{{\phi }_i(x,y)}{{\ell }_z}{e}_z,$$

$${\mathbf{K}}_{\nu }(\vartheta )\mathbf{a}(\vartheta )=\mathbf{f}(\vartheta )$$

$$\left[\begin{array}{ccc}{\mathbf{K}}_{\nu }^{ss}& 0& {\mathbf{K}}_{\nu }^{sI}\\ 0& {\mathbf{K}}_{\nu }^{rr}& {\mathbf{K}}_{\nu }^{rI}(\vartheta )\\ ({\mathbf{K}}_{\nu }^{sI}{)}^{\top }& ({\mathbf{K}}_{\nu }^{rI}{)}^{\top }(\vartheta )& {\mathbf{K}}_{\nu }^{II}(\vartheta )\end{array}\right]\left[\begin{array}{c}{\mathbf{a}}_s\\ {\mathbf{a}}_r\\ {\mathbf{a}}_I\end{array}\right]=\left[\begin{array}{c}{\mathbf{f}}^s\\ {\mathbf{f}}^r\\ {\mathbf{f}}^I(\vartheta )\end{array}\right],$$

$${\mathbf{K}}_{\nu }(\vartheta )=\sum _{j=1}^L\Theta (\ell ){\mathbf{K}}_{\nu }^j\text{mbox}and\mathbf{f}(\vartheta )=\sum _{j=1}^L\Theta (\ell ){\mathbf{f}}^j,$$

$$\left\{\begin{array}[]{l}\displaystyle\min_{\psi_{1},\ldots,\psi_{\mathcal{N}}\in\mathbb{R}^{\mathrm{N}}}\sum_{k\in\mathcal{K}}\,\Big{|}\mathbf{a}^{k}(\mu)-\sum_{i=1}^{\mathcal{N}}{\langle\mathbf{a}^{k}(\mu),\psi_{i}\rangle}_{\mathrm{W}}\,\psi_{i}\Big{|}_{\mathrm{W}}^{2}\vspace{2mm}\\ \hskip 18.49428pt\text{s.t. }{\langle\psi_{i},\psi_{j}\rangle}_{\mathrm{W}}=\delta_{ij}\text{ for }1\leq i,j\leq\mathcal{N},\end{array}\right.$$

$$\mathbf{R}\psi=\sum_{k\in\mathcal{K}}\langle\mathbf{a}^{k}(\mu),\psi\rangle_{\mathrm{W}}\mathbf{a}^{k}=\mathbf{A}_{\mathcal{K}}\big{(}\mathbf{A}_{\mathcal{K}}\big{)}^{\top}\mathrm{W}\psi\quad\text{for }\psi\in\mathbb{R}^{\mathrm{N}}.$$

$${\mathbf{A}}_{\mathcal{K}}^{\top }\mathrm{W}{\mathbf{A}}_{\mathcal{K}}{v}_i={\lambda }_i{v}_i,i=1,\dots ,\mathcal{N}$$

$$\sum_{k\in\mathcal{K}}\,\Big{|}\mathbf{a}^{k}-\sum_{i=1}^{\mathcal{N}}{\langle\mathbf{a}^{k},\psi_{i}\rangle}_{\mathrm{W}}\,\psi_{i}\Big{|}_{\mathrm{W}}^{2}=\sum_{i=\mathcal{N}+1}^{d}\lambda_{i},$$

$${\mathbf{a}}_s^{\mathcal{N}}=\sum _{i=1}^{{\mathcal{N}}_s}{\psi }_i^s{\bar{\mathbf{a}}}_{s,i}={\Psi }^s{\bar{\mathbf{a}}}_s\text{and}{\mathbf{a}}_r^{\mathcal{N}}=\sum _{i=1}^{{\mathcal{N}}_r}{\psi }_i^r{\bar{\mathbf{a}}}_{r,i}={\Psi }^r{\bar{\mathbf{a}}}_r,$$

$$\left[\begin{array}{ccc}({\Psi }^s{)}^{\top }{\mathbf{K}}_{\nu }^{ss}{\Psi }^s& 0& ({\Psi }^s{)}^{\top }{\mathbf{K}}_{\nu }^{sI}\\ 0& ({\Psi }^r{)}^{\top }{\mathbf{K}}_{\nu }^{rr}{\Psi }^r& ({\Psi }^r{)}^{\top }{\mathbf{K}}_{\nu }^{rI}(\vartheta )\\ ({\mathbf{K}}_{\nu }^{sI}{)}^{\top }{\Psi }^s& ({\mathbf{K}}_{\nu }^{rI}{)}^{\top }(\vartheta ){\Psi }^r& {\mathbf{K}}_{\nu }^{II}(\vartheta )\end{array}\right]\left[\begin{array}{c}{\bar{\mathbf{a}}}_s\\ {\bar{\mathbf{a}}}_r\\ {\mathbf{a}}_I\end{array}\right]=\left[\begin{array}{c}({\Psi }^s{)}^{\top }{\mathbf{f}}^s\\ ({\Psi }^r{)}^{\top }{\mathbf{f}}^r\\ {\mathbf{f}}^I(\vartheta )\end{array}\right].$$

$${\mathbf{K}}_{\nu }^{\mathcal{N}}(\vartheta )\bar{\mathbf{a}}={\mathbf{f}}^{\mathcal{N}}(\vartheta )$$

$$\frac{{\sum }_{i=1}^{\mathcal{N}}{\lambda }_i}{{\sum }_{i=1}^d{\lambda }_i}\le {\epsilon }_{rel}$$

$$\alpha (\vartheta )=\underset{v\in {\mathbb{R}}^{\mathrm{N}}\setminus \{0\}}{\inf }\frac{{\mathbf{v}}^{\top }{\mathbf{K}}_{\nu }(\vartheta )\mathbf{v}}{{\mathbf{v}}^{\top }\mathbf{Wv}}.$$

$${\mathbf{K}}_{\nu }(\vartheta )\mathbf{v}=\lambda \mathbf{Wv}$$

$$\Vert \mathbf{a}-{\mathbf{a}}^{\mathcal{N}}{\Vert }_{\mathbf{W}}\le {\Delta }_{\mathbf{a}}(\vartheta ):=\frac{\Vert r({\mathbf{a}}^{\mathcal{N}};\vartheta ){\Vert }_{{\mathbf{W}}^{-1}}}{\alpha (\vartheta )}.$$

$$\frac{\Vert \mathbf{a}-{\mathbf{a}}^{\mathcal{N}}{\Vert }_{\mathbf{W}}}{\Vert {\mathbf{a}}^{\mathcal{N}}{\Vert }_{\mathbf{W}}}\le {\Delta }_{\mathbf{a}}^{\mathrm{rel}}(\vartheta ):=2\frac{\Vert r({\mathbf{a}}^{\mathcal{N}};\vartheta ){\Vert }_{{\mathbf{W}}^{-1}}}{\alpha (\vartheta )\Vert {\mathbf{a}}^{\mathcal{N}}{\Vert }_{\mathbf{W}}}.$$

$${\mathcal{K}}_{ij}=150(i-1)+[j,j+12,j+24,\dots ],i=1,\dots ,6\text{mbox}andj=1,\dots ,12.$$

$${\mathcal{M}}_i=[i,i+72,i+144,\dots ],i=1,\dots ,72.$$

Full text

1 Introduction

Model order reduction for partial differential equations is a very active field in applied mathematics. When performing simulations in $2$D or $3$D using the finite element method (FEM) one arrives at large scale systems that have to be solved. Projection based model order reduction methods have shown to significantly reduce the computational complexity when applied carefully. While being applied to many different fields in physics, the application to rotating electrical machines is more recent HC14 ; LU16 ; MHCG16 ; SSSI16 . We will focus especially on the setting of non-symmetric machines. While the perfect machine is symmetric and simulations are usually carried out exploiting these properties, in real life the symmetry is often lost. This is due to perturbation in the geometry (e.g., eccentricity) and material properties and requires that the whole machine is simulated and not only a small portion (e.g., one pole). Hence this leads to an increase in the computational cost. The aim is to develop an adaptive strategy that is able to collect the required information systematically. Ideally, the algorithm is able to detect symmetries and exploits them if present. The greedy algorithm introduced in the context of the reduced basis method is a possible candidate PR06 ; QMN16 . A commonly used method in engineering and applied mathematics is the method of snapshots, or proper orthogonal decomposition (POD) Cha00 ; HLBR12 ; Sir87 . We opt for a combination of the two methodologies. The goal is a fast and efficient variant that avoids the necessity of an online-offline decomposition. This allows a broader application since no expensive offline costs have to be invested. Hence the method will already pay off after one simulation and not only in the many query context.

Additionally, the developed strategy has to be able to handle the motion of the rotor. While there are a number of methods to treat the rotation DRL85 ; PC95 ; SLDP06 , we will assume a constant rotational speed which allows us to utilize the locked step method PRS88 . Hence we can avoid the remeshing which is computationally prohibitive. Moreover, the application of other approaches should be straight forward.

Efficient simulation tools are a key ingredient when performing optimization or uncertainty quantification. The combination of model order reduction and optimization has caught a lot of attention DH15 ; GV17 ; NGMQ13 ; QGVW17 ; ZF15 . Especially in the many query context, where models have to be evaluated repeatedly and the need for accurate, fast and reliable reduced order models is high. While we will not look into the application of the reduced order models, we will develop a strategy that fulfills these needs. By using existing simulation tools in the adaptive procedure it is possible to insert the developed strategy into an existing framework and utilize the benefits of the reduced order model, as shown in LU16 .

The article is structured as follows: In Section 2 the permanent magnet synchronous machine (PMSM) is introduced. The discretization and the proposed model order reduction strategy are introduced in Section 3. Then in Section 4 numerical experiments are presented. Lastly, a conclusion is drawn in the last section.

2 Problem Description

The PMSM under investigation has six slots per pole and a double layered winding scheme with two slots per pole per phase. The geometry of the full machine is shown in Figure 1 (left) together with a detailed view on one pole (right). In each pole there is one buried permanent magnet indicated in gray. The machine has depth $\ell_{z}=10\,\mathrm{mm}$. It is operated at $50\,\mathrm{Hz}$, resulting in a rotational speed of $1000$ revolutions per minute (RPM). The machine is composed of laminated steel with a relative permeability $\mu_{r}=500$.

In the following $\Omega_{\mathrm{s}}$ and $\Omega_{\mathrm{r}}$ refer to the stator and rotor domain, respectively. The whole domain is then referred to as $\Omega=\Omega_{\mathrm{s}}\cup\Omega_{\mathrm{r}}$ and will be used for simpler notation when appropriate. Additionally, let us define the interface $\Gamma_{\mathrm{I}}=\partial\Omega_{\mathrm{s}}\cap\partial\Omega_{\mathrm{r}}$ (dashed line) in the airgap between the rotor and the stator. Furthermore, we introduce the boundaries $\Gamma_{\mathrm{s}}=\partial\Omega_{\mathrm{s}}\setminus\Gamma_{\mathrm{I}}$ and $\Gamma_{\mathrm{r}}=\partial\Omega_{\mathrm{r}}\setminus\Gamma_{\mathrm{I}}$ of the stator and rotor, correspondingly. In the simulation we will account for the movement of the rotor, hence we introduce the angle $\vartheta$ that describes the position of the rotor with respect to the stator. For clarity we will append $\vartheta$ to the components related to the rotor.

To calculate magnetic vector potential of the machine, the magnetostatic approximation of the Maxwell’s equations has to be solved for both domains. This implies that the eddy and displacement currents are neglected and one obtains the semi-elliptic partial differential equations

[TABLE]

with Dirichlet boundary conditions $\vec{A}\times\mathbf{n}=0$ on $\Gamma_{\mathrm{s}}$, with $\mathbf{n}$ the outer unit normal. The reluctivity is depicted by a scalar $\nu$ since only linear isotropic materials are considered and nonlinearity is disregarded since a linearization at a working point is assumed. $\vec{A}(\vartheta)$ is the magnetic vector potential, $\vec{J}_{\rm src}(\vartheta)$ represents the imposed source current density, which is related to the applied currents in the coils, and $\mathbf{B}_{\rm rem}$ the remanence of the permanent magnets. The applied current density is aligned with the $z$-direction, whereas the remanence is in the $xy$-plane. It is generally accepted that machines are adequately modeled in $2$D, meaning that the magnetic field has no $z$-component: $\vec{B}=(B_{x},B_{y},0)$. Since $\vec{B}=\nabla\times\vec{A}$, one can write $\vec{A}=(0,0,A_{z})$. Hence we end up with the linear elliptic equation

[TABLE]

where $j_{\mathrm{src}}$ and $j_{\mathrm{pm}}$ are the $z$-component of $\vec{J}_{\rm src}$ and $\nabla\times\left(\nu\mathbf{B}_{\rm rem}\right)$, respectively. Further, the boundary conditions are given as previously introduced by $A_{z}=0$ on $\Gamma_{\mathrm{s}}$.

Let us next have a look at possible imperfections in the presented geometry and model. These are introduced during mass production of the PMSMs and are given by small deviations in the geometry or material properties. While there are many possible properties that can occur, we focus on two types. On the one hand we look at imperfection in the material of the permanent magnet, more precisely we consider deviations in the magnetic field angle $\phi$ of the permanent magnet OH13 . The second imperfection we consider is the length $\ell$ of the teeth in the stator Cle13 ; OMNCDH15 . Both quantities are depicted in Figure 1 (right). Perturbations in these quantities may lead to underperformance of the PMSM.

3 Model order reduction

In this section we discuss a method to accelerate the simulation of the PMSM. We will start by introducing the finite element discretization. Further, the realization of the rotation is outlined for the discrete setting. The resulting linear system of equations are of large scale and hence expensive to solve. In a next step we will present a model order reduction method based on proper orthogonal decomposition to speedup the simulations.

3.1 Finite element discretization

We obtain the discrete version of (2) by utilizing the finite element method (FEM). Discretizing $A_{z}$ by linear edge shape functions $w_{i}(x,y)$ one makes the Ansatz

[TABLE]

where $\varphi_{i}(x,y)$ depicts the nodal finite elements which are associated with the triangulation of the machine’s cross-section and $e_{z}$ is the unit vector in $z$-direction. Using the Ritz-Galerkin approach, the $\mathrm{N}$-dimensional linear discrete system

[TABLE]

is obtained, where $\mathbf{K}_{\nu}$ are the finite element system matrices, $\mathbf{a}$ depict the degrees of freedom (DoFs) and $\mathbf{f}=\mathbf{j}_{\mathrm{src}}-\mathbf{j}_{\mathrm{pm}}$ the discretized versions of the current densities and permanent magnets.

To take the motion of the rotor into account in the simulation, we utilize the locked step method PRS88 ; SLDP06 . For the implementation, a contour in the airgap is defined ($\Gamma_{\mathrm{I}}$) which splits the full domain into two parts: a fixed outer domain connected to the stator $\Omega_{\mathrm{s}}$ and an inner domain connected to the rotor $\Omega_{\mathrm{r}}$, where the mesh will be rotated. At the contour the nodes are distributed equidistantly. The angular increments $\Delta\vartheta$ are chosen so that the mesh of the stator and rotor will match on the interface. The nodes on the interface are then reconnected leading to the mesh for the next computation. Using this technique, the rotation angle $\vartheta$ is given by $\vartheta^{k}=k\Delta\vartheta$ with $k\in\mathbb{N}_{0}$. We can hence partition the discrete unknown $\mathbf{a}$ in into a static part, a rotating part and the interface, with dimensions $\mathrm{N}_{s}$, $\mathrm{N}_{r}$ and $\mathrm{N}_{I}$, respectively. This idea is a particularization of non-overlapping domain decomposition AHH11 ; TW05 . The linear system (3) can then be written as

[TABLE]

where $\mathbf{K}^{ss}_{\nu}$, $\mathbf{K}^{rr}_{\nu}$, $\mathbf{f}^{s}$ and $\mathbf{f}^{r}$ are the stiffness matrices and right hand sides on the static and moving part and do no longer depend on the angle $\vartheta$. For the points on the interface there are two cases. The interface of the static part is independent of the angle $\vartheta$ and hence we get the corresponding stiffness matrix $\mathbf{K}^{sI}_{\nu}$. For the rotor side we have to perform the shift, this is indicated with $\vartheta$ in the corresponding stiffness matrix $\mathbf{K}^{rI}_{\nu}$. On the interface also a shift has to be performed hence also here the corresponding stiffness matrix $\mathbf{K}^{II}_{\nu}$ and right hand side $\mathbf{f}^{I}$ are dependent on $\vartheta$. Let us note that it is not required to reassemble matrices. All of these shifts can be performed by index shift and hence allow a very efficient implementation. The size of the system does not change, i.e., we have $\mathrm{N}=\mathrm{N}_{s}+\mathrm{N}_{r}+\mathrm{N}_{I}$.

For completeness let us remark on the domain deformation introduced by $\ell$. Since we are in the setting of a parametrized shape transformation, i.e., one $\ell$ for each tooth in the stator we can utilize an affine geometry preconditioning to describe the transformation. This allows to carry out the computation on a reference domain and hence avoids again expensive remeshing and its undesired consequences as e.g., mesh noise. The stiffness matrices and right hand sides can be written in the form

[TABLE]

where $\mathbf{K}_{\nu}^{j}$ and $\mathbf{f}^{j}$ are local matrices. For a detailed description including computational procedures, we refer the reader to RHP08 .

3.2 Proper orthogonal decomposition

The goal is to generate a reduced order model to accelerate the simulation of (4). The simulation of the rotation is computationally expensive since the discretization of (2) leads to very large linear systems that need to be solved for every angular position. While in symmetric machines this can be avoided, in the case of non-symmetric machines a whole revolution has to be simulated. Hence we require an efficient strategy to overcome this problem. For this we investigate an adaptive strategy that builds a surrogate model while performing the simulation and switches to it when the required accuracy is reached. We want to use information collected over the rotational angle $\vartheta$ and generate a projection based reduced order model. In the past, model order reduction methods based on proper orthogonal decomposition (POD) Cha00 ; HLBR12 ; KV01 , balanced truncation AHH11 ; HSS08 and the reduced basis method PR06 ; QMN16 ; RHP08 have been developed to speedup the computation. More recently, the POD method has been successfully applied to rotating machines HC14 ; LU16 ; MHCG16 ; SSSI16 . In this work we consider a combination of POD and the reduced basis method. We will not pursue an online-offline decomposition but rather see the reduced order model as an accelerator for the simulation.

Let us start by recalling the POD method so we can develop an extension suitable for the application presented. Let the solution to (2) be given in the discrete form, i.e., let $\mathbf{a}(\vartheta)\in\mathbb{R}^{\mathrm{N}}$ be the solution to (3) for a fixed angle $\vartheta$. The snapshots are then given by $\mathbb{R}^{\mathrm{N}}\ni\mathbf{a}^{k}\approx\mathbf{a}(k\Delta\vartheta)$ for $k\in\mathcal{K}$, where $\mathcal{K}$ is an index set with elements in $\mathbb{N}_{0}$ for which (3) is solved. A POD basis $\{\psi_{i}\}_{i=1}^{\ell}$ is then computed from these snapshots by solving the following optimization problem:

[TABLE]

where $\langle\cdot\,,\cdot\rangle_{\mathrm{W}}$ stands for the weighted inner product in $\mathbb{R}^{\mathrm{N}}$ with a positive definite, symmetric matrix $\mathrm{W}\in\mathbb{R}^{\mathrm{N}\times\mathrm{N}}$. The goal is to minimize the mean projection error of the given snapshots projected onto the subspace spanned by the POD basis $\psi_{i}$. By introducing the matrix $\mathbf{A}_{\mathcal{K}}$ as the collection of the snapshots $\mathbf{a}^{k}$ with $k\in\mathcal{K}$ we can write the operator $\mathbf{R}$ arising from the optimization problem (5) as

[TABLE]

Then the unique solution to (5) is given by the eigenvectors corresponding to the $\mathcal{N}$ largest eigenvalues of $\mathbf{R}$, i.e., $\mathbf{R}\psi_{i}=\lambda_{i}\psi_{i}$ with $\lambda_{i}>0$ GV17 . The operator $\mathbf{R}$ is largesince it is of dimension $\mathrm{N}$ which we want to reduce. Hence it might be better in many cases to set up and solve the eigenvalue problem

[TABLE]

and obtain the POD basis by $\psi_{i}=1/\sqrt{\lambda_{i}}\mathbf{A}_{\mathcal{K}}v_{i}$. Note that both approach are equivalent and are related by the singular value decomposition (SVD) of the matrix $\mathrm{W}^{1/2}\mathbf{A}_{\mathcal{K}}$. While the latter is computationally more efficient, the singular value decomposition is numerically more stable. A comparison of the different computations was carried out in LV13 . For completeness let us state the POD approximation error given by

[TABLE]

where $d$ is the rank of $\mathbf{A}_{\mathcal{K}}$. For easier notation we collect the POD basis $\psi_{i}$ in the matrix $\Psi=[\psi_{1},\ldots,\psi_{\mathcal{N}}]\in\mathbb{R}^{\mathrm{N}\times\mathcal{N}}$.

After introducing the computation of the POD basis we will now outline the adaptive approach utilized in this work. It is crucial to minimize the number of solves involving the FEM discretization in order to obtain a speedup of the computation. The goal is to push most of the computations in the simulation onto the reduced order models. However, we have to guarantee that the reduced order models are accurate in order to obtain reliable results. We will present a strategy that does not require precomputation as for example in the reduced basis method but performs the model order reduction during the simulation. First let us outline how the POD basis is applied to (4). In the second step we give the details on how to obtain the basis efficiently.

We generate for each part of the machine an individual POD basis. Hence we have one basis for the stator and one basis for the rotor. The interface between the stator and the rotor is not reduced, i.e., we work with the FEM ansatz space on the interface. This is motivated by the observation that the decay of the eigenvalues on $\Gamma_{\mathrm{I}}$ is very slow, which would result in a large POD basis. Since the FEM space for the interface is usually of moderate dimension, the gain of using POD would be negligible. We compute the POD basis as solution to (5) utilizing the snapshots $\mathbf{a}_{s}^{k}$ and $\mathbf{a}_{r}^{k}$ to obtain $\Psi^{s}$ and $\Psi^{r}$, respectively. We then make the ansatz

[TABLE]

where the POD coefficients are indicated with a bar. By projecting (4) onto the subspace spanned by the POD basis, we obtain the reduced order model

[TABLE]

In short notation the system will be written as

[TABLE]

with $\bar{\mathbf{a}}$ the vector of POD coefficients. This system is of dimension $\mathcal{N}_{s}+\mathcal{N}_{r}+\mathrm{N}_{I}$ and of much smaller dimension as the original system (4) which was of dimension $\mathrm{N}$.

Next we introduce the strategy on how to determine the POD basis adaptively. The goal is to reduce the computational cost with respect to the rotation. A full revolution requires $\mathrm{N}_{I}$ solves of the system (4), i.e., for all $\vartheta_{k}$ with $k\in\mathcal{K}:=\{0,1,\ldots,\mathrm{N}_{I}-1\}$. In the symmetric case it is not required to solve a full rotation but only for angles that cover one pole, for our particular example this means one sixth, i.e., $\mathrm{N}_{I}/6$ solutions are needed. Note, we assume that $\mathrm{N}_{I}$ is divisible by $\mathrm{N}_{p}$, where $\mathrm{N}_{p}$ is the number of poles of the machine. In the non-symmetric case this is not possible. The idea is to generate a sequence of disjoint sets $\mathcal{K}_{i}$ with $i=0,\ldots,\mathrm{N}_{\mathcal{K}}$. Note that the sets can also be chosen arbitrarily as long as they fulfill the subset property. This is required to be able to reuse the already computed snapshots and hence minimize computational overhead. Ideally the sets $\mathcal{K}_{i}$ are not too large but are large enough to cover the most important features. The strategy is then as follows: We start by choosing $\mathcal{K}_{0}$ and evaluate (4) for $\vartheta_{k}$ and $k\in\mathcal{K}_{0}$. From the computed solutions/snapshots $A_{\mathcal{K}}=[\mathbf{a}^{k}]_{k\in\mathcal{K}_{0}}$ a POD basis is computed. Then an error estimator $\Delta\mathbf{a}(\vartheta_{k})$, $k=1,\ldots,N_{I}$, is evaluated to determine the maximum error. If the error is larger than a given tolerance the index $k\in\{1,\ldots,\mathrm{N}_{I}\}$ is determined where the maximum error is attained. Determine the set $\mathcal{K}_{i}$ that contains the index $k$ and the snapshot set is enlarged by adding the new solutions corresponding to $\mathcal{K}_{i}$ to the old ones, i.e., $A_{\mathcal{K}}=[A_{\mathcal{K}},[\mathbf{a}^{k}]_{\mathcal{K}_{i}}]$. This procedure is repeated until the error $\Delta\mathbf{a}(\vartheta_{k})$, $k=1,\ldots,N_{I}$, is below the desired threshold. For stability reasons the sets $\mathcal{K}_{i}$ that have been already used are removed from the list. In this case the index for the second largest error is used. During the numerical experiments this scenario never occurred and hence will not be investigated further. In Algorithm 1 the strategy is summarized. This sampling of the sets $\mathcal{K}_{i}$ is similar to the greedy algorithm from the reduced basis method. The decision to add more than one solution $\mathbf{a}^{k}$ at a time to the snapshot set is to minimize the overhead of evaluating the error estimator and generating the reduced order models. Since we do not introduce an online-offline decomposition the computations of error estimator and the generation of the reduced order models are included in the computational costs. In our numerical tests the proposed strategy converges very fast, in at most $6$ iterations. Depending on the asymmetry different sets $\mathcal{K}_{i}$ may preform better. This will be illustrated with numerical experiments. For the symmetric case we observe that usually only one set is visited (if chosen appropriately). In other words: the symmetry is correctly detected by the algorithm. In Section 4 the specific choices of the sets $\mathcal{K}_{i}$ are described and a comparison of different strategies is performed. The dimension $\mathcal{N}$ of the computation of the POD basis is chosen such that

[TABLE]

holds for the stator and rotor independently. This is a popular choice, where a typical value is $\varepsilon_{rel}=0.9999$. Note that the denominator can be computed by $\mathrm{trace}(\mathbf{A}_{\mathcal{K}}^{\top}\mathrm{W}\mathbf{A}_{\mathcal{K}})$ and hence not all $d$ eigenvalues have to be computed.

Let us now shortly have a look at the error estimator. For this let us recall some basic quantities. We define the discrete coercivity constant by

[TABLE]

Hence the the coercivity constant is given by the smallest eigenvalue such that

[TABLE]

is satisfied for $(\lambda,\mathbf{v})\in\mathbb{R}_{+}\times\mathbb{R}^{\mathrm{N}}$ and $\mathbf{v}\neq 0$ QMN16 . Further, we define the residual $r(\mathbf{a}^{\mathcal{N}};\vartheta)=\mathbf{f}(\vartheta)-\mathbf{K}_{\nu}(\vartheta)\mathbf{a}^{\mathcal{N}}$. Then the error introduced by the reduced order model in the variable $\mathbf{a}$ can be characterized by

[TABLE]

Additionally, we look at the relative error which might be more interesting in many applications. The corresponding error estimator reads

[TABLE]

This is a standard result and can be found in PR06 ; RHP08 . In the numerical realization we set the weight matrix $\mathbf{W}$ to the $\mathrm{N}$ dimensional identity matrix. Other choices (e.g., $\mathbf{W}=\mathbf{K}_{\nu}$) are possible but not investigated at this point. Note that the rotation as introduced in this work does not influence the coercivity constant and hence the dependence can be omitted. This can also be seen in the numerics, where only small deviations can be observed which are in the order of discretization. Hence a very efficient realization is possible, since only one eigenvalue problem has to be solved. Let us remark that the error is measured with respect to the finite element solution. It is assumed that the finite element solution is accurate enough to approximate the solution of the continuous problem. Let us remark that the computation of the residual norms can be performed very efficiently using the introduced affine decomposition DH15 .

4 Numerical results

We will now present different numerical results. For this let us specify the settings. The geometry (Figure 1) is discretized using a triangular mesh with $56297$ nodes and the interface $\Gamma_{\mathrm{I}}$ is discretized by $900$ equidistant points. Hence one revolution requires the linear system (3) to be solved $900$ times. Here is where model order reduction will come into play and speedup the simulation significantly. All computation are performed on a standard desktop PC using Matlab R2016b. Throughout our numerical tests we will consider three settings:

sym     

Symmetric machine.

rot      

Perturbation of $\phi$ in one permanent magnet by $5^{\circ}$.

stat     

Perturbation of $\ell$ in one tooth by $0.3$mm.

rot_stat

Perturbation in both $\phi$ and $\ell$.

To start we will have a look at the performance of the POD method. For this we compute a full revolution for each of the four setting and have a look at the decay of the eigenvalues. Let us recall that a fast decay is essential for the POD method to perform well. In Figure 2 the decay of the normalized eigenvalues is shown. As can be seen the eigenvalues decay very fast for rotor and stator. Only for the interface the decay is much slower. This underlines the decision to not perform a model order reduction for the interface. Additionally it can be observed that a perturbation in the magnetic field angle $\phi$ has less influence on the decay of the eigenvalues than the perturbation in the length $\ell$ of one stator tooth. While the perturbation in $\phi$ causes a slight change in the decay of the eigenvalues related to the stator, the perturbation in $\ell$ dramatically influences the eigenvalues related to the rotor.

Next we depict the first three POD basis vectors for each setting (Figure 3-6). For the setting sym and rot very similar basis vectors are obtained. The last basis vector for both settings is zero on wide areas and only adds contributions at the interface (green is zero). In the settings stat and rot_stat this is very different. There the third basis vector for the rotor still contributes a lot of information to the reduced order model. What was already observed in the decay of the eigenvalues was verified again in the POD basis vectors. The settings sym and rot exhibit similar behaviors as well as stat and rot_stat.

From the decay of the eigenvalues it can be expected that we will be able to generate reduced order models of very low dimension. The plots of POD basis verify this for the sym and rot since already the third basis vector is almost zero.

Next we have a look at the performance of the presented adaptive algorithm. For this let us introduce the index sets $\mathcal{K}_{i}$. We will look at two settings. To be able to differentiate the two settings we will call sets $\mathcal{M}$ and $\mathcal{K}$. In the first we choose the sets by dividing the full rotation into six sections (6 poles). within each section we have a hierarchy of sets for refinement. Hence the sets can be written as follows

[TABLE]

Note that the values for $\mathcal{K}_{ij}$ are always less than $150i$. Using these sets we always compute snapshots only related to one pole. If the error is still to large in the particular pole a refinement is achieved.

In the second setting distributed index sets are considered. This allows for a broader information collection since the snapshots are immediately distributed over the whole range of $\vartheta$. We define the corresponding sets by

[TABLE]

With these particular choices we get in both settings $72$ sets with each $12$ indices and hence we can perform a fair comparison of the two approaches. In our experiment it turned out that these settings are a good trade between performance and accuracy. As the tolerance for our adaptive strategy we use $\varepsilon=10^{-3}$ which is sufficient for most applications. The error estimator overestimated the real error by at most one order of magnitude. This is a good result since it guarantees that the model is not refined unnecessarily often.

The results for the two different index sets $\mathcal{K}$ and $\mathcal{M}$ are shown in Figure 7-8. We plot the relative error $\Delta^{\mathrm{rel}}_{\mathbf{a}}$ of the $900$ angular position for every iteration of the adaptive algorithm. The actual error is omitted in the plots for visual clarity. By a dashed line the tolerance is indicated as a visual aide. It can be seen how each of the two index sets have very different behaviors.

The set $\mathcal{K}$, which uses a kind of local refinement of the snapshots, performs well for the sym and rot. In the behaviour of the error it can be clearly seen in which region new snapshots were added. The error of the region drops significantly. As was expected, the symmetric case only requires one iteration since already a few snapshots contain enough information to compute the full rotation. For the perturbation in the rotor a second set of snapshots is required to push the error below the tolerance $\varepsilon$. A clear accuracy difference between the poles can be observed which is reflected by steps in the error. The settings stat and rot_stat are not handled too well by the local nature of the set $\mathcal{K}$. For each pole a snapshot set is selected resulting in six iterations for both settings.

On the other hand the index set $\mathcal{M}$ shows the advantages of the global nature of the snapshot selection. While for sym and rot this results in more iterations as for the set $\mathcal{K}$ in the other two settings, the benefits can be observed. In particular for stat, a much faster convergence of the adaptive algorithm is obtained. The error also has a different nature. In each iteration the error is reduced uniformly over the whole rotation and not only locally.

Next we have a look at the performance of the adaptive algorithm. For this we compare computational time (wall clock time) of the different approaches. The computational time is determined by the average over $10$ runs to flatten out irregularities in the numerical realization. Note that the mesh and the constant finite element matrices are precomputed, and hence are not reflected in the run time. This does not influence the result since it is required by both approaches and gives a more clear indication of the computational costs, which is the main focus. In Table 4-4 the results are summarized.