How to define a simulation to call FEMM

This tutorial shows the different steps to compute magnetic flux and electromagnetic torque with pyleecan automated coupling with FEMM.

The notebook related to this tutorial is available on GitHub.

Every electrical machine defined in Pyleecan can be automatically drawn in FEMM to compute torque, airgap flux and electromotive force. To do so, the tutorial is divided into four parts:
- defining or loading the machine
- defining the simulation inputs
- setting up and running of the magnetic solver
- plotting of the magnetic flux for the first time step

Defining or loading the machine

The first step is to define the machine to simulate. For this tutorial we use the Toyota Prius 2004 machine defined in this tutorial.

# Change of directory to have pyleecan in the path
from os import chdir

from pyleecan.Functions.load import load

# Import the machine from a script
IPMSM_A = load('pyleecan/Data/Machine/IPMSM_A.json')

# Plot the machine
%matplotlib notebook

Simulation definition


The simulation is defined with a Simu1 object. This object correspond to a simulation with 5 sequential physics (or modules): - electrical
- magnetic
- force
- structural
- acoustic

Simu1 object enforce a weak coupling between each physics: the input of each physic is the output of the previous one.

In this tutorial we will focus only on the magnetic module. The Magnetic physic is defined with the object MagFEMM and the other physics are desactivated (set to None).

We define the starting point of the simulation with an InputCurrent object to enforce the electrical module output with: - angular and the time discretization
- rotor speed
- stator currents
from numpy import ones, pi, array, linspace
from pyleecan.Classes.Simu1 import Simu1
from pyleecan.Classes.InputCurrent import InputCurrent
from pyleecan.Classes.MagFEMM import MagFEMM

rotor_speed = 2000 # [rpm]

# Create the Simulation
mySimu = Simu1(name="EM_SIPMSM_AL_001", machine=IPMSM_A)

# Defining Simulation Input
mySimu.input = InputCurrent()

# time discretization [s]
mySimu.input.time.value= linspace(start=0, stop=60/rotor_speed, num=16, endpoint=False)# 16 timesteps

# Angular discretization along the airgap circonference for flux density calculation
mySimu.input.angle.value = linspace(start = 0, stop = 2*pi, num=2048, endpoint=False) # 2048 steps

# Rotor speed as a function of time [rpm]
mySimu.input.Nr.value = ones(16) * rotor_speed

