1 Introduction↩︎

Electric machines are widely used in modern applications and their efficient numerical simulation remains an important topic in computational engineering. In this work, we study permanent magnet synchronous machines (PMSM) and exploit their periodic structure by restricting the analysis to a single pole (cf. [1]). Due to their layered, rotationally symmetric design and the clear separation into rotor, stator and air gaps, electric motors are naturally well-suited to domain decomposition approaches such as hybrid mixed domain decomposition (HMDD), as is evident from Fig. 1.

The main contribution of this work is the application of HMDD to a magnetostatic rotor-stator coupling in a geometry-respecting high-order finite element setting. In addition, we transfer the well-posedness results and error estimates from [3] and present an academic proof-of-concept example.

2 Model Problem↩︎

We consider a two dimensional mixed Poisson equation derived from Maxwell’s equations for magnetostatics \[\begin{align} \mathop{\mathrm{\boldsymbol{curl}}}\vec{H} &=& \vec{J},\\ \mathop{\mathrm{div}}\vec{B} &=& 0,\label{Seibel:eq:MXWdivB0} \end{align}\tag{1}\] where \(\vec{H}\) is the magnetic field strength, \(\vec{B}\) the magnetic flux density and \(\vec{J}\) the electric current density. Further, we assume the affine material law \[\vec{B} = \mu\vec{H} + \mu_0\vec{M},\] where \(\vec{M}\) is the remanent magnetization of the permanent magnet [4]. The magnetic vector potential \(\vec{A}\) implicitly fulfilling 1 is introduced as \[\vec{B} = \mathop{\mathrm{\boldsymbol{curl}}}\vec{A}.\] Assuming that electric current density and vector potential point in \(z\)-direction only, we can reduce the problem to two dimensions [1]. The resulting equations read

\[\begin{align} \mu\vec{h} - \mathop{\mathrm{\boldsymbol{grad}}}a_z &= -\mu_0\vec{m},\label{Seibel:eq:MPP}\\ -\mathop{\mathrm{div}}\vec{h} &= j_z,\nonumber \end{align}\tag{2}\] with the \(z\)-component of the magnetic vector potential \(a_z\), the rotated magnetic field strength \(\vec{h}=(-H_y,H_x)^\top\), the permeability \(\mu\), the vacuum permeability \(\mu_0\), the rotated magnetization \(\vec{m}=(-M_y,M_x)^\top\) and the \(z\)-component of the electric current density \(j_z\). Note, first, that we used the identity \[\mathop{\mathrm{\text{curl}_\text{2D}}}\vec{v} = \partial_y v_x - \partial_x v_y = -\mathop{\mathrm{div}}\vec{v}^\perp,\] where \(\vec{v}^\perp = (v_y,-v_x)^\top\). Second, in contrast to the formulation as a single second order PDE (e.g. (3.8) [1]), we do not apply any differential operator to the magnetization \(\vec{m}\). This is beneficial as \(\vec{m}\) is not necessarily continuous and, hence, we circumvent the appearance of distributional derivatives.

2.1 Domain Decomposition and Boundary Conditions↩︎

The computational domain is split by the interface \(\Gamma\), that is located in the middle of the air gap, into two disjoint subdomains \(\Omega_\text{R}\) and \(\Omega_\text{S}\) for the rotor and stator, respectively. The outer boundary of the electric machine and, hence, the stator is denoted by \(\Sigma_\text{S}\) while the boundary between rotor and shaft is denoted by \(\Sigma_\text{R}\). The boundaries towards neighboring poles are represented by \(\Sigma^-\) and \(\Sigma^+\). See Fig. 2 for a sketch of the geometry.

On the interface \(\Gamma\) we choose boundary conditions respecting the tangential continuity of the magnetic field strength. Note, that \(\vec{h}\) is rotated by 90 ° and, thus, being continuous in the normal but not tangential direction. On the outer boundaries \(\Sigma_\text{S}\) and \(\Sigma_\text{R}\) the magnetic vector potential \(a_z\) fulfills homogeneous Dirichlet boundary conditions, while \(a_z\) and \(\vec{h}\cdot\vec{n}\) are anti-periodic on \(\Sigma^-\) and \(\Sigma^+\):

