Interactive and educational simulation tool for permanent magnet synchronous machines

Abstract

Purpose

This paper aims to present an interactive design and simulation tool for permanent magnet synchronous machines based on the finite-element-method. The tool is intended for education and research on electrical machines.

Design/methodology/approach

A coupling between the software MATLAB and finite element method magnetics is used. Several functionalities are included as modular scripts and represented in the form of a graphical user interface. Included are fully parametrized motor models, automatic winding generations and the evaluation of torque waveforms, core losses and speed-torque-diagrams. A survey was conducted to determine how the motivation of students concerning the covered topics is influenced by using the tool.

Findings

Due to its simplicity and the intuitive visualization of the results, the tool provides direct access to the topic of electrical machines without having to deal with separate scripts. The modular structure of the software allows simple extensions with new functions. Because students can directly contribute to the tool with their own work, their motivation for using and extending it increases.

Originality/value

The presented tool offers more functionalities compared to similar free software packages, e.g. the calculation of core losses and speed-torque diagrams. Also, it is designed in such a way that it can be easily understood and extended by students.

Citation

Schugardt, A., Kaiser, L., Avcilar, F. and Schäfer, U. (2024), "Interactive and educational simulation tool for permanent magnet synchronous machines", COMPEL - The international journal for computation and mathematics in electrical and electronic engineering, Vol. ahead-of-print No. ahead-of-print. https://doi.org/10.1108/COMPEL-11-2023-0554

1. Introduction

Nowadays, the design and optimization of permanent magnet synchronous machines (PMSM) is done mainly with software based on the finite-element-method (FEM). One of these is finite element method magnetics (FEMM) (Meeker, 2019), which is a free and open-source alternative to commercially available programs but lacks a comprehensive user interface and parametrized machine models. The connection to Octave, MATLAB or Python through predefined Lua script commands offers a wide range of possibilities for the design of electrical machines. The presented simulation tool is based on MATLAB functions that are coupled with FEMM and structured regarding their corresponding tasks. Similar programs (Bonneel et al., 2018; de Andrade et al., 2019; Kuptsov et al., 2019) and libraries (Lehikoinen et al., 2018) have been created before, but they either do not offer as many options as the presented tool or are not as accessible for beginners. The machine geometry is fully parametrized using MATLAB scripts. Furthermore, the winding table is automatically generated, and a series of magnetostatic simulations can be executed to create (discrete) time-dependent solutions. All results are evaluated and displayed.

A graphical user interface (GUI) based on a MATLAB application was created to support and simplify the usage of the scripts. First, this allows easier access to the topic and introduces students to the benefits of field simulations. The calculation methods for obtaining all needed quantities are introduced as part of lectures. Subsequently, the modular structure of the scripts allows a wide range of extensions, for example, in the context of student projects. The learning objective is to acquire knowledge about the calculation of e.g. counter-electromotive force (back-EMF), torque and core losses based on the vector potentials and flux densities from the FEM solution.

2. Simulation tool

To simulate a PMSM using the interactive GUI, the user must specify the geometry parameters. Furthermore, information about materials (conductors, iron core and magnets), current waveform and speed is required. The geometry can be plotted any time inside the GUI. Thus, changes in geometry parameters can be seen instantly.

2.1 Geometry and materials

The stator geometry is fully parametrized. It is possible to change the stator type from parallel slots to parallel tooths. This is useful for the simulation of rectangular conductors (e.g. in hairpin windings). Furthermore, details of the slot or tooth head geometry can be changed (e.g. the radius of an arc or the form of the tooth head).

The rotor geometry is also fully parametrized. Currently, the following topologies can be used: surface-mounted magnets, interior magnets with rectangular shapes (as shown in Figure 1) and V-shaped magnets. Switched reluctance machines and synchronous reluctance machines are included as well. It is also possible to add more parametrized rotor topologies as independent MATLAB functions.

In the next step, the materials can be selected. These include conductor material, magnet material and electrical steel. The unmodified library already offers a wide range of materials. Custom materials can be added in a comfortable and intuitive way inside the GUI. A spline interpolation is used for nonlinear B(H) characteristics of e.g. electrical steel. To avoid deviations from the given data, enough points for especially low and high flux densities must be specified.

2.2 Winding and current

The winding design is automatically performed based on the user input for the number of slots N, pole pair number p and winding layers by application of a star of slots. The symmetry conditions for m phases and t = gcd⁡(N, p) (greatest common divisor) are as follows:

(1) one-layer-winding: Nt⋅m∈ℕ  and  N2m∈ℕ