# Stator currents as a function of time, each column correspond to one phase [A]
mySimu.input.Is.value = array(
        [ 1.77000000e+02, -8.85000000e+01, -8.85000000e+01],
        [ 5.01400192e-14, -1.53286496e+02,  1.53286496e+02],
        [-1.77000000e+02,  8.85000000e+01,  8.85000000e+01],
        [-3.25143725e-14,  1.53286496e+02, -1.53286496e+02],
        [ 1.77000000e+02, -8.85000000e+01, -8.85000000e+01],
        [ 2.11398201e-13, -1.53286496e+02,  1.53286496e+02],
        [-1.77000000e+02,  8.85000000e+01,  8.85000000e+01],
        [-3.90282030e-13,  1.53286496e+02, -1.53286496e+02],
        [ 1.77000000e+02, -8.85000000e+01, -8.85000000e+01],
        [ 9.75431176e-14, -1.53286496e+02,  1.53286496e+02],
        [-1.77000000e+02,  8.85000000e+01,  8.85000000e+01],
        [-4.33634526e-13,  1.53286496e+02, -1.53286496e+02],
        [ 1.77000000e+02, -8.85000000e+01, -8.85000000e+01],
        [ 4.55310775e-13, -1.53286496e+02,  1.53286496e+02],
        [-1.77000000e+02,  8.85000000e+01,  8.85000000e+01],
        [-4.76987023e-13,  1.53286496e+02, -1.53286496e+02]

MagFEMM configuration

For the configuration of the Magnetic module, we use the object MagFEMM that compute the airgap flux density by calling FEMM. The model parameters are set though the properties of the MagFEMM object. In this tutorial we will present the main ones, the complete list is available by looking at Magnetics and MagFEMM classes documentation.

type_BH_stator and type_BH_rotor enable to select how to model the B(H) curve of the laminations in FEMM. The material parameter and in particular the B(H) curve are setup directly in the machine.

from pyleecan.Classes.MagFEMM import MagFEMM
# Definition of the magnetic simulation (is_mmfr=False => no flux from the magnets)
mySimu.mag = MagFEMM(
    type_BH_stator=0, # 0 to use the B(H) curve,
                           # 1 to use linear B(H) curve according to mur_lin,
                           # 2 to enforce infinite permeability (mur_lin =100000)
    type_BH_rotor=0,  # 0 to use the B(H) curve,
                           # 1 to use linear B(H) curve according to mur_lin,
                           # 2 to enforce infinite permeability (mur_lin =100000)
    angle_stator=0,  # Angular position shift of the stator
    file_name = "", # Name of the file to save the FEMM model

mySimu.struct = None # We only use the magnetic part

Pyleecan coupling with FEMM enables to define the machine with symmetry and with sliding band to optimize the computation time. Here we will draw and mesh only 1/8 of the machine (4 symmetry + antiperiodicity):

mySimu.mag.is_symmetry_a=True   # 0 Compute on the complete machine, 1 compute according to sym_a and is_antiper_a
mySimu.mag.sym_a = 4 # Number of symmetry for the angle vector
mySimu.mag.is_antiper_a=True # To add an antiperiodicity to the angle vector

At the end of the simulation, the mesh and the solution can be saved in the Output object with:

mySimu.mag.is_get_mesh = True # To get FEA mesh for latter post-procesing
mySimu.mag.is_save_FEA = False # To save FEA results in a dat file

Run simulation

To run the simulation, we first have to set the Output to store the results.

from pyleecan.Classes.Output import Output
myResults = Output(simu=mySimu)

When running the simulation, a FEMM window should open so you can see pyleecan drawing the machine and defining the surfaces. image0 The simulation will compute 16 different timesteps by updating the current and the sliding band boundary condition.

Once the simulation is finished, the results are stored in the magnetic part of the output (i.e. myResults.mag ) and we can call different plots. This object contains:
- time: magnetic time vector without symmetry
- angle: magnetic position vector without symmetry
- Br: radial airgap flux density
- Bt: tangential airgap flux density
- Tem: electromagnetic torque
- Tem_av: average electromagnetic torque
- Tem_rip: torque ripple
- Phi_wind_stator: stator winding flux
- emf: electromotive force

Plot results

Output object embbed different plot to visualize results easily. You can find a dedicated tutorial here.

For instance, we can plot the radial and tangential magnetic flux in the airgap at a specific timestep. We also plot their space Fourier Transform by setting is_fft to True.

%matplotlib notebook
myResults.plot_A_space("mag.Br",t_index=1, is_fft=True, r_max=76)
%matplotlib notebook
myResults.plot_A_space("mag.Bt",t_index=1, is_fft=True, r_max=76)

You can also define your own plot. The following plot requires plotly to display the radial flux density in the airgap over time and angle.

#%run -m pip install plotly # Uncomment this line to install plotly
import plotly.graph_objects as go
from plotly.offline import init_notebook_mode

x = myResults.mag.angle*180/pi # rad -> °
y = myResults.mag.time
z = myResults.mag.Br.values
fig = go.Figure(data=[go.Surface(z=z, x=x, y=y)])
fig.update_layout( )
fig.update_layout(title='Radial flux density in the airgap over time and angle',
                  scene = dict(
                      xaxis_title='Angle [°]',
                      yaxis_title='Time [s]',
                      zaxis_title='Flux [T]'
                  margin=dict(r=20, b=100, l=10, t=100),
                 ) = {"displaylogo":False})