---
author:
- 
title: 'Shape Optimization of Rotating Electric Machines using Isogeometric Analysis and Harmonic Stator-Rotor Coupling'
---

Introduction
============

As a result of the energy transition, the simulation of electromechanical energy converters, in particular electric machines, is becoming increasingly important to obtain efficient and robust designs. Commonly, a workflow based on analytical estimates and the Finite Element Method (FEM) is used for 2D and finally 3D domains. In most classical approaches, the geometry is only approximated, e.g., with an accuracy that depends on the mesh refinement. These errors can be avoided when using Isogeometric Analysis (IGA) [@Hughes_2005aa], [@Buffa_2010aa], [@Bontinck_2017ag], [@Dolz_2018aa] which uses B-splines and/or Non-Uniform Rational B-splines (NURBS) as basis for geometry and solution space. In IGA, the geometry can be easily and smoothly transformed by moving the control points of the splines such that there is no need to remesh the domain when the geometry is modified. This makes IGA very well suited for shape optimization [@1809.03377]. Numerical optimization based on magnetic equivalent circuits or finite element models has lead to large improvements in the designs of technical applications, e.g., of permanent magnet synchronous machines. During the last 30 years there has been a lot of research on finite element based optimization methods (see, e.g., [@Di-Barba_2010aa], [@Duan_2013aa] and the references therein). Originally, mainly gradient-based optimization methods were used (see, e.g., [@Russenschuck_1990aa; @Weeber_1992aa; @Takorabet_1997aa]) but stochastic optimization algorithms became more popular during the last 20 years (see, e.g., [@Hameyer_1993aa], [@Lok_2017aa]). Most of the proposed algorithms use stochastic or population-based optimization, e.g., genetic algorithms and particle swarm optimization (see, e.g., [@Ma_2015aa]) which have also been extended to multi-objective optimization problems (see, e.g., [@Baumgartner_2004aa; @Ho_2005aa]). For permanent magnet synchronous machines, stochastic optimization methods are commonly used (see, e.g.,[@Cassimere_2009aa; @Bash_2011aa; @Sizov_2013aa]).

This contribution deals for the first time with shape optimization of a rotating electrical machine discretized with IGA. We use shape calculus to obtain shape derivatives such that gradient-based optimization becomes feasible. The rotor and stator domains are discretized separately and are coupled using harmonic stator-rotor coupling [@Bontinck_2018ac].

The paper is structured as follows: In Section \[sec:model\] the model of the electric machine and its quantities of interest are introduced. We introduce the shape optimization problem and derive the formula of the shape derivative using shape calculus in Section \[sec:opti\]. In Section \[sec:discretization\] the discretization of the machine using harmonic stator-rotor coupling and isogeometric analysis is explained and the resulting saddle-point problem is presented. We conclude the paper by explaining our gradient-based shape optimization algorithm and present numerical results for the minimization of the total harmonic distortion of the electromotive force in Section \[sec:num\_shape\_opti\].

