Shape Optimization of Rotating Electric Machines using Isogeometric Analysis and Harmonic Stator-Rotor Coupling

TL;DR

This paper presents a shape optimization method for rotating electric machines utilizing isogeometric analysis and harmonic coupling, enabling efficient geometry modifications and improved electromagnetic performance.

Contribution

It introduces a novel multipatch isogeometric modeling approach with harmonic stator-rotor coupling for optimized electric machine design.

Findings

Reduced harmonic distortion in electromotive force

Efficient geometry modification without remeshing

Validated approach on a 6-pole permanent magnet machine

Abstract

This work deals with shape optimization of electric machines using isogeometric analysis. Isogeometric analysis is particularly well suited for shape optimization as it allows to easily modify the geometry without remeshing the domain. A 6-pole permanent magnet synchronous machine is modeled using a multipatch isogeometric approach and rotation of the machine is realized by modeling the stator and rotor domain separately and coupling them at the interface using harmonic basis functions. Shape optimization is applied to the model minimizing the total harmonic distortion of the electromotive force as a goal functional.

Tables

Table 1. TABLE I : 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.

	$\mathrm{THD}\left(ℰ\right)$	$\mathrm{THD}\left(ℰ\right)$
	original design	optimized design
IGA linear	0.099268	0.024793
JMAG linear	0.099754	0.02234
JMAG nonlinear	0.10639	0.034106

Equations

$$\nabla \times \left(\nu \nabla \times \mathbf{A}\right)={\mathbf{J}}_{\mathrm{src}}+\nabla \times \mathbf{M},$$

$${\mathrm{THD}}_I(f)=\frac{\sqrt{{\sum }_{n\in I,n\ne 1}|c_n{|}^2}}{|c_1|},$$

$$f(t)=\sum _{n=-\infty }^{\infty }{c}_n{e}^{ınt},$$

$$f(t)=\frac{A_0}{2}+\sum _{n=1}^{\infty }{A}_n\cos (nt)+B_n\sin (nt)$$

$$\mathrm{THD}(f)=\sqrt{\frac{{\sum }_{n\in I,n\ne 1}{A}_n^2+B_n^2}{A_1^2+B_1^2}}.$$

$${\mathcal{E}}_k(\mathbf{A}(t))={\partial }_t{\Psi }_k(\mathbf{A}(t)),$$

$${\Psi }_k(\mathbf{A}(t))=N_p{\int }_{D_{3D}}{\mathbf{\chi }}_k\cdot \mathbf{A}(t)\mathrm{d}{D}_{3D},$$

$$-\nabla \cdot \left(\nu \nabla u\right)=J_{\mathrm{src},z}+\nabla \times \mathbf{M}\cdot {\mathbf{e}}_z,$$

$$\underset{\Omega \in \mathcal{A}}{\min }\mathcal{J}(\Omega ):={\mathrm{THD}}_I(\mathcal{E}(u(t,\Omega ))),$$

$$\underset{=:a[\alpha (t)](\Omega ;u,v)}{\underbrace{{\int }_D{\nu }_{\Omega }(\alpha (t))\nabla u\cdot \nabla v\mathrm{d}D}}=⟨r(t,\alpha (t)),v⟩$$

$$\begin{array}{rl}⟨r(t,& \alpha (t)),v⟩={\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))}\left(\genfrac{}{}{0ex}{}{-M_2}{M_1}\right)\cdot \nabla v(x,y)\mathrm{d}{D}_{\mathrm{pm}}(\alpha (t))\end{array}$$

$$\underset{\Omega }{\text{mbox}min}\mathcal{J}(u_1,\dots ,u_{N_{\alpha }})$$