(2) two-layer-winding: Nt⋅m∈ℕ  and  Nm∈ℕ

The positions of the positive and negative phases inside the slots are determined based on the phasor diagram of the slot voltages. To demonstrate the concept, a two-layer-winding with N = 9, p = 4 and m = 3 is going to be analyzed. The phasors are drawn with a phase shift of Δφ = 360°⋅p/N = 160° to one another (left-hand-side of Figure 2). Starting with slot 1, the positive phases U+ (current out of the plane) are put into all slots over an angle of 360°/(2m) = 60°. Then, the next positive phase V+ starts after 360°/m = 120°. The negative phases are put into the opposite slots of the phasor diagram. The result is shown in Figure 2 on the right (Pyrhönen et al., 2013).

Based on this phasor diagram, a generated winding table shows the position of all positive and negative phases with their corresponding slot numbers (see Figure 3 at the bottom). Furthermore, the number of slots per pole per phase q and the winding factors are calculated using the following formula (Müller et al., 2008):

(3) slots per pole per phase: q=N2⋅p⋅m=qZqN

(4) distribution factor: ξZ,k=sin⁡(k⋅π2m)qZ⋅sin⁡(k⋅π2m⋅qZ)

(5) chording factor: ξS,k=sin⁡(k⋅π2⋅sτP)

The winding factor is calculated by ξ_k = ξ_Z,k⋅ξ_S,k for the spatial harmonic k. The pole pitch equals τP=π⋅di,Sta2p using the inner diameter of the stator d_i,Sta. The skewing s is derived from the winding table. It is the distance between a positive phase and its negative counterpart as an arc length. The winding factor is plotted as a function of the harmonic orders (see Figure 3). Also, it is automatically displayed whether one- or two-layer windings are possible (Burress et al., 2008; Hsu et al., 2004).

2.3 Meshing and field simulation

For the field simulation, a partial differential equation (PDE) is needed. At first, Maxwell’s equations and the material law (Pyrhönen et al., 2013) are considered for the magnetostatic case (time t = const. and the electrical field is neglected):

(6) rot H→=J→

(7) div B→=0

(8) B→=μ⋅H→

H→ is the magnetic field strength, B→ the flux density, J→ the current density and μ the permeability. The magnetic vector potential A→ is defined as B→=rot A→. The derivation of the PDE is based on the given equations, and the Coulomb gauge div A→=0 is considered (Müller et al., 2008):

(9) rot H→=rotB→μ=1μrot rot A→=1μ⋅(graddiv A→︸= 0−div grad A→︸Δ)=J→

(10) ↔1μ⋅ΔA→=−J→ 

The second step assumes that the permeability does not depend on the spatial coordinates, and Δ is the Laplace operator. A formulation for the 2D case is enough to model the basic aspects of radial flux machines. Therefore, the PDE can be simplified for an infinitely extended z-axis:

(11) 1μ⋅(∂2Az∂x2+∂2Az∂y2)=Jz

The FEM is used to solve the PDE. At first, the 2D geometry of the machine is meshed in FEMM using first-order triangular elements with three nodes (see Figure 4). The meshing algorithm is based on the Delaunay method (Meeker, 2019).

The z-component of the vector potential Aze(x,y) inside an element e is calculated using the area of the element Ω^e and interpolation constants a_n, b_n and c_n. The interpolation function N_n(x, y) is 1 at the position of node n and 0 at the position of the other nodes (Polycarpou, 2006). This results in:

(12) Nn(x,y)=12Ωe⋅(an+bn⋅x+cn⋅y)

(13) Aze(x,y)=A1,z⋅N1(x,y)+A2,z⋅N2(x,y)+A3,z⋅N3(x,y)

Therefore, the vector potentials A_1,z, A_2,z and A_3,z of the three nodes must be calculated first. Using the method of weighted residuals and partial integration, the given 2D PDE [equation (11)] is converted to the FEM matrix equation. Because of the nonlinear permeability, the equation is solved using the fixed-point algorithm or the Newton–Raphson method (Aliprantis and Wasynczuk, 2023). The solution contains the vector potentials of all nodes inside the given mesh. The three corresponding vector potentials of an element e (see Figure 4) are used to calculate the flux density:

(14) B→e=rot A→e=rot(00Aze)=(∂Aze∂y−∂Aze∂x0)=12Ωe(A1,z⋅c1+A2,z⋅c2+A3,z⋅c3−(A1,z⋅b1+A2,z⋅b2+A3,z⋅b3)0)

Afterward, the flux densities are smoothed using the nearest neighbor algorithm. This results in flux densities that are not constant inside an element (Meeker, 2019; Polycarpou, 2006).