Model of the Electric Machine {#sec:model}
=============================

Electromagnetic phenomena are described by Maxwell’s equations. For many applications it is sufficient to consider the magnetostatic approximation in a working point, i.e., neglecting displacement and eddy currents and nonlinearity of the materials. In a domain $D_{\mathrm{3D}}$ the magnetostatic formulation of Maxwell’s equations is given by $$\begin{aligned}
\nabla\times\left(\nu\nabla \times \mathbf{A}\right) = \mathbf{J}_{\mathrm{src}} + \nabla \times \mathbf{M} , \label{eq:curlcurl}\end{aligned}$$ with the piecewise constant reluctivity $\nu$, the magnetic vector potential $\mathbf{A}$, the current density $\mathbf{J}_{\mathrm{src}}=\sum_k\pmb{\chi}_ki_k$ given by winding functions $\pmb{\chi}_k$ and currents $i_k$ [@Schops_2013aa], the magnetization of the permanent magnets $\mathbf{M}$ and homogeneous Dirichlet boundary conditions $\mathbf{A}\times\mathbf{n}=0$ on $\partial D_{\mathrm{3D}}$, where $\mathbf{n}$ is the normal vector. The source current density $\mathbf{J}_{\mathrm{src}}$ and the permanent magnetization $\mathbf{M}$ vanish outside the coil (${D}_{\mathrm{c}}$) and permanent magnet (${D}_{\mathrm{pm}}$) regions, respectively, see Fig. \[fig:kartoffel\].

plot \[smooth cycle, tension=2\] coordinates [(2,1.5) (2,2.5) (1,2.7)]{}; plot \[smooth cycle, tension=2\] coordinates [(2,1.5) (2,2.5) (1,2.7)]{};

plot \[smooth cycle, tension=2\] coordinates [(3,3) (3.5,3.5) (4,2.5)]{}; plot \[smooth cycle, tension=2\] coordinates [(3,3) (3.5,3.5) (4,2.5)]{};

plot \[smooth cycle, tension=2\] coordinates [(6,3) (6,4) (4.5,2.5)]{}; plot \[smooth cycle, tension=2\] coordinates [(6,3) (6,4) (4.5,2.5)]{};

plot \[smooth cycle, tension=2, pattern color=blue, pattern=north east lines\] coordinates [(0,1) (4,4.5) (6,2)]{}; plot \[smooth cycle, tension=2\] coordinates [(0,1) (4,4.5) (6,2)]{};

at (3.5,1.5) ; (dc) at (2.7,2.5) ; at (5.5,3.2) ; at (1.3,4) [$\partial{D}$]{}; (dc) – (3.5,3); (dc) – (1.8,2.2);

![Multipatch model of one pole of a 6-pole permanent magnet synchronous machine. The magnetic flux density $\mathbf{B}$ is shown which was computed using isogeometric analysis.[]{data-label="fig:pmsm"}](images/B_labels.pdf){width="35.00000%"}

A common quantity of interest when designing electric machines is the total harmonic distortion ($\operatorname{THD}$) of the electromotive force (EMF) ${\mathcal E}$. The total harmonic distortion of a function $f=f(t)$ is defined as

$$\begin{aligned}
\operatorname{THD}_{\operatorname{\mathcal{I}}}(f) = \frac{\sqrt{\sum_{n\in\operatorname{\mathcal{I}},n\neq1}|c_n|^2}}{|c_1|},\end{aligned}$$

where $c_n$ are the coefficients of the Fourier series of $f(t)$, i.e. $$\label{expFourier1}
    f(t) = \sum_{n=-\infty}^\infty c_n e^{\imath n t},$$ and $\operatorname{\mathcal{I}}\subset\operatorname{\mathbb{N}}$ is an index set of frequencies to consider, e.g., one may disregard frequencies that cannot be diminished by shape optimization. Alternatively, the expansion can be written as $$\begin{aligned}
    f(t) = \frac{A_0}{2} + \sum_{n=1}^\infty A_n \cos (nt) + B_n \sin(nt)\end{aligned}$$ with coefficients $A_n = c_n + c_{-n}$ for $n \geq 0$ and $ B_n = \imath(c_n - c_{-n})$ for $n \geq 1$. Note that it holds that $|c_n| = \sqrt{A_n^2+B_n^2}/2$. Using these coefficients, the THD can be written as $$\begin{aligned}
 \label{eq_thd_AB}
    \operatorname{THD}(f)=\sqrt{ \frac{\sum_{n\in \operatorname{\mathcal{I}}, n\neq 1} A_n^2+B_n^2}{A_1^2+B_1^2} }.\end{aligned}$$

We are interested in the THD of the EMF. Therefore, the amplitudes $c_n$ are obtained by Fourier analysis of the voltages $$\label{eq:voltages} 
\mathcal{E}_k(\mathbf{A}(t)) =\partial_t{\Psi}_k(\mathbf{A}(t)),$$ with the flux linkage $$\label{eq:fluxlinkage} 
{\Psi}_k(\mathbf{A}(t)) =N_p\int_{D_{\mathrm{3D}}}\pmb{\chi}_k \cdot \mathbf{A}(t)\; {\mathrm{d}} D_{\mathrm{3D}},$$ where $N_p$ is the number of pole pairs. The voltages $\mathcal{E}_k$ are induced by the electromagnetic force in the windings under the assumption that the coils of the poles are connected in series.

In the planar 2D case is considered on a cross section ${D} \subset \mathbb R^2$ and reduces to $$\begin{aligned}
-\nabla\cdot\left(\nu\nabla u\right) = J_{\mathrm{src},z} + \nabla \times \mathbf{M} \cdot  \mathbf{e}_z, \label{eq:reducedmaxwell}\end{aligned}$$ where $u(x,y)$ is the $z$-component of the magnetic vector potential $\mathbf{A}=[0,0,u]^\top$, $J_{\mathrm{src},z}=\sum_k\chi_ki_k$ is the $z$-component of the current density given by winding functions $\chi_k$ and currents $i_k$ and $\mathbf{e}_z$ is the unit vector in $z$-direction.

Optimization {#sec:opti}
============

Problem statement
-----------------

We want to optimize the $2D$ shape $\Omega$ of the $3$-phase permanent magnet synchronous machine from Fig. \[fig:pmsm\], which has $N_p=3$ pole pairs, in terms of the total harmonic distortion ($\operatorname{THD}$) of the electromotive force (EMF) as introduced in . Due to periodicity of the stator windings, we will restrict ourselves to the optimization of the first voltage and set $\mathcal E := \mathcal E_1$ and $\chi = \chi_1$. Note the dependence of $\mathcal E$ on the shape $\Omega$ via the solution $u$ to , i.e., $\mathcal E = \mathcal E(u(t,\Omega))$. Thus, we consider the optimization problem $$\min_{\Omega\in\mathcal{A}} \mathcal{J}(\Omega) := \operatorname{THD}_{\operatorname{\mathcal{I}}}({\mathcal E}(u(t,\Omega))), \label{eq_minTHD}$$ that is constrained by the machine model , which, for $t \in  [0,T]$, reads in its weak form: Find $u = u(\alpha(t))$ such that $$\label{pde_constraint_weak}
    \underbrace{
    \int_{{D}} \nu_{\Omega}(\alpha(t))\nabla u \cdot\nabla v\,{\mathrm{d}} {D}}
    _{=:a[\alpha(t)](\Omega;u,v)}
     = 
     \left\langle r(t, \alpha(t)), v \right\rangle$$ for all test functions $v$. Here $T=2\pi/\omega$ is the (electrical) period length, $\alpha(t)$ the rotor angle, which may be given by the equation of motion. We allow all geometries from an admissible set $\mathcal{A}$. Both $u$ and $v$ satisfy homogeneous Dirichlet boundary conditions on the inner and outer circular parts of the boundary and periodic boundary conditions on the left and right parts. Finally, the right-hand-side is given by $$\begin{aligned}
\begin{aligned} 
 \langle r(t, &\alpha(t)), v \rangle = \int_{D_{\mathrm{c}}}\sum_{k=1}^{3}\chi_k(x,y)i_k(t) \, v(x,y) \; {\mathrm{d}} D_{\mathrm{c}} \\
 &+ \int_{D_{\mathrm{pm}}(\alpha(t))} \binom{-M_2}{M_1} \cdot \nabla v(x,y) \; {\mathrm{d}} D_{\mathrm{pm}}(\alpha(t)) \label{eq:rhs}
\end{aligned}\end{aligned}$$ and contains the excitations due to permanent magnets and the coils. Here, $M_1$ and $M_2$ are the first and second component of the magnetization vector $\mathbf{M}$, respectively, $D_{\mathrm{pm}}(\alpha(t))$ denotes the permanent magnet region after rotation by the angle $\alpha(t)$ and we exploit for a more compact notation that the permanent magnets lie in the rotating part and the coils in the stator. We indicate the dependence of the reluctivity function on the current shape $\Omega$ which is subject to the shape optimization by writing $\nu_{\Omega}$.

Moreover, we introduce a temporal discretization $\lbrace t_1, \dots, t_{N_\alpha} \rbrace$ into ${N_{\alpha}}=120$ points in time and a corresponding discretization of the range of angular displacements into ${N_{\alpha}}$ rotor positions, $\boldsymbol \alpha := (\alpha_1,\dots,\alpha_{{N_{\alpha}}})$. For $j \in \{ 1,\dots, {N_{\alpha}} \}$ and a given shape $\Omega$, let $a_j(\Omega;\cdot, \cdot) := a[\alpha_j](\Omega;\cdot,\cdot)$ and $r_j:=r(t_j,\alpha_j)$ according to the definitions in and , respectively. Approximating $u(t,\Omega)$ for $t\in [0,T]$ by $\lbrace u_1, \dots, u_{N_\alpha} \rbrace$ and introducing the notation $\mathcal{J}(u_1, \dots, u_{N_\alpha}) := \operatorname{THD}_{\operatorname{\mathcal{I}}}({\mathcal E}(u_1, \dots, u_{N_\alpha}))$, problem – after discretization with respect to rotor positions can be written as the optimization problem $$\qquad\underset{\Omega}{\mbox{min }} \mathcal{J}(u_1, \dots, u_{{N_{\alpha}}}) \label{minJ}$$ $$\begin{aligned}
 \label{pdeconstraints}
    \mbox{s.t. } \left\lbrace \begin{array}{c}
                    \;\quad a_1(\Omega;u_1,v_1) = \langle r_1, v_1 \rangle \; \forall v_1,\\
                    \;\;\qquad\vdots\\
                    a_{{N_{\alpha}}}(\Omega;u_{{N_{\alpha}}},v_{{N_{\alpha}}}) = \langle r_{{N_{\alpha}}}, v_{{N_{\alpha}}} \rangle  \; \forall  v_{{N_{\alpha}}},
                              \end{array} \right. \end{aligned}$$ which is constrained by ${N_{\alpha}}$ boundary value problems corresponding to the ${N_{\alpha}}$ rotor positions under consideration.

Shape Sensitivity Analysis
--------------------------

The shape derivative $d \mathcal{J}(\Omega; \mathbf{W})$ of a domain-dependent functional $\mathcal{J}= \mathcal{J}(\Omega)$ represents the sensitivity of the functional with respect to a perturbation of the domain in the direction of a given vector field $\mathbf{W}$. The shape derivative is defined as $$d \mathcal{J}(\Omega; \mathbf{W}) = \underset{\delta \searrow 0}{\mbox{lim}} \frac{\mathcal{J}(T_\delta^{\mathbf{W}}(\Omega)) - \mathcal{J}(\Omega)}{\delta},$$ if this limit exists and the mapping $\mathbf{W} \mapsto d \mathcal{J}(\Omega; \mathbf{W})$ is linear and continuous on the space of smooth vector fields. Here, $T_\delta^{\mathbf{W}}$ represents a transformation which moves every point $\operatorname{\mathbf{x}}$ a distance $\delta>0$ in the direction given by the vector field $\mathbf{W}$, $T_\delta^{\mathbf{W}}(\operatorname{\mathbf{x}}) = \operatorname{\mathbf{x}}+ \delta \mathbf{W}(\operatorname{\mathbf{x}})$.

For deriving the shape derivative $d\mathcal{J}(\Omega; \mathbf{W})$ for the optimization problem –, we follow the steps taken in [@GanglLangerLaurainMeftahiSturm2015]. First, we introduce the Lagrangian $$\begin{aligned}
    \mathcal L(\Omega, &\varphi_1,\dots,\varphi_{{N_{\alpha}}},\psi_1,\dots, \psi_{{N_{\alpha}}}) := \\
    & \mathcal{J}(\varphi_1, \dots, \varphi_{{N_{\alpha}}}) + \sum_{k=1}^{{N_{\alpha}}} \left(a_k(\varphi_k,\psi_k) - \langle r_k, \psi_k \rangle \right).\end{aligned}$$ For ease of notation, we will use the notation $\mathbf{u} = (u_1,\dots, u_{{N_{\alpha}}})$ and similar for other quantities. Note that, for $u_i$ satisfying the $i$-th equation of , it holds that $$\begin{aligned}
    \frac{\partial}{\partial \psi_i} \mathcal L(\Omega, \mathbf{u},\boldsymbol \psi)(q_i) 
     = a_i(u_i,q_i) - \langle r_i,q_i\rangle = 0\end{aligned}$$ for any test function $q_i$. Similarly, we introduce the adjoint states $p_i$, $i=1,\dots,{N_{\alpha}}$, as the solutions to $$\begin{aligned}
 \label{eq_adjoint}
    0= \frac{\partial}{\partial \varphi_i} \mathcal L(\Omega, \mathbf{u},\mathbf{p})(v_i) =  \frac{\partial \mathcal{J}}{\partial \varphi_i}(\mathbf{u})(v_i) + a_i(v_i,p_i)\end{aligned}$$ for all test functions $v_i$.

Let us now consider the functional $\mathcal{J}= \mathcal{J}(u_1, \dots, u_{{N_{\alpha}}})$ more closely. We will use the representation for the THD. [Let $(C_n)_{n \in \operatorname{\mathcal{I}}}$ the Fourier coefficients of the function ${\Psi}(u_1, \dots, u_{{N_{\alpha}}})$, i.e., $$\label{eq:flux_fourier}
    {\Psi}(u_1, \dots, u_{{N_{\alpha}}}) =  \sum_{n=-{N_{\alpha}}}^{N_{\alpha}} C_n e^{\imath n t}.$$ Using the Fourier representation of the flux linkage and , the electromotive force can be written as $${\mathcal E}(u_1, \dots, u_{{N_{\alpha}}}) =  \sum_{n=-{N_{\alpha}}}^{N_{\alpha}} \underbrace{\imath n C_n}_{c_n} e^{\imath n t},$$ where $(c_n)_{n \in \operatorname{\mathcal{I}}}$ are the Fourier coefficients of the electromotive force ${\mathcal E}$. This Fourier series can be rewritten as $${\mathcal E}(u_1, \dots, u_{{N_{\alpha}}}) = \frac{A_0}{2} + \sum_{n=1}^{N_{\alpha}} A_n \cos (nt) + B_n \sin(nt). \label{eq:fourierSeries}$$ The coefficients $A_n$ and $B_n$ of the Fourier representation of $\partial_{t}{\Psi}$ are obtained by time derivation and Fourier transform, which we will denote by $A_n=[\mathcal{F'}_a({\Psi}(u_1, \dots, u_{{N_{\alpha}}}))]_n$ and $B_n=[\mathcal{F'}_b({\Psi}(u_1, \dots, u_{{N_{\alpha}}}))]_n$, respectively. Exploiting the linearity of the discrete Fourier transform, the vectors $\mathbf{A}$ and $\mathbf{B}$ can also be written in terms of transformation matrices $\mathbf{M}_a = \mathcal{F'}_a({\mathbf{I}})$ and $\mathbf{M}_b=\mathcal{F'}_b({\mathbf{I}})$, where ${\mathbf{I}}$ is the identity matrix, i.e., $$\begin{aligned}
    A_k(u_1,\dots,u_{{N_{\alpha}}}) &= \sum_{j=1}^{N_{\alpha}}\left( \left(\mathbf{M}_a\right)_{k,j} N_p l_z\int_{{D}}\chi u_j\; {\mathrm{d}} {D} \right),\label{eq:def_Ak}\\
    B_k(u_1,\dots,u_{{N_{\alpha}}}) &= \sum_{j=1}^{N_{\alpha}} \left( \left(\mathbf{M}_a\right)_{k,j} N_p l_z\int_{{D}}\chi u_j\; {\mathrm{d}} {D} \right),\label{eq:def_Bk}\end{aligned}$$ where $l_z$ is the length of the machine in $z$-direction. ]{}

In order to solve the adjoint state equation , we need to differentiate the functional $\mathcal{J} = \operatorname{THD}({\mathcal E}(u_1,\dots,u_{{N_{\alpha}}}))$ with respect to $u_i$ for $i \in \{1,\dots,{N_{\alpha}}\}$. Using the relation , we get $$\begin{aligned}
    &\frac{d\, \operatorname{THD}({\mathcal E}(\mathbf{u}))}{d u_i}(\mathbf{u})(v_i) = \\
    &\frac{1}{\sqrt{A_1^2+B_1^2}}  \frac{1}{ \sqrt{\sum_{k\in \operatorname{\mathcal{I}}, k\neq 1} A_k^2 + B_k^2} } \left(\sum_{k\in \operatorname{\mathcal{I}}, k\neq 1} A_k A_k' + B_k B_k'\right) \\
    &- \frac{1}{(A_1^2+B_1^2)^{3/2}} \left( \sqrt{\sum_{k\in \operatorname{\mathcal{I}}, k\neq 1} A_k^2+B_k^2} \left( A_1 A_1' + B_1 B_1'\right) \right).\end{aligned}$$ Here, we used the abbreviations $A_k' := \frac{d A_k}{d u_i}(\mathbf{u})(v_i)$ and $B_k' := \frac{d B_k}{d u_i}(\mathbf{u})(v_i)$. It can be seen from and that $$\begin{aligned}
    \frac{d A_k}{d u_i}(\mathbf{u})(v_i) &= \left(\mathbf{M}_a\right)_{k,i} N_p l_z \int_D \chi v_i \; {\mathrm{d}} {D}, \\
    \frac{d B_k}{d u_i}(\mathbf{u})(v_i) &= \left(\mathbf{M}_b\right)_{k,i} N_p l_z \int_D \chi v_i \; {\mathrm{d}} {D}.\end{aligned}$$ Given the solutions to the forward problem for all rotor positions $l \in \{1,\dots {N_{\alpha}}\}$, we obtain for the adjoint problem defining the adjoint variable at rotor position $l \in \{1, \dots {N_{\alpha}}\}$: Find $p_i$ such that $$\begin{aligned}
    a_i(v_i, p_i) = - \frac{d\, \operatorname{THD}({\mathcal E}(\mathbf{u}))}{d u_i}(\mathbf{u})(v_i) \end{aligned}$$ for all test functions $v_i$.

Finally, in a similar way as it was proposed in [@GanglLangerLaurainMeftahiSturm2015], assuming that the deformation vector field $\mathbf{W}$ vanishes on the interface $\Gamma$ and on the stator, and that the permanent magnet region remains unchanged, we obtain the formula for the shape derivative in the direction of a smooth vector field $\mathbf{W} \in C^1(D,\mathbb R^2)$: $$\begin{aligned}
 \label{eq_shape_der}
    \begin{aligned}
        d\mathcal{J}&(\Omega; \mathbf{W}) = \sum_{l=1}^{N_{\alpha}} \int_{\Omega_{\mathrm{pm}}} (\nabla \cdot (\mathbf{W}) \mathbf{I} - \mathbf{\mathbf{DW}}^\top) \nabla p_l \cdot \binom{-M_2}{M_1} \;  {\mathrm{d}} x \\
        &+ \sum_{l=1}^{N_{\alpha}}\int_D \nu \, (\nabla \cdot \mathbf{W} \, \mathbf{I} - \mathbf{\mathbf{DW}}^\top - \mathbf{\mathbf{DW}}) \nabla u_l \cdot \nabla p_l \; {\mathrm{d}} x.       
    \end{aligned}\end{aligned}$$ Here, $\mathbf{I} \in \mathbb R^{2\times 2}$ denotes the two-dimensional identity matrix.

Solution of the Forward Problem by IGA {#sec:discretization}
======================================

In this section, we will treat the numerical solution of the forward problem using Isogeometric Analysis and harmonic stator-rotor coupling allowing for a flexible treatment of the rotation.

Harmonic Stator-Rotor Coupling
------------------------------

As proposed in [@Bontinck_2018ac], stator and rotor domains ${D}_q$, where $q\in \{ \text{st},\text{rt} \}$ distinguishes between stator and rotor and $\overline{ {D}}=  \overline {D}_{\mathrm{rt}} \cup \overline  {D}_{\mathrm{st}}$, are considered separately. They are coupled at the stator-rotor interface $\Gamma={\overline{{D}}_{\mathrm{rt}}} \cap {\overline{{D}}_{\mathrm{st}}}$ in the air gap by enforcing the continuity of the magnetic vector potential $u$ and of the azimuthal component of the magnetic field strength $H_{\theta}^{(\mathrm{st})}|_{\Gamma}(\theta_{\mathrm{st}})=H_{\theta}^{(\mathrm{rt})}|_{\Gamma}(\theta_{\mathrm{rt}})$, where $\theta_{\mathrm{st}}$ and $\theta_{\mathrm{rt}}$ are the angular coordinates attached to stator and rotor domain, respectively. The angular displacement between the domains is $\alpha=\theta_{\mathrm{st}}-\theta_{\mathrm{rt}}$. In its weak form the problem can be formulated as:

Find $({{u}_{\mathrm{st}}},{{u}_{\mathrm{rt}}},\lambda) \in {V_{\mathrm{st}}} \times {V_{\mathrm{rt}}} \times{\Lambda}$ such that $$\begin{aligned}
 \sum_q \int\limits_{{D}_{q}} \nu_{q}{\nabla u_{q}} \cdot {\nabla v_{q}} \operatorname{d} {D}_{q} + \int\limits_{\Gamma} \lambda \left\llbracket v\right\rrbracket \operatorname{d}\Gamma&= {\sum_{q} \langle r_q, v_q \rangle}\\
\intertext{with interface condition} \int_\Gamma \left\llbracket u\right\rrbracket \mu \operatorname{d}\Gamma &= 0, \end{aligned}$$ for all $v_q \in V_q$ and $\mu \in \Lambda$, $q\in \{ \mathrm{st},\mathrm{rt} \}$, where $V_q$ are suitable function spaces which contain the necessary Dirichlet boundary conditions, $v_q$ are the test functions, and $r_q$, $q \in \{ \mathrm{st}, \mathrm{rt} \}$ are the contributions of the right hand side on the stator and rotor as in , $$\begin{aligned}
    \langle r_{\mathrm{st}}, v_{\mathrm{st}} \rangle &= \int_{D_{\mathrm{c}}} \sum_k\chi_{k}(x,y)i_k(t) \, v_{\mathrm{st}}(x,y) \;{\mathrm{d}} D_{\mathrm{c}}, \\
    \langle r_{\mathrm{rt}}, v_{\mathrm{rt}} \rangle &= \int_{D_{\text{pm}}(\alpha)}  \binom{-M_2}{M_1} \cdot \nabla v(x,y) \;{\mathrm{d}} D_{\text{pm}}(\alpha) .\end{aligned}$$

Moreover, $\llbracket v\rrbracket =(v_{\textrm{st}}-v_{\textrm{rt}})\big |_{\Gamma}$ denotes the jump of $v$ across the interface $\Gamma$ and $\lambda = \nu{\nabla u^{(\mathrm{st})}}\cdot{\bm{\mathbf{e}}}_{r} = \nu{\nabla u^{(\mathrm{rt})}}\cdot{\bm{\mathbf{e}}}_{r} = H_\theta^{(q)}$, where ${\bm{\mathbf{e}}}_{r}$ is the unit vector in radial direction, is the $\theta$-component of the magnetic field at the interface and can be interpreted as a Lagrange multiplier [@Hansbo_2005aa].

Discretization
--------------

The $z$-component of the magnetic vector potential $u$ is discretized by a linear combination of the same scalar basis functions $w^{(q)}_{j}$ used for weighting and $H_{\theta}$ is expressed by a superposition of ${N_\Gamma}$ harmonic basis functions [@De-Gersem_2004ad] as $$\begin{aligned}
\label{eq:basis}
 u^{(q)}\approx \sum \limits_{j=1}^{N_q} u^{(q)}_{j} w^{(q)}_{j}, \quad
 H_{\theta}^{(q)}(\theta_q) \approx \sum_{k=1}^{{N_\Gamma}} \lambda_{k}^{(q)} e^{-\imath\ell_k\theta_q},\end{aligned}$$ where $\bm{\lambda}^{(q)}\in \mathbb{R}^{{N_\Gamma}}$ is the vector of Fourier coefficients for the harmonic basis functions, $\ell_k$ are the harmonic orders and $\mathbf{u}^{(q)}\in\mathbb{R}^{N_{q}}$ is the vector of degrees of freedom for each domain. This allows for a more efficient simulation of angular displacements $\alpha=\theta_{\mathrm{st}}-\theta_{\mathrm{rt}}$ when compared to classical mortaring, moving-band or sliding-surface methods [@De-Gersem_2004af]. Note that $u(\alpha)=u(x,y;\alpha)$ with $(x,y)\in{D}$, and the degrees of freedom $\mathbf{u}^{(q)}(\alpha)$ depend on the rotor position. The approach of harmonic stator-rotor coupling in an IGA framework leads to the Mortar-type saddle-point problem [@Bontinck_2018ac], $$\begin{bmatrix}
        {{\mathbf{K}}_{\mathrm{st}}}    &   0   &   -{{\mathbf{G}}_{\mathrm{st}}}\mathbf{R}(\alpha)\\
        0   &   {{\mathbf{K}}_{\mathrm{rt}}}    &   {{\mathbf{G}}_{\mathrm{rt}}}\\
-\mathbf{R}(\alpha){{\mathbf{G}}_{\mathrm{st}}^H}   &   {{\mathbf{G}}_{\mathrm{rt}}^H}  &   0 
\end{bmatrix}
\hspace{-0.1cm}
\begin{bmatrix}
{{\mathbf{u}}^{(\mathrm{st})}}\\
{{\mathbf{u}}^{(\mathrm{rt})}}\\
\bm{\lambda}
\end{bmatrix}=
\begin{bmatrix}
{{\mathbf{j}}_{\mathrm{st}}}\\
{{\mathbf{j}}_{\mathrm{rt}}}\\
0
\end{bmatrix}. \label{eq:system}$$ Here, $\bm{\lambda}=\bm{\lambda}^{(\mathrm{st})}=\bm{\lambda}^{(\mathrm{rt})}$ is enforced strongly by testing with $v=w$, $\mathbf{K}_q \in \mathbb{R}^{N_{q} \times N_{q}}$ are the stiffness matrices with $$(\mathbf K_q)_{ij} = \int_{{D}_{q}} \left(\nu{\frac{\partial w_i}{\partial x}}{\frac{\partial w_j}{\partial x}}+\nu{\frac{\partial w_i}{\partial y}}{\frac{\partial w_j}{\partial y}} \right)\;\text{d}{D}_{q},$$ $\mathbf{G}_q \in \mathbb{C}^{N_{q} \times {N_\Gamma}}$ are the coupling matrices with $$\begin{aligned}
(\mathbf{G}_{\mathrm{st}})_{ik} &=-\int_{0}^{2\pi}e^{-\imath\ell_k\theta_{\mathrm{st}}} w_i\left(\mathbf{r}\left(\theta_{\mathrm{st}}\right)\right)R_{\Gamma}\;{\mathrm{d}}\theta_{\mathrm{st}}, \\
(\mathbf{G}_{\mathrm{rt}})_{ik} &= \int_{0}^{2\pi}e^{-\imath\ell_k\theta_{\mathrm{rt}}}w_i\left(\mathbf{r}\left(\theta_{\mathrm{rt}}\right)\right)R_{\Gamma}\;{\mathrm{d}}\theta_{\mathrm{rt}},\end{aligned}$$ where $R_{\Gamma}$ is the radius of the interface $\Gamma$, $\mathbf{r}(\theta)=\left(R_{\Gamma}\cos(\theta), R_{\Gamma}\sin(\theta)\right)^{\top}$ is a mapping from an angle $\theta$ to the point on the interface at this angle in Cartesian coordinates and $\mathbf{R}(\alpha) \in \mathbb{C}^{{N_\Gamma} \times {N_\Gamma}}$ is the diagonal rotation matrix with $\left(\mathbf{R}(\alpha)\right)_{kk}=e^{\imath\ell_k\alpha}$. The saddle-point problem is stable if the numbers $N_q$ and ${N_\Gamma}$ of basis functions are chosen consistently [@Bontinck_2018ac].

Elimination of Inner Degrees of Freedom
---------------------------------------

The system can be rewritten as $$\begin{aligned}
\begin{bmatrix}
\mathbf{I} & 0 & -{{\mathbf{K}}_{\mathrm{st}}}^{-1}{{\mathbf{G}}_{\mathrm{st}}}\mathbf{R}(\alpha)\\
0 & \mathbf{I} & {{\mathbf{K}}_{\mathrm{rt}}}^{-1}{{\mathbf{G}}_{\mathrm{rt}}}\\
-\mathbf{R}(\alpha){{\mathbf{G}}_{\mathrm{st}}^H} & {{\mathbf{G}}_{\mathrm{rt}}}^H & 0\\
\end{bmatrix}
\begin{bmatrix}
{{\mathbf{u}}^{(\mathrm{st})}}\\
{{\mathbf{u}}^{(\mathrm{rt})}}\\
\bm{\lambda}
\end{bmatrix}
&=
\begin{bmatrix}
{{\mathbf{K}}_{\mathrm{st}}}^{-1}{{\mathbf{j}}_{\mathrm{st}}}\\{{\mathbf{K}}_{\mathrm{rt}}}^{-1}{{\mathbf{j}}_{\mathrm{rt}}}\\0
\end{bmatrix}.\end{aligned}$$ The internal degrees of freedom can be eliminated using the *Schur-complement*, giving rise to the low-dimensional *interface problem* $$\begin{aligned}
\mathbf{K}_{\textrm{int}}(\alpha)
\bm{\lambda}^{(\mathrm{rt})}
&=
\mathbf{f}_{\textrm{int}}(\alpha) \label{eq:sysint}\end{aligned}$$ with $$\begin{aligned}
\mathbf{K}_{\textrm{int}}(\alpha)&={{\mathbf{G}}_{\mathrm{rt}}^H} {{{\mathbf{K}}_{\mathrm{rt}}}^{-1}{{\mathbf{G}}_{\mathrm{rt}}}} + \mathbf{R}(\alpha){{\mathbf{G}}_{\mathrm{st}}^H} {{{\mathbf{K}}_{\mathrm{st}}}^{-1}{{\mathbf{G}}_{\mathrm{st}}}}\mathbf{R}(\alpha), \label{eq:kint}
    \\
    \mathbf{f}_{\textrm{int}}(\alpha)&={{\mathbf{G}}_{\mathrm{rt}}}^H {{{\mathbf{K}}_{\mathrm{rt}}}^{-1}{{\mathbf{j}}_{\mathrm{rt}}}} - \mathbf{R}(\alpha){{\mathbf{G}}_{\mathrm{st}}^H} {{{\mathbf{K}}_{\mathrm{st}}}^{-1}{{\mathbf{j}}_{\mathrm{st}}}}. \label{eq:fint} \end{aligned}$$ The inverses in and are not needed explicitly. Instead, one factorization and a few forward/backward substitutions (for each spectral basis at the interface) can be used to precompute the necessary expressions, e.g., $\mathbf{K}_{*}^{-1}\mathbf{G}_{*}$, which are independent of $\alpha$. Thus, only the small system with the system matrix $\mathbf{K}_{\textrm{int}} \in \mathbb{C}^{{N_\Gamma} \times {N_\Gamma}}$ has to be solved in the online phase for different rotation angles.

The internal degrees of freedom can be cheaply reconstructed, i.e., $$\begin{aligned}
    {{\mathbf{u}}^{(\mathrm{st})}}&={{{\mathbf{K}}_{\mathrm{st}}}^{-1}{{\mathbf{j}}_{\mathrm{st}}}}+{{{\mathbf{K}}_{\mathrm{st}}}^{-1}{{\mathbf{G}}_{\mathrm{st}}}}\mathbf{R}(\alpha)\bm{\lambda}^{(\mathrm{rt})},\label{eq:recostructionust}\\
    {{\mathbf{u}}^{(\mathrm{rt})}}&={{{\mathbf{K}}_{\mathrm{rt}}}^{-1}{{\mathbf{j}}_{\mathrm{rt}}}}-{{{\mathbf{K}}_{\mathrm{rt}}}^{-1}{{\mathbf{G}}_{\mathrm{rt}}}}\bm{\lambda}^{(\mathrm{rt})}.\label{eq:recostructionurt} \end{aligned}$$ This is computationally convenient when dealing with rotation since only the low-dimensional matrix $\mathbf{R}(\alpha)$ depends on $\alpha$. This leads to a significant reduction in computational cost of the shape optimization, e.g., during the computation of the total harmonic distortion, where has to be solved for $N_{\alpha}$ rotor positions. The computational cost of the simulation of a rotor rotation of $\SI{120}{\degree}$ with $N_{\alpha}=120$ discretized points in time using the interface problem is compared to the computational cost using the full system in \[fig:schur\_vs\_full\]. It can be seen that the computational time needed to solve the full system \[eq:system\] depends on the number of the IGA degrees of freedom whereas the computational time needed to solve the interface problem \[eq:sysint\] is almost independent of the number of IGA degrees of freedom and depends on the number of harmonics that are used for the coupling. The rotor is rotated $\SI{120}{\degree}$ using $N_{\alpha}=120$ time steps. The simulation is carried out for different spatial refinements and for different numbers of harmonics ${N_\Gamma}$ for the coupling. The preprocessing, i.e., computing the inverses in \[eq:kint,eq:fint\], took at most in the case of ${N_\Gamma}=50$ and $25746$ IGA DoFs and the postprocessing, i.e., \[eq:recostructionust,eq:recostructionurt\], took at most $\SI{1.3012}{s}$. The computation is is repeated 100 times and the computational time is averaged. The computation is carried out in Matlab^^ R2019a on a 6-core machine (Intel^^ Core i7-5820K CPU) with RAM.

table \[x index=0, y index=1, col sep=comma\] [images/dofs\_timeFullNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=2, col sep=comma\] [images/dofs\_timeFullNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=3, col sep=comma\] [images/dofs\_timeFullNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=4, col sep=comma\] [images/dofs\_timeFullNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=1, col sep=comma\] [images/dofs\_timeSchurNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=2, col sep=comma\] [images/dofs\_timeSchurNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=3, col sep=comma\] [images/dofs\_timeSchurNharm\_6\_20\_36\_50\_std.csv]{}; table \[x index=0, y index=4, col sep=comma\] [images/dofs\_timeSchurNharm\_6\_20\_36\_50\_std.csv]{};

B-splines and NURBS {#sec:splines}
-------------------

Isogeometric analysis is based on the idea of using freeform curves for an exact geometry representation of CAD models. This avoids errors in the geometry induced by mesh generation. For this, the notion of *B-splines*, i.e., piecewise polynomial functions, is of central importance. B-splines are defined via *knot vector* $\Xi=\left\{\xi_{1}, \xi_{2}, \dots, \xi_{n+p+1} \right\}$ where the knots $\xi_{i} \in [0,1]$ for all $i \in \{1,\dots,n+p+1\}$ are coordinates in the parametric space, $p$ is the degree of the B-spline and $n$ will be the dimension of the space. B-splines are defined by the recursive Cox-de Boor formula [@Piegl_1997aa] $$\begin{aligned}
 B_{i,p}(\xi) = \frac{\xi-\xi_{i}}{\xi_{i+p}-\xi_{i}}B_{i,p-1}(\xi)
              + \frac{\xi_{i+p+1}-\xi}{\xi_{i+p+1}-\xi_{i+1}}B_{i+1,p-1}(\xi), \end{aligned}$$ for all $p \geq 1$ and for $p=0$ via $$\begin{aligned}
 B_{i,0}(\xi) = \begin{cases}
            1 \quad \text{if } \xi_{i} \leq \xi < \xi_{i+1},\\
            0 \quad \text{else},
           \end{cases}\end{aligned}$$ where $0/0=0$ is formally assumed. *NURBS* (Non-Uniform Rational B-splines) basis functions are then defined by $$\begin{aligned}
 N_{i,p}(\xi) = \frac{w_{i} B_{i,p}(\xi)}{\sum_{k=1}^{n} w_{k} B_{k,p}(\xi)},\end{aligned}$$ where $w_{k}>0$ for all $k=1,2,\dots,n$ are so-called *weights*. NURBS curves can then be defined by the NURBS basis functions and control points $\mathbf{P}_i$ as $$\begin{aligned}
 \mathbf{C}(\xi) = \sum_{i=1}^{n}\mathbf{P}_i N_{i,p}(\xi).\end{aligned}$$ A NURBS curve can be locally modified by moving the control points. This is visualized in Fig. \[fig:controlpoints\].

[0.49]{}

[0.49]{}

Isogeometric Analysis {#sec:iga}
---------------------

Using the same functions for the representation of the geometry as in CAD software, i.e., NURBS, has the advantage that there is no need for the construction of a finite element geometry, i.e., the mesh. In isogeometric analysis the geometry is represented by a smooth mapping $$\begin{aligned}
 \mathbf{F}: \hat{{D}}\rightarrow {D},\end{aligned}$$ using NURBS as basis functions, where $\hat{{D}}$ is the reference domain, i.e., in the 2D case the unit square, and ${D}$ is the physical domain. The mapping $\mathbf{F}$ can be directly obtained from CAD software. However, not all geometries can be represented by a regular transformation of the reference domain $\hat{{D}}=[0,1]^n$, e.g., domains with a hole. In this case, a multipatch approach can be used to represent the physical domain. The physical domain ${D}$ is subdivided into $k$ patches ${D}_k$ which can each be represented by a regular transformation of the reference domain ${D}_k=\mathbf{F}_k(\hat{{D}})$. A visualization is given in Fig. \[fig:multipatch\]. The patches have a consistent discretization with a one-to-one matching of the degrees of freedom at the interfaces as in classical FEM, leading to a $C^0$ smoothness of the solution at the patch interfaces.

(0,0) rectangle (1,1) node (ref) \[pos = 0.5\] [$\hat{{D}}$]{}; (0,0)–(0,1.5); (0,0)–(1.5,0); (0) at (0,0) \[anchor=north east\] [$0$]{}; (1,2pt)–(1,-2pt) node \[anchor=north\] [$1$]{}; (2pt,1)–(-2pt,1) node \[anchor=east\] [$1$]{}; () at (0,0) \[anchor=north east\] [$0$]{};

(0,0) node (middle) circle (1.5);

(-1.5,0)–(1.5,0); (0,-1.5)–(0,1.5);

(p1) at (1,1) [${D}_1$]{}; (p2) at (-1,1) [${D}_2$]{}; (p3) at (-1,-1) [${D}_3$]{}; (p4) at (1,-1) [${D}_4$]{};

(0,0) circle (0.5);

(ref) to\[in=110, out=30\] node \[pos=0.5, anchor=south east\] [$\mathbf{F}_1(\hat{{D}})$]{} (p1); (ref) to\[in=170, out=20\] node \[pos=0.5, anchor=north\] [$\mathbf{F}_2(\hat{{D}})$]{} (p2); (ref) to\[in=-170, out=-20\] node \[pos=0.5, anchor=south\] [$\mathbf{F}_3(\hat{{D}})$]{} (p3); (ref) to\[in=-110, out=-30\] node \[pos=0.5, anchor=north east\] [$\mathbf{F}_4(\hat{{D}})$]{} (p4);

Thus, the geometry of the CAD models can be exactly represented by a multipatch model [@Buffa_2015aa] and there is no need for the generation of a computational mesh that approximates the geometry, e.g., a triangulation. For the analysis, IGA uses the same setting as the classical finite element analysis with the exception of using B-splines or NURBS as basis and test functions $w_j$ for $j=1,\dots , N_q$.

Numerical Shape Optimization {#sec:num_shape_opti}
============================

We solve the optimization problem – by means of a gradient-based shape optimization algorithm based on the shape derivative . By solving an auxiliary boundary value problem of the form $$\begin{aligned}
 \label{eq_aux_prob}
    b(\mathbf{W}, \mathbf{Z}) = d \mathcal{J}(\Omega; \mathbf{Z}) \quad \forall \mathbf{Z}\end{aligned}$$ with some positive definite bilinear form $b(\cdot, \cdot)$ satisfying $b(\mathbf{Z},\mathbf{Z})>0$ for all vector fields $\mathbf{Z}$, we can extract a shape gradient $\mathbf{W}$, which satisfies $d \mathcal{J}(\Omega; \mathbf{W}) = b(\mathbf{W},\mathbf{W}) >0$. Thus, moving the control points of the motor geometry a small distance $\delta$ into the direction of $\mathbf{W}$ will yield an increase of the objective function $\mathcal{J}$. Likewise, since $d\mathcal{J}(\Omega;\mathbf{W})$ is linear in $\mathbf{W}$, a decrease can be achieved by moving the control points into the direction of the negative shape gradient $-\mathbf{W}$. The auxiliary boundary value problem can be interpreted as finding a Riesz representative $\mathbf{W}$ of the functional $d \mathcal{J}(\Omega, \cdot)$ with respect to the metric given by $b(\cdot, \cdot)$. Of course, here different bilinear forms $b(\cdot, \cdot)$ can be chosen which amount to shape gradients in different metrics. In our algorithm, we choose $$\begin{aligned}
 \label{eq_b_W_Z}
    b(\mathbf{W}, \mathbf{Z} ) = \int_D \mathbf{DW} : \mathbf{DZ} + \mathbf{W} \cdot \mathbf{Z} \; {\mathrm{d}} x\end{aligned}$$ where we used $\mathbf{A}:\mathbf{B} = \sum_{i,j=1}^n A_{ij} B_{ij}$ denotes the Frobenius inner product for two matrices $\mathbf{A}, \mathbf{B} \in \mathbb R^n$ and $\mathbf{DW}$, $\mathbf{DZ}$ denote the Jacobi matrices of $\mathbf{W}$, $\mathbf{Z}$, respectively. Our algorithm consists in iteratively determining a descent vector field $\mathbf{W}$ by solving the auxiliary boundary value problem with the bilinear form $b(\cdot, \cdot)$ given by and moving the domain a distance $\delta$ into this direction. The step size $\delta$ is chosen as the maximum of the set $\{1, 1/2, 1/4 , \dots \}$ such that, no intersections of the patches occur and the objective value is decreased. When no further improvement can be achieved, the algorithm terminates. An overview of the shape optimization algorithm is given in Fig. \[fig:shape\_opti\_overview\]. = \[diamond, aspect=4, draw, fill=green!7, text badly centered, inner sep=0pt\] = \[rectangle, draw, fill=green!3, text width=0.5, text centered, rounded corners, minimum height=2em\] = \[draw, -latex’\]

Results of the Optimization {#sec:results}
===========================

The method described in Section \[sec:num\_shape\_opti\] is used to minimize the total harmonic distortion of the electromotive force by shape optimizing the rotor of a 6-pole permanent magnet synchronous machine (for the description of geometry and material coefficients see [@Bontinck_2018af Chapter V.A] and to download the geometry of the machine as an IGES file see [@Merkel_2020git]) in generator mode under no load condition. The shape and position of the permanent magnets in the rotor are fixed. The machine is discretized by $N_{\mathrm{DoF}}=4354$ degrees of freedom using B-splines of degree $p=2$ as basis functions and linearized material laws. The implementation is based on GeoPDEs [@Falco_2011aa]. The level of refinement is chosen such that the relative error in the $L^2$-norm of the magnetic vector potential $u$ in the domain is smaller than $\si{10^{-3}}$. The coupling of stator and rotor domain is realized using ${N_\Gamma} = 36$ harmonic basis functions. The simulation of a rotation of $\SI{120}{\degree}$ with $N_{\alpha}=120$ takes about $\SI{10}{\second}$. The optimization process takes about $\SI{2}{h}$ and is carried out in Matlab^^ R2019a on a 6-core machine (Intel^^ Core i7-5820K CPU) with RAM. The original and optimized design of the rotor can be seen in Fig. \[fig:rotor\_design\_opt\]. The total harmonic distortion is reduced by more than $75\,\%$ in 59 iterations. The total harmonic distortion in the iterations of the shape optimization algorithm and the Fourier coefficients of the original design compared to the Fourier coefficients of the optimized design can be seen in \[fig:THD\_opt,fig:fourier\_coef\_opt\].

[0.49]{}

at (0,0) [![Results of shape optimization in generator mode under no load condition. The shape of the rotor was optimized, minimizing the total harmonic distortion. Shape and position of the permanent magnet was fixed. Figure (a) shows the original design, figure (b) shows the optimized design. The black line shows the changes in the rotor amplified by a factor of 5. \[fig:rotor\_design\_opt\]](images/original_rotor_iter0.png "fig:"){width="\textwidth"}]{};

[0.49]{}

at (0,0) [![Results of shape optimization in generator mode under no load condition. The shape of the rotor was optimized, minimizing the total harmonic distortion. Shape and position of the permanent magnet was fixed. Figure (a) shows the original design, figure (b) shows the optimized design. The black line shows the changes in the rotor amplified by a factor of 5. \[fig:rotor\_design\_opt\]](images/optimized_rotor_iter59.png "fig:"){width="\textwidth"}]{}; [ at (0,0) [![Results of shape optimization in generator mode under no load condition. The shape of the rotor was optimized, minimizing the total harmonic distortion. Shape and position of the permanent magnet was fixed. Figure (a) shows the original design, figure (b) shows the optimized design. The black line shows the changes in the rotor amplified by a factor of 5. \[fig:rotor\_design\_opt\]](images/exaggerated_rotor_iter59x5outline_thick_airgap.png "fig:"){width="\textwidth"}]{}; ]{}

table \[x index=0, y index=1, col sep=comma\] [images/thd\_iter1-60.csv]{}; table \[x index=0, y index=1, col sep=comma\] [images/thd\_jmag.csv]{}; table \[x index=0, y index=1, col sep=comma\] [images/thd\_jmag\_nonlinear.csv]{};

table \[restrict y to domain\*=0:1.16e-2,x index=0, y index=1, col sep=comma\] [images/harmonics\_fourierCoefficients\_iter1-60.csv]{}; table \[restrict y to domain\*=0:1.19e-2,x index=0, y index=1, col sep=comma\] [images/harmonics\_fouriercoeffs\_jmag\_iter0.csv]{}; table \[restrict y to domain\*=0:1.18e-2,x index=0, y index=60, col sep=comma\] [images/harmonics\_fourierCoefficients\_iter1-60.csv]{}; table \[restrict y to domain\*=0:1.19e-2,x index=0, y index=1, col sep=comma\] [images/harmonics\_fouriercoeffs\_jmag\_iter59.csv]{}; (0) at (axis cs:0.4,0.0126) \[rotate=90\] [[0.132]{}]{}; (0) at (axis cs:0.8,0.0129) \[rotate=90\] [[0.135]{}]{}; (0) at (axis cs:1.2,0.0128) \[rotate=90\] [[0.134]{}]{}; (0) at (axis cs:1.6,0.0129) \[rotate=90\] [[0.135]{}]{};

  ---------------- ------------------------------------ ------------------------------------
                    $\operatorname{THD}({\mathcal E})$   $\operatorname{THD}({\mathcal E})$
                             original design                      optimized design
     IGA linear                  0.099268                             0.024793
    JMAG linear                  0.099754                             0.02234
   JMAG nonlinear                0.10639                              0.034106
  ---------------- ------------------------------------ ------------------------------------

  : Total harmonic distortion of the electromotive force ${\mathcal E}$ for original and the optimized rotor design computed with IGA and with JMAG^^ for linear and nonlinear material laws.[]{data-label="table:thd"}

Validation
==========

Exploiting the geometry representation which is given directly in NURBS, i.e., the representation commonly used in computer-aided design (CAD), the model can be easily imported into computer-aided-engineering (CAE) software, e.g., using the IGES (Initial Graphics Exchange Specification) format. The results of the optimization in \[sec:results\] is validated by using JMAG^^ as reference solution. The JMAG^^ model has 281198 elements of order $1$ and the meshing takes about $\SI{15}{\second}$. The simulation of a single rotation of $\SI{120}{\degree}$ with $N_{\alpha}=120$ takes about $\SI{30}{\minute}$. The results of the JMAG^^ simulation using linearized material laws is shown in \[fig:THD\_opt,fig:fourier\_coef\_opt\]. They are in very good agreement with the results from GeoPDEs for both geometries. To investigate the impact of saturation on the results, the material of the rotor is modeled in JMAG^^ by the nonlinear material curve M530-50A for both, the original and the optimized geometry. The results for the THD of the EMF can be seen in Fig. \[fig:THD\_opt\] and \[table:thd\] where still a reduction of $68\,\%$ is achieved. The comparison of the Fourier coefficients of the electromotive force for the original and the optimized machine considering nonlinear material laws is depicted in Fig. \[fig:fourier\_coef\_nonlinear\].

table \[restrict y to domain\*=0:1.05e-2,x index=0, y index=1, col sep=comma\] [images/harmonics\_fouriercoeffs\_jmag\_iter0\_nonlinear.csv]{}; table \[restrict y to domain\*=0:1.03e-2,x index=0, y index=1, col sep=comma\] [images/harmonics\_fouriercoeffs\_jmag\_iter59\_nonlinear.csv]{}; (0) at (axis cs:0.6,0.0117) \[rotate=90\] [[0.086]{}]{}; (0) at (axis cs:1.4,0.0115) \[rotate=90\] [[0.083]{}]{};

In conclusion, the consideration of fully nonlinear material laws increases the THD but the optimization using the linear model is still very effective.

Conclusion
==========

This paper has proposed a freeform shape optimization methodology based on isogeometric analysis. The optimized geometry exhibits a significantly reduced total harmonic distortion, i.e., by a factor of four. It has been demonstrated that the spline-based geometry representation allows easy integration into existing workflows using proprietary software.

Acknowledgements {#acknowledgements .unnumbered}
================

This work is supported by the German BMBF by the PASIROM project (grant nr. 05M2018RDA), the ‘Excellence Initiative’ of the German Federal and State Governments and by the Graduate School of Computational Engineering at Technische Universität Darmstadt. The authors thank C. Mellak, H. De Gersem, Z. Bontinck and J. Corno for their help and the fruitful discussions.

[1]{}

T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs, “Isogeometric analysis: [CAD]{}, finite elements, [NURBS]{}, exact geometry and mesh refinement,” *Comput. Meth. Appl. Mech. Eng.*, vol. 194, pp. 4135–4195, 2005. DOI: 10.1016/j.cma.2004.10.008.

A. Buffa, G. Sangalli, and R. Vázquez, “Isogeometric analysis in electromagnetics: B-splines approximation,” *Comput. Meth. Appl. Mech. Eng.*, vol. 199, pp. 1143–1152, 2010. DOI: http://dx.doi.org/10.1016/j.cma.2009.12.002.

Z. Bontinck, J. Corno, H. De Gersem, S. Kurz, A. Pels, S. Schöps, F. Wolf, C. de Falco, J. Dölz, R. Vázquez, and U. Römer, “Recent advances of isogeometric analysis in computational electromagnetics,” *ICS Newsletter (International Compumag Society)*, vol. 3, 2017.

J. Dölz, H. Harbrecht, S. Kurz, S. Schöps, and F. Wolf, “A fast isogeometric [BEM]{} for the three dimensional [Laplace]{}- and [Helmholtz]{} problems,” *Comput. Meth. Appl. Mech. Eng.*, vol. 330, pp. 83–101, 2018. DOI: 10.1016/j.cma.2017.10.020.

P. Gangl, U. Langer, A. Mantzaflaris, and R. Schneckenleitner, “Isogeometric simulation and shape optimization with applications to electrical machines,” 2018, arXiv e-print.

P. Di Barba, *Multiobjective Shape Design in Electricity and Magnetism*, ser. Lecture Notes in Electrical Engineering.1em plus 0.5em minus 0.4emSpringer, 2010.

Y. Duan and D. M. Ionel, “A review of recent developments in electrical machine design optimization methods with a permanent-magnet synchronous motor benchmark study,” *IEEE Trans. Ind. Appl.*, vol. 49, no. 3, pp. 1268–1275, 2013. DOI: 10.1109/TIA.2013.2252597.

S. Russenschuck, “Mathematical optimization techniques for the design of permanent magnet synchronous machines based on numerical field calculation,” *IEEE Trans. Magn.*, vol. 26, no. 2, pp. 638–641, 1990.

K. Weeber and S. R. H. Hoole, “Geometric parametrization and constrained optimization techniques in the design of salient pole synchronous machines,” *IEEE Trans. Magn.*, vol. 28, no. 4, pp. 1948–1960, 1992.

N. Takorabet, B. Laporte, and G. Vinsard, “On the optimization of linear induction devices,” *Electr. Eng.*, vol. 80, pp. 221–226, 1997.

K. Hameyer and M. Kasper, “Shape optimization of a fractional horse-power [DC]{}-motor by stochastic methods,” in *Computer Aided Optimimum Design of Structures [III]{}: Optimization of Structural Systems and Applications*, S. Hernandez and C. Brebbia, Eds., 1993, pp. 15–30.

C. L. Lok, B. Vengadaesvaran, and S. Ramesh, “Implementation of hybrid pattern search–genetic algorithm into optimizing axial-flux permanent magnet coreless generator [(AFPMG)]{},” *Electr. Eng.*, vol. 99, pp. 751–761, 2017. DOI: 10.1007/s00202-016-0443-9.

C. Ma and L. Qu, “Multiobjective optimization of switched reluctance motors based on design of experiments and particle swarm optimization,” *IEEE Trans. Energ. Convers.*, vol. 30, no. 3, pp. 1144–1153, 2015.

U. Baumgartner, C. Magele, and W. Renhart, “Pareto optimality and particle swarm optimization,” *IEEE Trans. Magn.*, vol. 40, no. 2, p. 1172, 2004.

S. L. Ho, S. Yang, L. E. W. C. Ni, Guangzheng, and H. C. Wong, “A particle swarm optimization-based method for multiobjective design optimizations,” *IEEE Trans. Magn.*, vol. 41, no. 5, pp. 1756–1759, 2005. DOI: 10.1109/TMAG.2005.846033.

B. N. Cassimere and S. D. Sudhoff, “Population-based design of surface-mounted permanent-magnet synchronous machines,” *IEEE Trans. Energ. Convers.*, vol. 24, no. 2, pp. 338–346, 2009.

M. L. Bash and S. D. Pekarek, “Modeling of salient-pole wound-rotor synchronous machines for population-based design,” *IEEE Trans. Energ. Convers.*, vol. 26, no. 2, pp. 381–392, 2011.

G. Y. Sizov, P. Zhang, D. M. Ionel, N. A. O. Demerdash, and M. Rosu, “Automated multi-objective design optimization of [PM]{} [AC]{} machines using computationally efficient [FEA]{} and differential evolution,” *IEEE Trans. Ind. Appl.*, vol. 49, no. 5, pp. 2086–2096, 2013.

Z. Bontinck, J. Corno, S. Schöps, and H. De Gersem, “Isogeometric analysis and harmonic stator-rotor coupling for simulating electric machines,” *Comput. Meth. Appl. Mech. Eng.*, vol. 334, pp. 40–55, 2018. DOI: 10.1016/j.cma.2018.01.047.

S. Schöps, H. De Gersem, and T. Weiland, “Winding functions in transient magnetoquasistatic field-circuit coupled simulations,” *COMPEL*, vol. 32, no. 6, pp. 2063–2083, 2013. DOI: 10.1108/COMPEL-01-2013-0004.

P. Gangl, U. Langer, A. Laurain, H. Meftahi, and K. Sturm, “Shape optimization of an electric motor subject to nonlinear magnetostatics,” *SIAM J. Sci. Comput.*, vol. 37, no. 6, pp. B1002–B1025, 2015. DOI: 10.1137/15100477X.

P. Hansbo, C. Lovadina, I. Perugia, and G. Sangalli, “A lagrange multiplier method for the finite element solution of elliptic interface problems using non-matching meshes,” *Numer. Math.*, vol. 100, no. 1, pp. 91–115, 2005. DOI: 10.1007/s00211-005-0587-4.

H. De Gersem and T. Weiland, “Harmonic weighting functions at the sliding interface of a finite-element machine model incorporating angular displacement,” *IEEE Trans. Magn.*, vol. 40, no. 2, pp. 545–548, 2004. DOI: 10.1109/TMAG.2004.824616.

H. De Gersem, J. Gyselinck, P. Dular, K. Hameyer, and T. Weiland, “Comparison of sliding-surface and moving-band techniques in frequency-domain finite-element models of rotating machines,” *COMPEL*, vol. 23, no. 4, pp. 1006–1014, 2004.

L. Piegl and W. Tiller, *The [NURBS]{} Book*, 2nd ed.1em plus 0.5em minus 0.4emSpringer, 1997.

A. Buffa, R. H. Vázquez, G. Sangalli, and L. B. da Veiga, “Approximation estimates for isogeometric spaces in multipatch geometries,” *Numer. Meth. Part. Differ. Equat.*, vol. 31, no. 2, pp. 422–438, 2015. DOI: 10.1002/num.21943.

Z. Bontinck, “Numerical methods for the estimation of the impact of geometric uncertainties on the performance of electromagnetic devices,” 2018.

M. Merkel, P. Gangl, and S. Schöps, “Results of the shape optimization of rotating electric machines using isogeometric analysis,” 2020. DOI: 10.5281/zenodo.3688479.

C. de Falco, A. Reali, and R. Vázquez, “[GeoPDEs]{}: A research tool for isogeometric analysis of [PDE]{}s,” *Advances in Engineering Software*, vol. 42, pp. 1020–1034, 2011. DOI: 10.1016/j.advengsoft.2011.06.010.
