An air-gap element for the isogeometric space-time simulation of electric machines

In the design of electric machines, simulation plays an important role in predicting machine performance and guiding design choices. Conventionally, these simulations are performed using the finite element method (FEM) with low-order elements (Monk, 2003). When considering rotating electric machines, space-time FEMs (Gangl et al., 2025) can resolve the air gap between the rotor and the stator by using simplicial space-time finite elements. However, an accurate representation of the cylindrical geometry poses challenges as such shapes are costly, e.g. in terms of number of elements, in particular in the space-time approach. Since rotor and stator are well separated, and not changing in time, it is sufficient to consider fixed spatial meshes for both, which are coupled via the air-gap domain within a time stepping approach. To avoid remeshing, techniques such as locked step, sliding surface, moving band and mortar methods have been developed, which introduce a (nonconforming) interface to couple rotating subdomains (Rodger et al., 1990; Davat et al., 1985; De Gersem et al., 2004). In particular, the mortar setting was recently extended to be applied in the isogeometric setting (Bontinck et al., 2018; Kapidani et al., 2022). Another approach to realize rotation is to replace the air-gap region by another field description, e.g. the boundary element method (Kurz et al., 1999) or an analytical formulation, i.e. the air-gap method (Abdel-Razek et al., 1982).

Isogeometric analysis (IGA) is a well-suited FEM variant for simulating electric machines, particularly due to its ability to represent geometries exactly. Introduced by Hughes et al. (Hughes et al., 2005), IGA bridges the gap between computer-aided design (CAD) and finite element analysis by using CAD basis functions, i.e. B-splines and nonuniform rational B-splines, for both geometry representation and the basis functions for the simulation. It allows for an exact geometric representation of cylindrical structures and the smoothness of IGA basis functions can improve the accuracy of simulations for the (smooth) magnetic fields in electric machine analysis.

The combination of the air-gap method with IGA is especially well suited as both methods represent the rotor and stator interfaces of the air gap geometrically exactly. This paper introduces, for the first time, a method that integrates a flexible rotor movement enabled by a scaled air-gap element in the IGA context, combining those approaches to overcome these difficulties and enhance simulation efficiency. However, a naive implementation of the air-gap element leads to badly conditioned algebraic systems. Therefore, we introduce a new scaling which is inspired by (Gangl et al., 2022).

The paper is structured as follows. Sections 2 and 3 introduce the model and in particular the air-gap element. Sections 4 and 5 discuss the discretization in space and time. The quantity of interest (torque) is defined in Section 6. A machine model is introduced in Section 7 and the simulation results are shown in Section 8. The paper closes with conclusions in Section 9.

Two-dimensional simulations of electric machines, e.g. (Salon, 1995), deal with complex variants of the principal geometry Ω depicted in Figure 1, where the interior domain Ω_R is the rotor (e.g. including permanent magnets), Ω_S denotes the region of the stator (e.g. containing slots with the windings) and Ω_A denotes the region of the air gap. We use the isogeometric variant of the FEM (IGA) (Hughes et al., 2005), i.e. those domains are exactly represented by patches of spline-based mappings.

An appropriate low-frequency model of Maxwell’s equations (Jackson, 1998) is the eddy current approximation (Schmidt et al., 2008), which reads in two dimensions:

with conductivity σ, reluctivity ν and the z-components of magnetic vector potential A, excitation current density J and suitable boundary conditions. As there is no excitation current density in the air gap, we split the domain Ω into nonoverlapping subdomains, yielding:

with the constraint that A_z and its conormal derivative are continuous across the interface boundaries Γ₁ and Γ₂. If we denote the solutions in the rotor, air and stator region by $AzR,AzA,AzS$⁠, respectively, the transmission conditions are:

where the subscript of ν indicates the possibly different material behavior of the respective regions toward the interface. If the rotor moves, the transmission conditions are angle dependent. Note that the normal vectors n_R, n_S and n_A are defined exterior to the respective domains.

For the solution in the air-gap region, we follow the ideas of Abdel-Razek et al. (1982). Due to the simple geometry, the air-gap solution of equation (4) can be computed analytically. In polar coordinates, we have:

where the coefficients $α0,α′0,αk,α′k,βk,β′k$ are to be determined. As there is no time derivative in equation (4), and hence the solution is only dependent on boundary conditions at every point in time; for better readability, we suppress the time dependency of the coefficients in this section. We impose Dirichlet boundary conditions:

on both interfaces, where ν_i is assumed to be given. For comparison of coefficients, the boundary conditions are evolved into a Fourier series, i.e.

with:

where R_i denotes the radius of the interface and for brevity of notation, we make the values of the boundary data dependent on the arc length on the respective interface. With the Fourier coefficients (10)–(12) of the boundary data, enforcing Dirichlet conditions for the zeroth mode leads to the 2 × 2 system:

and for all higher modes, we have to solve the two independent 2 × 2 systems:

and:

If one interface is rotated by an angle δ, the calculation of the boundary data’s Fourier coefficients changes and one can use properties of trigonometric functions to obtain:

and analogously:

from the data ν_i in the reference domain. So, the rotated coefficients are just weighted sums of the Fourier coefficients in the reference domain. Noticing, that the integrals on the right hand side are the standard Fourier coefficients, this can be put into matrix form yielding:

where the occurring matrix is a rotation matrix for the angle kδ.

We consider the remaining equation (3) on the reference domain with angle-dependent transmission conditions. As in this viewpoint, the observer is fixed to the respective domain, the total time derivative becomes a partial one. Using a variational approach in space and applying integration by parts, yields the problem to find $Az∈HΓ1(ΩD)$⁠, such that:

for all $V∈HΓ1(ΩD)$⁠, where Ω_D = Ω_R ∪ Ω_S and $HΓ1(ΩD)$ denotes H¹-functions with vanishing trace on Γ. Further S is the rotation-dependent Dirichlet-to-Neumann mapping, which assigns to every function in $HΓ1(ΩD)$ and angle δ the according Neumann trace arising from the air-gap solution. We discretize equation (19) by a Galerkin approach, i.e.:

where ${ψn}n=1N$ is the basis of the IGA space comprising rotor and stator, yielding:

∀j = 1, …, N. We collect all IGA degrees of freedom in a vector A. The discretized Dirichlet-to-Neumann mapping S_h(δ) is realized via the air-gap solution. For that purpose, we abort the infinite series in equation (7) at a finite number K. The necessary integrals for the Fourier coefficients are computed by:

If we now collect all Fourier coefficients ν_i,k and $v′i,k$ for both interfaces into a vector F, then there is a linear relation between F and A, i.e.

where the row $r¯⊤$ of C₂ corresponding to equation (24) is given by:

The other rows follow similarly. Note that due to equation (18), the associated matrix can be efficiently adapted to different rotation angles and quadrature is only needed once for the reference domain, yielding the angle-dependent matrix C₂(δ). The evaluation of ∂_nA_z at the interfaces, taking the air-gap solution to determine the Neumann trace, yields similar integrals and a matrix denoted by C₁(δ). Note that these matrices have the same structure as the coupling matrices in Egger et al. (2022). When assembling the degrees of freedom of the IGA solution A together with the air-gap solution coefficients $α¯=(α0,α′0,…,β1,β′1…)⊤$ into a vector:

we obtain the following block system of algebraic equations:

which is a differential algebraic system in time. M and K are the standard mass and stiffness matrices of the IGA discretization, J is the standard load vector, K_F comprises the 2 × 2 equation systems for the air-gap solution. C₁(δ) and C₂(δ) are defined above. Note that the discrete Dirichlet-to-Neumann mapping:

is implicitly contained in the differential-algebraic system of equations. For brevity, we define:

For the computation of the coefficients of the analytical solution (7), one has to solve 2 × 2 systems [equations (13)–(15)], where the occurring matrices only depend on R₁, R² and the mode number k. In the later used machine model, the radii are R₁ = 44.3 × 10^–3 m and R² = 44.7 × 10^–3 m. For these values, the condition number κ₂ of the resulting 2 × 2 systems grows prohibitively as depicted in Figure 2(a). Therefore, we choose to evaluate the air-gap element on a scaled domain with R₁ = 1 and use the scaling laws for the magnetic fields to meet the coupling conditions (Stipetic et al., 2015). This leads to the scaled 2 × 2 systems:

with $R˜=R2R1$⁠. As shown in Figure 2(b), the condition number remains reasonable, even for higher mode numbers.

In addition to standard time-stepping schemes, we solve equation (28) using a variational approach in time. For this purpose, let I = (0, T) be the time interval of interest, which is divided into N_t (not necessarily equidistant) intervals, yielding a one-dimensional mesh T_h on which we define piecewise linear and continuous ansatz functions ${φi(t)}i=0Nt$ and piecewise constant test functions ${ψj(t)}j=1Nt$⁠. The solution in space-time is then given by the ansatz:

Inserting this into equation (28) and testing with ψ_j yields:

This is a discontinuous Galerkin formulation in time, introduced e.g. in Gander and Neumüller (2016). If we use midpoint rule for quadrature and denote the mesh width of the i-th element by h_i and its midpoint by t_i, this leads to the block equation system:

with

Note that the block matrix (33) is rectangular having more columns than rows. Given an initial condition U₀, the homogenized system is given by:

Note that this system is sub-diagonal and hence can be solved via forward substitution. This would correspond to a time-stepping scheme. However, the interpretation as a space-time system opens possibility for parallel solvers in the temporal dimension, e.g. via multigrid (Gander and Neumüller, 2016).

In electric machines, one of the key performance indicators is the electromagnetic torque, which describes the rotational force generated by the machine. The torque acting on a volume can be determined by integrating the Maxwell stress tensor over a surface that encloses this volume (Salon, 1995). For electric machines in a 2D context, the torque is given by:

with the machine length L, the vacuum permeability μ₀ and the radius of the integration surface/line r. In conventional machine simulations using low-order finite elements, torque evaluation often suffers from inaccuracies. These issues primarily arise from the numerical differentiation of A_z to obtain the required components of B, which amplifies numerical errors in the vector potential. The angular component B_φ is especially prone to inaccuracies because its computation involves the normal derivative of A_z across element boundaries. To reduce these errors, various techniques have been proposed, such as the eggshell method or Arkkio’s method, which use a volume integral instead of a surface integral (or a surface integral in place of a line integral in 2D) (Arkkio, 1987; Henrotte et al., 2004; Henrotte et al., 2010). Using the air-gap method resolves these accuracy challenges, as the solution in the air gap is known in closed form, facilitating the derivation of the magnetic vector potential and eliminating errors from numerical differentiation. The B field can then be calculated from A_z via:

where the matrix entries are with respect to Cartesian coordinates. The gradient in polar coordinates is given by:

and straightforward calculations yield the components of the B field in polar coordinates:

For any radius r inside the air gap (Abdel-Razek et al., 1981), this yields:

Using orthogonality of the occurring functions and exploiting:

yields

To validate the described approach, we consider the permanent magnet synchronous motor given in Figure 3. The model consists of a three-phase driven stator with homogenized copper coils in red, the linear iron yoke in gray and air shown in blue. The rotor comprises a permanent magnet in green with additional air slits to compensate for nonlinear saturation effects. This justifies to consider only linear material behavior in the following. The air-gap element for the coupling is shown as orange band that connects rotor and stator. Due to symmetry, only one-sixth of the motor is simulated with appropriate boundary conditions. Eddy currents are assumed to be fully prevalent in the magnet, and in a reduced way in the iron due to lamination which is modeled using the formulation from Gyselinck et al. (1999).

The current density is given by the three-phase current $Jz=J0∑k=13sinωt+2/3πk$ with the pole-pair number p = 3 and ω = 2π⋅50 Hz. For synchronous operation, the rotation angle $δ=ωtp$ is prescribed by a constant rotation velocity.

The rotor and stator domains are discretized using quadratic B-Splines with 6044 and 5196 degrees of freedom, respectively. The analytical solution of the air-gap element is expressed using the first 35 nonzero terms in equation (7), resulting in 140 unknowns for the coupling. Since only one-sixth of the motor is simulated and the solution is antiperiodic, this corresponds to the modes k ∈ {3, 9, 15, …, 207}.

We perform transient simulations on the presented machine for the first 20° rotational degrees. After that, the torque will be periodic in the stationary case. The time interval is discretized with N_t = 40 elements. This corresponds to 40 time-steps of a classical time-stepping scheme. At this point, the resulting block matrix system is solved using a direct solver. The initial conditions for the problem are computed by a static simulation at time zero. The static and transient simulation results are shown in Figure 4. The torque values from the static solution are given as reference. As expected, the transient torque is slightly reduced because of the energy losses due to eddy-currents in the magnet and iron which are only captured in the transient case. The solution from the space-time setting is in very good agreement with the solutions obtained from the time-stepping scheme.

In this contribution, we coupled the analytical solution of the air-gap element to numerical solutions discretized with IGA and solved the transient problem with a space-time formulation. The combination of these methods is advantageous for several reasons: First, the boundary integrals are evaluated precisely due to the exact geometry representation. Second, the rotor movement can be incorporated without the need of remeshing or reevaluation of boundary integrals. Third, the discretization of the time-domain can be carried out a priori even without known rotor movement. Fourth, the space-time setting allows for new possibilities regarding the parallel solving of the equation system, e.g. using Parareal (Schöps et al., 2018; Bast et al., 2020), and multigrid techniques (Bolten et al., 2020).

Future steps include the extension of the simulations to the nonlinear case, i.e. nonlinear material behavior and unknown rotor motion. Finally, benchmarking must show if the space-time approach has computational benefits.

The work is supported by the joint DFG/FWF Collaborative Research Centre CREATOR (CRC – TRR361/10.55776/F90) at TU Darmstadt, TU Graz and JKU Linz.

Funding: This study received funding from Austrian Science Fund, No. TRR361/F90; Deutsche Forschungsgemeinschaft, No. TRR361/F90.