2.4 Evaluation and results

The flux densities are used to calculate different quantities and waveforms of the machine. A magnetic flux Φ in one part of the machine is derived from the integration of the flux density over the relevant cross-section area. All flux that goes through the coils of one phase contributes to the corresponding flux linkage. For phase U the flux linkage is as follows:

(15) ΨU=w⋅ΦU

Here, w is the number of turns. The flux linkage is evaluated for different mechanical rotation angles φ_mech of the rotor. Because of the periodicity of machines with p > 1, the waveforms are plotted over a full electrical period with φ_el = p⋅φ_mech. Note, that these (temporal) angles can be converted to corresponding time values based on the electrical frequency f₁ of the machine. An example solution for all three phases can be seen in Figure 5 at the bottom left.

Next, the back-EMF u_i of a phase can be derived based on Faraday’s law of induction (Pyrhönen et al., 2013):

(16) ui=−dΨdt=−2π⋅n⋅p⋅dΨdφel

The speed of the rotor n is used. Figure 5 shows example waveforms at the upper-right.

To calculate the torque of the machine, Maxwell’s stress tensor for the two-dimensional case in cylindrical coordinates (r, φ, z) is used:

(17) Tm=12μ⋅(Br2−Bφ22⋅Br⋅Bφ2⋅Br⋅BφBφ2−Br2)

This term is evaluated for every element inside the airgap. Then, the force vector can be calculated by integration over the rotor volume V_r or the rotor surface A_r and using the outer radius of the rotor r_o,r (Binder, 2012):

(18) F→=(FrFφ)=∭VrdivTm⋅dV=∯∂VrTm⋅dA→r=∬Ar12μ(Br2−Bφ22⋅Br⋅Bφ)⋅ro,r dφ dz

The torque has only a z-component that can be derived from the given force:

(19) T→=r→×F→=(ro,r00)×(FrFφ0)=(00ro,r⋅Fφ)

The torque waveform for the example geometry is shown in Figure 5 at the bottomright. The maximum electrical angle for the simulation and the overall number of points per simulation are input values inside the GUI. A preview plot shows the expected waveforms of the current with one period and the torque with ν periods over 360° electrically, where:

(20) ν = lcm(2p, N)p

(Müller et al., 2008; Pyrhönen et al., 2013).

A core loss calculation based on loss separation (Bertotti, 1988; Zhou and Bowen, 2020) is implemented. The user can specify the regression coefficients C_h, C_e, C_ex, n_h and n_ex for the used material data. The flux density is evaluated in every mesh element over the simulation steps. Then, the flux density waves are decomposed into their harmonic content. The specific losses are calculated based on the harmonic flux density amplitude of every element and the corresponding harmonic frequency f_k:

(21) pc,ke= Ch⋅fk⋅(Be)nh+ Ce⋅fk2⋅(Be)2+Cex⋅(fk⋅Be)nex

For every element, the harmonic losses are summed. Afterward, a summation of all element losses multiplicated with the element volume V^e is performed:

(22) Pc=∑e((∑kpc,ke)⋅Ve)

2.5 Coordinate transformation

It is possible to change the amplitude and phase shift for the three-phase current. Alternatively, d- and q-currents can be set directly. The following transformation formula (Park, 1929) is used:

(23) (idiq)=23(cos⁡(φ)cos⁡(φ−2π3)cos⁡(φ−4π3)−sin⁡(φ)−sin⁡(φ−2π3)−sin⁡(φ−4π3))(iUiViW)

Therefore, the initial rotor position as a reference must be found. At first, the normal airgap flux density is calculated for a model without permanent magnets (material replaced with air) over the spatial angle. The phase shift of the fundamental component φrot′ is evaluated using a Fourier decomposition of the waveform. The corresponding plots can be seen in Figure 6.

The mechanical rotation angle for the initial position can then be calculated using the number of pole pairs:

(24) φrot=90°−φrot′p

An angle of 90° is used to put the initial position on the q-axis. Figure 7 shows an example of pure d- or q-currents without permanent magnets. The pure q-current is generated with a current phase shift of 0°. A phase shift of – 90° is needed for a pure positive d-current.

Then, the flux linkages and inductances in dq-coordinates can be calculated. The values can also be plotted over different d- and q-current pairs. Starting from this point, a torque-speed-diagram can be generated.

3. Tool validation