$$\displaystyle\mbox{s.t. }\left\{\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.$$

$$d\mathcal{J}(\Omega ;\mathbf{W})=\underset{\delta \searrow 0}{\text{mbox}lim}\frac{\mathcal{J}(T_{\delta }^{\mathbf{W}}(\Omega ))-\mathcal{J}(\Omega )}{\delta },$$

$$\mathcal{L}(\Omega ,$$

$$\mathcal{J}({\phi }_1,\dots ,{\phi }_{N_{\alpha }})+\sum _{k=1}^{N_{\alpha }}\left(a_k({\phi }_k,{\psi }_k)-⟨r_k,{\psi }_k⟩\right).$$

$$\frac{\partial }{\partial {\psi }_i}\mathcal{L}(\Omega ,\mathbf{u},\mathbf{\psi })(q_i)=a_i(u_i,q_i)-⟨r_i,q_i⟩=0$$

$$0=\frac{\partial }{\partial {\phi }_i}\mathcal{L}(\Omega ,\mathbf{u},\mathbf{p})(v_i)=\frac{\partial \mathcal{J}}{\partial {\phi }_i}(\mathbf{u})(v_i)+a_i(v_i,p_i)$$

$$\Psi (u_1,\dots ,u_{N_{\alpha }})=\sum _{n=-N_{\alpha }}^{N_{\alpha }}{C}_n{e}^{ınt}.$$

$$\mathcal{E}(u_1,\dots ,u_{N_{\alpha }})=\sum _{n=-N_{\alpha }}^{N_{\alpha }}\underset{c_n}{\underbrace{ın{C}_n}}{e}^{ınt},$$

$$\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).$$

$$A_k(u_1,\dots ,u_{N_{\alpha }})$$

$$B_k(u_1,\dots ,u_{N_{\alpha }})$$

$$\frac{d\mathrm{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 I,k\ne 1}{A}_k^2+B_k^2}}\left(\sum _{k\in I,k\ne 1}{A}_k{A}_k^{\prime }+B_k{B}_k^{\prime }\right)$$

$$-\frac{1}{(A_1^2+B_1^2{)}^{3/2}}\left(\sqrt{\sum _{k\in I,k\ne 1}{A}_k^2+B_k^2}\left(A_1{A}_1^{\prime }+B_1{B}_1^{\prime }\right)\right).$$

$$\frac{d{A}_k}{d{u}_i}(\mathbf{u})(v_i)$$

$$\frac{d{B}_k}{d{u}_i}(\mathbf{u})(v_i)$$

$$a_i(v_i,p_i)=-\frac{d\mathrm{THD}(\mathcal{E}(\mathbf{u}))}{d{u}_i}(\mathbf{u})(v_i)$$

$$\begin{array}{rl}d\mathcal{J}& (\Omega ;\mathbf{W})=\sum _{l=1}^{N_{\alpha }}{\int }_{{\Omega }_{\mathrm{pm}}}(\nabla \cdot (\mathbf{W})\mathbf{I}-{\mathbf{DW}}^{\top })\nabla p_l\cdot \left(\genfrac{}{}{0ex}{}{-M_2}{M_1}\right)\mathrm{d}x\\ & +\sum _{l=1}^{N_{\alpha }}{\int }_D\nu (\nabla \cdot \mathbf{W}\mathbf{I}-{\mathbf{DW}}^{\top }-\mathbf{DW})\nabla u_l\cdot \nabla p_l\mathrm{d}x.\end{array}$$

$$\displaystyle\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$$

$$\displaystyle\int_{\Gamma}\left\llbracket u\right\rrbracket\mu\operatorname{d}\Gamma$$

$$⟨r_{\mathrm{st}},v_{\mathrm{st}}⟩$$

Full text

I 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) [1], [2], [3], [4] 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 [5]. 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., [6], [7] and the references therein). Originally, mainly gradient-based optimization methods were used (see, e.g., [8, 9, 10]) but stochastic optimization algorithms became more popular during the last 20 years (see, e.g., [11], [12]). Most of the proposed algorithms use stochastic or population-based optimization, e.g., genetic algorithms and particle swarm optimization (see, e.g., [13]) which have also been extended to multi-objective optimization problems (see, e.g., [14, 15]). For permanent magnet synchronous machines, stochastic optimization methods are commonly used (see, e.g.,[16, 17, 18]).

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 [19].

The paper is structured as follows: In Section II 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 III. In Section IV 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 V.

II Model of the Electric Machine

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

[TABLE]

with the piecewise constant reluctivity $\nu$, the magnetic vector potential $\mathbf{A}$, the current density $\mathbf{J}_{\mathrm{src}}=\sum_{k}\boldsymbol{\chi}_{k}i_{k}$ given by winding functions $\boldsymbol{\chi}_{k}$ and currents $i_{k}$ [20], 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. 1.

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