\[\begin{alignat*}{2} a_z^+ &= a_z^- &&\quad\text{on } \Gamma,\\ \vec{h}^+\cdot\boldsymbol{n}_\Gamma &= \vec{h}^-\cdot\boldsymbol{n}_\Gamma &&\quad\text{on } \Gamma,\\ a_z &= 0 &&\quad\text{on }\Sigma_\text{S},\\ a_z &= 0 &&\quad\text{on }\Sigma_\text{R},\\ a_z\left(\tfrac\pi3,r\right) &= -a_z\left(0,r\right) &&\quad\forall (0,r)\in\Sigma^-\\ \vec{h}\cdot\vec{n}_{\Sigma}^+\left(\tfrac\pi3,r\right) &= \vec{h}\cdot\vec{n}_\Sigma^-\left(0,r\right) &&\quad\forall (0,r)\in\Sigma^-. \end{alignat*}\] By \(a_z^\pm\) we denote the left- and right-handed traces of \(a_z\) on \(\Gamma\), \(\Sigma^+\) and \(\Sigma^-\), by \(\vec{n}_\Sigma^-\) and \(\vec{n}_{\Sigma}^+\) the outward normal vectors on \(\Sigma^-\) and \(\Sigma^+\), respectively.

3 Hybrid Mixed Domain Decomposition↩︎

The most relevant part in the field of domain decomposition is the (re)coupling of local solutions to obtain a global solution. In this work, we use the hybrid mixed domain decomposition (HMDD) method introduced in [3]. The formulation can be transferred to the present magnetostatic rotor-stator problem without modification. The advantage of this approach is the reduced number of globally coupled degrees of freedom and the absence of problems at cross sections of subdomains as they arise in other domain decomposition methods such as FETI.

3.1 HMDD Method↩︎

Inspired by hybridized discontinuous Galerkin (HDG) methods (see e.g. [5]) we introduce a hybrid variable \(\lambda\), approximating the trace of \(a_z\) on the interface \(\Gamma\), to couple the subdomain problems and a penalization parameter \(\tau\) is introduced to control the continuity of the finite element solution. Using the mixed problem 2 in combination with the boundary conditions discussed in Sec. 2.1, we obtain the following formulation:
Seek \((\vec{h}_h,a_{z,h},\lambda_h)\in W_h\times Q_h\times M_h\) such that \[\begin{align} \int_{{\Omega\setminus\Gamma}}\mu \vec{h}_h\cdot\vec{h}_h^\prime + a_{z,h} \mathop{\mathrm{div}}\vec{h}_h^\prime \,\mathrm{d}\vec{x} + \int_{\Gamma_\Sigma} \lambda_h \big\llbracket \vec{h}_h^\prime\cdot\vec{n} \big\rrbracket\,\mathrm{d}\sigma\nonumber &= -\mu_0\int_{{\Omega\setminus\Gamma}}\vec{m}\cdot\vec{h}_h^\prime\,\mathrm{d}\vec{x},\\ -\int_{{\Omega\setminus\Gamma}} \mathop{\mathrm{div}}\vec{h}_h a_{z,h}^\prime \,\mathrm{d}\vec{x} + \sum_\pm \int_{\Gamma_\Sigma}\tau(a_{z,h}^\pm-\lambda_h) a_{z,h}^{\prime\pm} \,\mathrm{d}\sigma &= \int_{{\Omega\setminus\Gamma}} j_z a_{z,h}^\prime \,\mathrm{d}\vec{x} ,\label{Seibel:eq:FEMFormulation}\\ -\int_{\Gamma_\Sigma}\big\llbracket \vec{h}_h\cdot\vec{n} \big\rrbracket\lambda_h^\prime\,\mathrm{d}\sigma - \sum_\pm \int_{\Gamma_\Sigma}\tau(a_{z,h}^\pm-\lambda_h)\lambda_h^\prime\,\mathrm{d}\sigma &= 0,\nonumber \end{align}\tag{3}\] for all \((\vec{h}_h^\prime,a_{z,h}^\prime,\lambda_h^\prime)\in W_h\times Q_h\times M_h\).

