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. This tutorial was tested with the release 21Apr2019 of FEMM. Please note that the coupling with FEMM is only available on Windows.

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.

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.

[1]:
%matplotlib notebook
# Load the machine
from os.path import join
from pyleecan.Functions.load import load
from pyleecan.definitions import DATA_DIR

IPMSM_A = load(join(DATA_DIR, "Machine", "IPMSM_A.json"))
IPMSM_A.plot()

Simulation definition

Inputs

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

For now pyleecan includes: - an Electrical module for PMSM machine with FEMM - a Magnetic module with FEMM for all machines - a Force module (Maxwell Tensor)

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

The Magnetic physics is defined with the object MagFEMM and the other physics are deactivated (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
[2]:
from os.path import join

from numpy import ones, pi, array, linspace, cos, sqrt

from pyleecan.Classes.Simu1 import Simu1
from pyleecan.Classes.InputCurrent import InputCurrent
from pyleecan.Classes.MagFEMM import MagFEMM

# Create the Simulation
simu_femm = Simu1(name="Webinar_1_MagFemm", machine=IPMSM_A)
p = simu_femm.machine.stator.winding.p
qs = simu_femm.machine.stator.winding.qs

# Defining Simulation Input
simu_femm.input = InputCurrent()

# Rotor speed [rpm]
simu_femm.input.N0 = 2000

# time discretization [s]
time = linspace(start=0, stop=60/simu_femm.input.N0, num=32*p, endpoint=False) # 32*p timesteps
simu_femm.input.time = time

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

# Stator currents as a function of time, each column correspond to one phase [A]
I0_rms = 250/sqrt(2)
felec = p * simu_femm.input.N0 /60 # [Hz]
rot_dir = simu_femm.machine.stator.comp_rot_dir()
Phi0 = 140*pi/180  # Maximum Torque Per Amp

Ia = (
    I0_rms
    * sqrt(2)
    * cos(2 * pi * felec * time + 0 * rot_dir * 2 * pi / qs + Phi0)
)
Ib = (
    I0_rms
    * sqrt(2)
    * cos(2 * pi * felec * time + 1 * rot_dir * 2 * pi / qs + Phi0)
)
Ic = (
    I0_rms
    * sqrt(2)
    * cos(2 * pi * felec * time + 2 * rot_dir * 2 * pi / qs + Phi0)
)
simu_femm.input.Is = array([Ia, Ib, Ic]).transpose()

In this example stator currents are enforced as a function of time for each phase.

MagFEMM configuration

For the configuration of the Magnetic module, we use the object MagFEMM that computes 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 parameters and in particular the B(H) curve are setup directly in the machine.

[3]:
from pyleecan.Classes.MagFEMM import MagFEMM

simu_femm.mag = MagFEMM(
    type_BH_stator=0, # 0 to use the material 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 material B(H) curve,
                      # 1 to use linear B(H) curve according to mur_lin,
                      # 2 to enforce infinite permeability (mur_lin =100000)
    file_name = "", # Name of the file to save the FEMM model
)

# We only use the magnetic part
simu_femm.force = None
simu_femm.struct = None

Pyleecan coupling with FEMM enables to define the machine with symmetry and with sliding band to optimize the computation time. The angular periodicity of the machine will be computed and (in the particular case) only 1/8 of the machine will be drawn (4 symmetry + antiperiodicity):

[4]:
simu_femm.mag.is_periodicity_a=True

The same is done for time periodicity only half of one electrical period is calculated (i.e: 1/8 of mechanical period):

[5]:
simu_femm.mag.is_periodicity_t=True

Pyleecan enable to parallelize the call to FEMM by simply setting:

[6]:
simu_femm.mag.nb_worker = 4  # Number of FEMM instances to run at the same time (1 by default)

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

[7]:
simu_femm.mag.is_get_mesh = True # To get FEA mesh for latter post-procesing
simu_femm.mag.is_save_FEA = False # To save FEA results in a dat file

Run simulation

[8]:
out_femm = simu_femm.run()
2020-11-16 23:07:29,230-INFO-Pyleecan.Simulation: Starting Magnetic module

When running the simulation, an FEMM window runs in background. You can open it to see pyleecan drawing the machine and defining the surfaces. image0 The simulation will compute 32*p/8 different timesteps by updating the current and the sliding band boundary condition.

Once the simulation is finished, an Output object is return. The results are stored in the magnetic part of the output (i.e. out_femm.mag ) and different plots can be called. This out_femm.mag contains:
- Time: magnetic time axis
- Angle: magnetic position
- B: airgap flux density (contains radial and tangential components)
- Tem: electromagnetic torque
- Tem_av: average electromagnetic torque - Tem_rip_pp : Peak to Peak Torque ripple - Tem_rip_norm: Peak to Peak Torque ripple normalized according to average torque
- Phi_wind_stator: stator winding flux
- emf: electromotive force

Some of these properties are “Data objects” from the SciDataTool project. These object enables to handle unit conversion, interpolation, fft, periodicity…

Plot results

Output object embbed different plots to visualize results easily. A dedicated tutorial is available here.

For instance, the radial and tangential magnetic flux in the airgap at a specific timestep can be plotted with:

[9]:
# Radial magnetic flux
out_femm.plot_2D_Data("mag.B","angle","time[1]",component_list=["radial"])
out_femm.plot_2D_Data("mag.B","wavenumber=[0,76]","time[1]",component_list=["radial"])
[10]:
# Tangential magnetic flux
out_femm.plot_2D_Data("mag.B","angle","time[1]",component_list=["tangential"])
out_femm.plot_2D_Data("mag.B","wavenumber=[0,76]","time[1]",component_list=["tangential"])

The torque can be plotted with:

[11]:
out_femm.plot_2D_Data("mag.Tem", "time")

One can notice that the torque matrix includes the periodicity (only the meaningful part is stored)

[12]:
print(out_femm.mag.Tem.values.shape)
print(simu_femm.input.Nt_tot)
(16,)
128

If the mesh was saved in the output object (mySimu.mag.is_get_mesh = True), it can be plotted with:

[13]:
out_femm.mag.meshsolution.plot_contour(label="B", group_names="stator")
WARNING:root:VTK requires 3D points, but 2D points given. Appending 0 third component.

Finally, it is possible to extend pyleecan by implementing new plot by using the results from output. For instance, the following plot requires plotly to display the radial flux density in the airgap over time and angle.

[15]:
#%run -m pip install plotly # Uncomment this line to install plotly
import plotly.graph_objects as go
from plotly.offline import init_notebook_mode
init_notebook_mode()

result = out_femm.mag.B.components["radial"].get_along("angle{°}", "time")
x = result["angle"]
y = result["time"]
z = result["B_r"]
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',
                  autosize=True,
                  scene = dict(
                      xaxis_title='Angle [°]',
                      yaxis_title='Time [s]',
                      zaxis_title='Flux [T]'
                  ),
                  width=700,
                  margin=dict(r=20, b=100, l=10, t=100),
                 )

fig.show(config = {"displaylogo":False})