[TABLE]

where $c_{n}$ are the coefficients of the Fourier series of $f(t)$, i.e.

[TABLE]

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 (3) can be written as

[TABLE]

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

[TABLE]

We are interested in the THD of the EMF. Therefore, the amplitudes $c_{n}$ are obtained by Fourier analysis of the voltages

[TABLE]

with the flux linkage

[TABLE]

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 (1) is considered on a cross section ${D}\subset\mathbb{R}^{2}$ and reduces to

[TABLE]

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_{k}i_{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.

III Optimization

III-A Problem statement

We want to optimize the $2D$ shape $\Omega$ of the $3$-phase permanent magnet synchronous machine from Fig. 2, which has $N_{p}=3$ pole pairs, in terms of the total harmonic distortion ($\operatorname{THD}$) of the electromotive force (EMF) as introduced in (5). 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 (7), i.e., $\mathcal{E}=\mathcal{E}(u(t,\Omega))$. Thus, we consider the optimization problem

[TABLE]

that is constrained by the machine model (7), which, for $t\in[0,T]$, reads in its weak form: Find $u=u(\alpha(t))$ such that

[TABLE]

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

[TABLE]

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 $\{t_{1},\dots,t_{N_{\alpha}}\}$ into ${N_{\alpha}}=120$ points in time and a corresponding discretization of the range of angular displacements into ${N_{\alpha}}$ rotor positions, $\bm{\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 (9) and (10), respectively. Approximating $u(t,\Omega)$ for $t\in[0,T]$ by $\{u_{1},\dots,u_{N_{\alpha}}\}$ 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 (8)–(9) after discretization with respect to rotor positions can be written as the optimization problem

[TABLE]

[TABLE]

which is constrained by ${N_{\alpha}}$ boundary value problems corresponding to the ${N_{\alpha}}$ rotor positions under consideration.

III-B 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

[TABLE]

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 (11)–(15), we follow the steps taken in [21]. First, we introduce the Lagrangian

[TABLE]

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 (15), it holds that

[TABLE]

for any test function $q_{i}$. Similarly, we introduce the adjoint states $p_{i}$, $i=1,\dots,{N_{\alpha}}$, as the solutions to

[TABLE]

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 (4) 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.,

[TABLE]

Using the Fourier representation (18) of the flux linkage and (5), the electromotive force can be written as

[TABLE]

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

[TABLE]

The coefficients $A_{n}$ and $B_{n}$ of the Fourier representation of $\partial_{t}{\Psi}$ (20) are obtained by time derivation and Fourier transform, which we will denote by $A_{n}=[\mathcal{F^{\prime}}_{a}({\Psi}(u_{1},\dots,u_{{N_{\alpha}}}))]_{n}$ and $B_{n}=[\mathcal{F^{\prime}}_{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^{\prime}}_{a}({\mathbf{I}})$ and $\mathbf{M}_{b}=\mathcal{F^{\prime}}_{b}({\mathbf{I}})$, where ${\mathbf{I}}$ is the identity matrix, i.e.,

[TABLE]

where $l_{z}$ is the length of the machine in $z$-direction.

In order to solve the adjoint state equation (17), 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 (4), we get

[TABLE]

Here, we used the abbreviations $A_{k}^{\prime}:=\frac{dA_{k}}{du_{i}}(\mathbf{u})(v_{i})$ and $B_{k}^{\prime}:=\frac{dB_{k}}{du_{i}}(\mathbf{u})(v_{i})$. It can be seen from (21) and (22) that

[TABLE]

Given the solutions to the forward problem (9) for all rotor positions $l\in\{1,\dots N_{\alpha}\}$, we obtain for the adjoint problem (17) defining the adjoint variable at rotor position $l\in\{1,\dots N_{\alpha}\}$: Find $p_{i}$ such that

[TABLE]

for all test functions $v_{i}$.

Finally, in a similar way as it was proposed in [21], 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})$:

[TABLE]

Here, $\mathbf{I}\in\mathbb{R}^{2\times 2}$ denotes the two-dimensional identity matrix.

IV Solution of the Forward Problem by IGA

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

IV-A Harmonic Stator-Rotor Coupling

As proposed in [19], 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 (7) can be formulated as:

Find $({{u}_{\mathrm{st}}},{{u}_{\mathrm{rt}}},\lambda)\in{V_{\mathrm{st}}}\times{V_{\mathrm{rt}}}\times{\Lambda}$ such that

[TABLE]

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 (10),

[TABLE]

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 [22].

IV-B 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 [23] as

[TABLE]

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 [24]. 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 [19],

[TABLE]

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

[TABLE]

$\mathbf{G}_{q}\in\mathbb{C}^{N_{q}\times{N_{\Gamma}}}$ are the coupling matrices with

[TABLE]

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 (27) is stable if the numbers $N_{q}$ and ${N_{\Gamma}}$ of basis functions are chosen consistently [19].

IV-C Elimination of Inner Degrees of Freedom

The system (27) can be rewritten as

[TABLE]

The internal degrees of freedom can be eliminated using the Schur-complement, giving rise to the low-dimensional interface problem

[TABLE]

with

[TABLE]

The inverses in (33) and (34) 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 (32) 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.,

[TABLE]

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 (27) has to be solved for $N_{\alpha}$ rotor positions. The computational cost of the simulation of a rotor rotation of $120\text{\,}\mathrm{\SIUnitSymbolDegree}$ with $N_{\alpha}=120$ discretized points in time using the interface problem is compared to the computational cost using the full system in Fig. 3. It can be seen that the computational time needed to solve the full system 27 depends on the number of the IGA degrees of freedom whereas the computational time needed to solve the interface problem 32 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 $120\text{\,}\mathrm{\SIUnitSymbolDegree}$ 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 33 and 34, took at most $0.678\text{\,}\mathrm{s}$ in the case of ${N_{\Gamma}}=50$ and $25746$ IGA DoFs and the postprocessing, i.e., 35 and 36, took at most $1.3012\text{\,}\mathrm{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 $16\text{\,}\mathrm{G}\mathrm{B}$ RAM.

IV-D B-splines and NURBS

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 [25]

[TABLE]

for all $p\geq 1$ and for $p=0$ via

[TABLE]

where $0/0=0$ is formally assumed. NURBS (Non-Uniform Rational B-splines) basis functions are then defined by

[TABLE]

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

[TABLE]

A NURBS curve can be locally modified by moving the control points. This is visualized in Fig. 4.

IV-E Isogeometric Analysis

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

[TABLE]

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

Thus, the geometry of the CAD models can be exactly represented by a multipatch model [26] 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}$.

V Numerical Shape Optimization

We solve the optimization problem (8)–(9) by means of a gradient-based shape optimization algorithm based on the shape derivative (23). By solving an auxiliary boundary value problem of the form

[TABLE]

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 (42) 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

[TABLE]

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 (42) with the bilinear form $b(\cdot,\cdot)$ given by (43) 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. 6.

VI Results of the Optimization

The method described in Section V 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 [27, Chapter V.A] and to download the geometry of the machine as an IGES file see [28]) 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 [29]. 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 $\mathrm{1}\mathrm{0}^{-3}$. The coupling of stator and rotor domain is realized using ${N_{\Gamma}}=36$ harmonic basis functions. The simulation of a rotation of $120\text{\,}\mathrm{\SIUnitSymbolDegree}$ with $N_{\alpha}=120$ takes about $10\text{\,}\mathrm{s}$. The optimization process takes about $2\text{\,}\mathrm{h}$ and is carried out in Matlab® R2019a on a 6-core machine (Intel® Core™ i7-5820K CPU) with $16\text{\,}\mathrm{G}\mathrm{B}$ RAM. The original and optimized design of the rotor can be seen in Fig. 7. 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 Figs. 8 and 9.

VII 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 Section VI is validated by using JMAG® as reference solution. The JMAG® model has 281198 elements of order $1$ and the meshing takes about $15\text{\,}\mathrm{s}$. The simulation of a single rotation of $120\text{\,}\mathrm{\SIUnitSymbolDegree}$ with $N_{\alpha}=120$ takes about $30\text{\,}\mathrm{min}$. The results of the JMAG® simulation using linearized material laws is shown in Figs. 8 and 9. 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. 8 and Table I 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. 10.

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

VIII 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

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.