Here, \({\Omega\setminus\Gamma}= \Omega_\text{S}\cup\Omega_\text{R}\), \(\Gamma_\Sigma = \Gamma\cup\Sigma^+\cup\Sigma^-\), and \(\big\llbracket \vec{h}_h\cdot\vec{n} \big\rrbracket\) and \(\big\llbracket \vec{h}_h^\prime\cdot\vec{n} \big\rrbracket\) denote the jump of \(\vec{h}_h\) and \(\vec{h}_h^\prime\) across \(\Gamma\) or \(\Sigma^\pm\), respectively. Formulation 3 fits into the HMDD framework [3], hence, we can transfer the properties of the HMDD method to 3 in the following.

3.2 Finite Element Spaces↩︎

In the finite element formulation 3 , the space \(W_h\) for the flux density is chosen to be the Raviart-Thomas space of order \(q\) on \({\Omega\setminus\Gamma}\), where discontinuities across \(\Gamma\) are permitted. The spaces \(Q_h\) and \(M_h\) are discontinuous, piecewise polynomial spaces of order \(q\) on \({\Omega\setminus\Gamma}\) and \(\Gamma\), respectively.

In fact, we choose the same spaces as for the finite element part in [3], where we also discussed spaces for a more general Galerkin or even continuous setting.

3.3 Well-Posedness↩︎

In [3] we were able to show that the HMDD method is well-posed, both in the continuous as well as in the discrete setting including finite elements:

For any \(f \in L^2({\Omega\setminus\Gamma})\) there exist unique solutions \((\vec{h}_h, a_{z,h}, \lambda_h) \in W_h \times Q_h \times M_h\) of 3 where it holds with a constant \(C\) independent of \(\tau\) that \[\tag{4} \begin{align} \big\|\vec{h}_h\big\|_{H(\mathop{\mathrm{div}},\Omega)} + \big\|a_{z,h}\big\|_{L^2({\Omega\setminus\Gamma})} + \big\|\lambda_h\big \|_{L^2(\Gamma)} &\leq C \big\|f \big\|_{L^2({\Omega\setminus\Gamma})}, \\ \big\|\big\llbracket \vec{h}_h\cdot\vec{n} \big\rrbracket\big\|_{L^2(\Gamma)} &\leq C \sqrt{\tau} \|f\|_{L^2({\Omega\setminus\Gamma})},\tag{5}\\ \sqrt{\tau} \left( \big\|\big\llbracket a_{z,h} \big\rrbracket\big\|_{L^2(\Gamma)} + \big\| \left\{a_{z,h}\right\} - \mu_h\big\|_{L^2(\Gamma)} \right) &\leq C \|f\|_{L^2({\Omega\setminus\Gamma})}\tag{6},\\ \sum_\pm \tfrac{1}{\sqrt{1+\tau}} \big\| \sqrt{\tau} a_{z,h}^\pm \big\|_{L^2(\Gamma)} &\leq C \big\|f \big\|_{L^2({\Omega\setminus\Gamma})}. \end{align}\] Note, that the estimates partially depend on the stabilization parameter \(\tau\). From 5 we deduce that the continuity of \(\vec{h}_h\cdot\vec{n}\) is increasingly enforced for decreasing \(\tau\), while 6 shows the continuity of \(a_{z,h}\) being increasingly enforced for increasing \(\tau\).

3.4 Error Estimates↩︎