The electrical motor from the Toyota Prius 2004 is used as a reference machine to validate the functionalities of the presented tool. The necessary data are taken from Burress et al. (2008) and Hsu et al. (2004). Figure 8 shows the motor geometry that is used inside of FEMM. Due to its symmetry, 1/8 of the machine is enough to characterize the magnetic field. Therefore, symmetry conditions must be defined. Here, only one pole is modeled, so an antiperiodic condition is needed. Dirichlet boundaries (A = 0) are used for the inner and outer radius of the motor.

After simulation, the results are compared to measured locked-rotor torque curves. The plot for a peak current of 250A is displayed in Figure 9. It can be shown that the results obtained using the tool qualitatively agree with the comparison values. In consequence, they can be considered valid. Quantitative deviations can be attributed to unknown measurement conditions and inaccuracies. The procedure for obtaining the locked-rotor simulation results is described in Katona and Orosz (2022).

4. Educational purposes

The presented tool offers a variety of possible applications for educational purposes in the field of electrical machines, which is supported by its modular and transparent structure. This allows different levels of detail that can be used for teaching (see Figure 10) to achieve a better understanding of electrical machines and the underlying theory (Stoev et al., 2017) and thus educational success:

• Level 1: Using the GUI to simulate and design PMSM. This includes, for example, the usage of parameter studies and to learn about winding design.

• Level 2: Looking inside the underlying MATLAB scripts. The students learn how to run FEMM and how the solution can be used to calculate the back-EMF, torque, dq-inductances, torque-speed-diagrams and efficiency maps.

• Level 3: Learning about the basic principles of FEM. It can be taught how FEMM calculates all field quantities based on the solution of a PDE. This knowledge can also be applied to other disciplines, like mechanical or thermal simulations.

The desired degree of detail can be chosen together with the students. This provides an opportunity to adapt the teaching to the interests of the students. In addition, different types of tasks for lectures, student projects or theses are possible. This includes a complete motor design using the GUI, an extension of the MATLAB scripts with new functionalities or the modeling of real motors and the comparison of simulation results with measurement data. In contrast to many other projects, students who are involved in expanding the GUI can be certain that their results will find application over the long term, for example, in the context of lectures and research tasks. Overall, this can contribute to greater motivation among the students. Several theses have already been assigned to extend the simulation tool, e.g. the addition of parametrized hairpin models. A survey was conducted with students after using the program for motor simulations. The following three statements were evaluated:

• statement 1: “My motivation for the topics is increased by using the program”;

• statement 2: “I would like to use the program in more lectures or projects”; and

• statement 3: “I would like to learn more about how the program works”.

The answers are based on a scale from −2 (= totally disagree) to 2 (= totally agree). Figure 11 shows the results for a sample size of nine students. The mean values are marked with red lines.

Due to the small number of participants, the survey is not representative. Because the participants were all electrical engineering students with a clear interest in the subject of electrical machines, it is still possible to derive trends from the results. Looking at the answers to the first statement, it is noticeable that the motivation of most participants in relation to the topic has increased using the program. For the second statement, it is remarkable that all students surveyed would like to use the software in future courses. Regarding the answers to the third statement, some students are very interested in learning how exactly the program works, others are not particularly concerned. Overall, using the tool in the context of teaching is very well received by the students surveyed, but not everyone necessarily wants to understand exactly how the underlying functions and calculations work.

5. Summary and outlook

This paper provides a comprehensive overview of the FEMM GUI, an easy-to-use tool for the simulation of electrical machines, which was developed primarily to be used in the context of education. Although a large variety of functions and models are already included, there is still much room for improvement. The full support of electrically excited synchronous machines and induction machines could be a possible extension in the future. Also, thermal simulations in FEMM and mechanical simulations using other free software could be included and linked to the already existing models. The presented tool will continue to offer students the opportunity to incorporate their own ideas.

Figures

Figure 1.

Rotor parametrization with interior magnets and hairpin winding

Figure 2.

Phasor diagram of slot voltages with positive and negative phases

Figure 3.

Winding tab of the GUI (data from Toyota Prius 2004)

Figure 4.

First-order triangular element

Figure 5.

Result plots for the Toyota Prius 2004 with 125A at 1,540 rpm

Figure 6.

Airgap flux density plot and harmonic waves (fundamental, 5. and 7.)

Figure 7.

Magnetic flux lines for a pure q-current (left) and a pure d-current (right)

Figure 8.

Toyota Prius 2004 motor model (left: gray = electrical steel, red = magnets | right: model in FEMM)

Figure 9.

Locked rotor torque of the Prius 2004 for 250A

Figure 10.

Educational levels of detail using the presented tool

Figure 11.

Results of the student survey (scale: −2 = totally disagree, −1 = rather disagree, 0 = neutral, 1 = rather agree and 2 = totally agree)