We, further, conjectured error estimates in the case of sufficiently smooth solutions [3], which we can assume to hold for 3 as well: \[\begin{align} \| a_{z} - a_{z,h} \|_{L^2({\Omega\setminus\Gamma})} + \| \mu - \mu_h \|_{L^2(\Gamma)} &\leqC_1 h^{q+1}, \tag{7}\\ \| \vec{h} - \vec{h}_h \|_{L^2({\Omega\setminus\Gamma})} &\leqC_1 h^{q+1} + \tfrac{C_2 h\tau}{C_3 + h\tau} h^{q+\frac{1}{2}}, \tag{8} \\ \| \mathop{\mathrm{div}}(\vec{h} - \vec{h}_h) \|_{L^2({\Omega\setminus\Gamma})} &\leq C_1 h^{q+1} + \tfrac{C_2 h\tau}{C_3 + h\tau} h^{q-\frac{1}{2}}, \tag{9} \\ \| \big\llbracket \vec{h}_h\cdot \vec{n} \big\rrbracket \|_{L^2(\Gamma)} &\leq \tfrac{C_2 h \tau}{C_3 + h\tau} h^{q}, \tag{10} \\ \| \big\llbracket a_{z,h}^\pm \big\rrbracket \|_{L^2(\Gamma)} &\leq \tfrac{1}{(C_1)^{-1} + h\tau} h^{q+1}. \tag{11} \end{align} \tag{12}\] We find that the convergence rates depend on the product of the mesh width and the stabilization parameter \(\tau\). For a more extensive discussion of this behavior we refer the reader to [3].

4 Numerical Implementation in Concepts↩︎

The experiments are performed using the numerical C++-library Concepts 2 [6], which provides interfaces to various direct solvers solvers – for this work we used SuperLU. It supports being run in parallel, and is suitable for high-order and \(hp\)-adaptive FEM [7]. In particular, Concepts supports quadrilateral mesh elements with circularly curved boundaries, which is beneficial for exact geometry representation. It can read gmsh mesh files and we used the interface to MATLAB for graphical representation. The initial mesh as depicted in Fig. 3.

We perform an experiment that serves as a proof of concept reinforcing the applicability of HMDD for the simulation of a PMSM. Therefore, we assume that adjacent slots pairwise share a phase in the current density. The phase between two neighboring pairs is shifted by 120 °. We, further, assume the current density to be \(\left|j_z\right|=\SI{5}{\ampere\per\milli\meter\squared}\) and the remanence being \(\left|\mu_0\vec{m}\right|=\SI{1}{\tesla}\). The results are compared with those obtained by an in-house iso-geometric analysis (IGA) code [8] to assess the approach’s accuracy.

First, we consider the magnitude of the (rotated) magnetic field density \(|\vec{b}|\). On the left-hand side of Fig.4 we see the results obtained by the in-house IGA code while the right-hand subfigure shows the result of the HMDD method implemented in Concepts. Second, we visualize the magnetic potential lines given by the iso-potential lines of \(a_{z,h}\) in Fig. 5. Matching the left- with the right-hand side of Fig. 4 and Fig. 5, respectively, we find the results to be overall in good agreement besides a visualization issue in the center of the magnet.

5 Conclusion & Outlook↩︎

We have derived a finite element method from a mixed magnetostatics problem in two dimensions that fits into the HMDD framework. The properties of the HMDD method are naturally transported to the considered example which then was implemented in Concepts. In an academic example, serving as proof of concept, we found good agreement with an in-house IGA code.
In a next step, we aim to advance to magnetoquasistatics with rotational motion by using a time-stepping procedure.
Regarding motion, we have the reasoned hope that the presented formulation is well-suited for mortaring as this is inherently supported by HDG methods for \(0~<~\tau~<~\infty\) [5].

Acknowledgements↩︎

This work is supported by the Graduate School CE within the Profile Topic Computational Engineering at the Technical University of Darmstadt. The authors thank Michael Wiesheu, funded by the Collaborative Research Centre – TRR361/F90: CREATOR, for the provision of Fig. 1 and the comparison results in Fig. 4 and Fig. 5. The first author thanks Christian Bergfried, also funded by CREATOR, for fruitful discussions and helpful insights on this topic.

References↩︎