Version information

%matplotlib notebook
from PySide2.QtWidgets import *
from datetime import date
print("Running date:", date.today().strftime("%B %d, %Y"))
import pyleecan
print("Pyleecan version:" + pyleecan.__version__)
import SciDataTool
print("SciDataTool version:" + SciDataTool.__version__)

Running date: April 27, 2023
Pyleecan version:1.5.0
SciDataTool version:2.5.0

How to set the Operating Point

This tutorial explains how to use the object OPdq, InputCurrent and VarLoadCurrent to run a magnetic simulation on several operating points by setting Id/Iq or I0/Phi0.

The reference used to validate this tutorial is: Z. Yang, M. Krishnamurthy and I. P. Brown, “Electromagnetic and vibrational characteristic of IPM over full torque-speed range,” 2013 International Electric Machines & Drives Conference, Chicago, IL, 2013, pp. 295-302, doi: 10.1109/IEMDC.2013.6556267.

Machine and Simulation definition

This tutorial use the machine Toyota_Prius (2004) defined in the “How to define a machine” tutorial. The magnetic module is the same as the symmetrical one from the tutorial “How to define a simulation to call FEMM”.

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

Toyota_Prius = load(join(DATA_DIR, "Machine", "Toyota_Prius.json"))
fig, ax = Toyota_Prius.plot()

from pyleecan.Classes.Simu1 import Simu1
from pyleecan.Classes.MagFEMM import MagFEMM
# Initialization of the Simulation
simu_op = Simu1(name="tuto_Id_Iq", machine=Toyota_Prius)

# Definition of the magnetic simulation (FEMM with symmetry and sliding band)
simu_op.mag = MagFEMM(
    type_BH_stator=0,
    type_BH_rotor=0,
    is_periodicity_a=True,
    is_periodicity_t=True,
    nb_worker=4,
    Kgeo_fineness=1,
)
# Run only Magnetic module
simu_op.elec = None
simu_op.force = None
simu_op.struct = None

Defining an Operating point with Id/Iq

The InputCurrent object enable to create an “OutElec” object that corresponds to the output of the Electrical module and the input of the Magnetic module. In this example, the Operating Point is set with an OPdq object to define the starting point with a sinusoidal current defined with Id_ref and Iq_ref:

The tutorial “How to define a simulation to call FEMM” uses the same InputCurrent object to enforce any current by directly setting Is

from pyleecan.Classes.InputCurrent import InputCurrent
from pyleecan.Classes.OPdq import OPdq
from numpy import sqrt, exp,pi

# Definition of a sinusoidal current
simu_op.input = InputCurrent()
# I0, Phi0 to set
I0_rms = 250/sqrt(2) # Maximum current [Arms]
Phi0 = 140*pi/180  # Maximum Torque Per Amp
# Compute corresponding Id/Iq
Id_ref = (I0_rms*exp(1j*Phi0)).real # [Arms]
Iq_ref = (I0_rms*exp(1j*Phi0)).imag # [Arms]
# Setting the values
simu_op.input.OP = OPdq(Id_ref=Id_ref, Iq_ref=Iq_ref)

(Id_ref,Iq_ref)

(-135.4188051049254, 113.62986941801093)

The discretization of the current and for the magnetic computation can be set with time and angle (as in “How to define a simulation to call FEMM” tutorial) or by setting the following parameters:

simu_op.input.Nt_tot = 128 # Number of time step
simu_op.input.Na_tot = 2048 # Spatial discretization
simu_op.input.OP.N0 = 2000 # Rotor speed [rpm]

When Nt_tot is defined, the time vector is automatically set to:

linspace(0, 60 / N0 * Nrev, Nt_tot)

With Nrev the number of revolution of the rotor (1 by default)

When Na_tot is defined, the angle vector is automatically set to:

linspace(0, 2*pi, Na_tot)

The input is now fully defined, the simulation can now be run:

from pyleecan.Functions.Plot import dict_2D

out_op = simu_op.run()

[15:57:36] Starting running simulation tuto_Id_Iq (machine=Toyota_Prius)
[15:57:36] Starting Magnetic module
[15:57:38] Computing Airgap Flux in FEMM
[15:57:45] End of simulation tuto_Id_Iq

# Plot the flux
out_op.mag.B.plot_2D_Data("angle", **dict_2D)
# Plot the torque
out_op.mag.Tem.plot_2D_Data("time", **dict_2D)
# Plot the current
out_op.elec.get_Is().plot_2D_Data("time", "phase[]", **dict_2D)

The Operating Point can also be defined directly with I0 / Phi0 with:

from numpy import pi

simu_op.input.OP.set_I0_Phi0(I0=I0_rms, Phi0=Phi0)
simu_op.input.OP.get_Id_Iq()

{'Id': -135.4188051049254, 'Iq': 113.62986941801093}

Iterating on several Operating Point

Each pyleecan simulation is assumed to be quasi-static and run on a single operating point (fixed speed). To run a simulation on several operating points two steps are needed: First define a simulation that run correctly on a single operating point (like the one defined above), then define a VarLoadCurrent object.

The VarLoadCurrent object is defined with a matrix with each line corresponding to an operating point (so each line will run a single speed simulation). The matrix columns order/containt are set by the user during the call to the method set_OP_matrix for more flexibility. Usually the first 3 columns are (N0, Id, Iq) or (N0 ,I0 ,Phi0)

A torque column can be added by adding “Tem” as a column name. It enables to define the reference torque for the Operating Point. The reference is stored in output.elec.Tem_av_ref, the real computed torque is stored in output.mag.Tem_av.

The following VarLoadCurrent object will run the previous simulation N_simu times by changing the value of Phi0.

Reference torque and current angle vector are:

from numpy import linspace, array, pi

Tem_av_ref = array([79, 125, 160, 192, 237, 281, 319, 343, 353, 332, 266, 164, 22]) # Yang et al, 2013
Phi0_ref = linspace(60 * pi / 180, 180 * pi / 180, Tem_av_ref.size)
N_simu = Tem_av_ref.size

from pyleecan.Classes.VarLoadCurrent import VarLoadCurrent
from pyleecan.Classes.DataKeeper import DataKeeper
from numpy import zeros, ones, linspace, array, sqrt, arange

# Create new Variable speed simulation
simu_vop = simu_op.copy()
simu_vop.mag.import_file = None
varload = VarLoadCurrent()
simu_vop.var_simu = varload

# Creating the Operating point matrix
OP_matrix = zeros((N_simu,4))

# Set N0 = 2000 [rpm] for all simulation
OP_matrix[:,0] = 2000 * ones((N_simu))

# Set I0 = 250 / sqrt(2) [A] (RMS) for all simulation
OP_matrix[:,1] = I0_rms * ones((N_simu))

# Set Phi0 from 60° to 180°
OP_matrix[:,2] = Phi0_ref

# Set reference torque from Yang et al, 2013
OP_matrix[:,3] = Tem_av_ref
print("Operating Point matrix:")
print(OP_matrix)
print("Original Operating point:")
print(simu_vop.input.OP.get_I0_Phi0())
# Column name are given in the order of OP_matrix
varload.set_OP_array(OP_matrix, "N0", "I0", "Phi0", "Tem")
# The operating point of the simulation is automatically set as
# the first one from the OP_matrix (to speed-up computation due to how pyleecan works)
print("Updated Operating point:")
print(simu_vop.input.OP.get_I0_Phi0())

# All the simulation use the same machine
# No need to draw the machine for all OP
varload.is_reuse_femm_file=True

# Add I0/Phi0 to the data available at the end of the simulation
I0_dk = DataKeeper(
    name="Stator current rms amplitude",
    symbol="I0",
    unit="Arms",
    keeper="lambda output: output.elec.OP.get_I0_Phi0()['I0']",
)
Phi0_dk = DataKeeper(
    name="Stator current phase",
    symbol="Phi0",
    unit="rad",
    keeper="lambda output: output.elec.OP.get_I0_Phi0()['Phi0']",
)
varload.datakeeper_list= [I0_dk, Phi0_dk]

Operating Point matrix:
[[2.00000000e+03 1.76776695e+02 1.04719755e+00 7.90000000e+01]
 [2.00000000e+03 1.76776695e+02 1.22173048e+00 1.25000000e+02]
 [2.00000000e+03 1.76776695e+02 1.39626340e+00 1.60000000e+02]
 [2.00000000e+03 1.76776695e+02 1.57079633e+00 1.92000000e+02]
 [2.00000000e+03 1.76776695e+02 1.74532925e+00 2.37000000e+02]
 [2.00000000e+03 1.76776695e+02 1.91986218e+00 2.81000000e+02]
 [2.00000000e+03 1.76776695e+02 2.09439510e+00 3.19000000e+02]
 [2.00000000e+03 1.76776695e+02 2.26892803e+00 3.43000000e+02]
 [2.00000000e+03 1.76776695e+02 2.44346095e+00 3.53000000e+02]
 [2.00000000e+03 1.76776695e+02 2.61799388e+00 3.32000000e+02]
 [2.00000000e+03 1.76776695e+02 2.79252680e+00 2.66000000e+02]
 [2.00000000e+03 1.76776695e+02 2.96705973e+00 1.64000000e+02]
 [2.00000000e+03 1.76776695e+02 3.14159265e+00 2.20000000e+01]]
Original Operating point:
{'I0': 176.77669529663686, 'Phi0': 2.443460952792061}
Updated Operating point:
{'I0': 176.77669529663686, 'Phi0': 1.0471975511965976}

The original simulation will be duplicated N_simu times with the value of InputCurrent updated according to the matrix.

# Run simulation is the same as for the fixed speed case
Xout = simu_vop.run()

[15:58:22] Starting running simulation tuto_Id_Iq (machine=Toyota_Prius)
[15:58:22] Computing reference simulation for Variable Load
[15:58:22] Starting Magnetic module
[15:58:23] Computing Airgap Flux in FEMM
[15:58:30] Variable Load Reference Results: I0=176.8 [Arms], Phi0=1.047 [rad], N0=2000 [rpm], Id=88.39 [Arms], Iq=153.1 [Arms], Tem_av_ref=79 [N.m], Tem_av=90.37 [N.m], Tem_rip_pp=30.72 [N.m], Tem_rip_norm=0.34 [-], Pem_av=1.893e+04 [W]
[===                                               ]   7%
[15:58:30] Variable Load: Using same FEMM file for all simulations (C:\Users\LAP21\Documents\pyleecan-public2\pyleecan\Results\tuto_Id_Iq\Femm\Toyota_Prius_model.fem)
[15:58:30] Running simulation 1/13 with N0=2000 [rpm], Id=88.39 [Arms], Iq=153.1 [Arms]
[15:58:30] Simulation matches reference one: Skipping computation
[15:58:30] Variable Load Results: I0=176.8 [Arms], Phi0=1.047 [rad], N0=2000 [rpm], Id=88.39 [Arms], Iq=153.1 [Arms], Tem_av_ref=79 [N.m], Tem_av=90.37 [N.m], Tem_rip_pp=30.72 [N.m], Tem_rip_norm=0.34 [-], Pem_av=1.893e+04 [W]
[=======                                           ]  14%
[15:58:30] Running simulation 2/13 with N0=2000 [rpm], Id=60.46 [Arms], Iq=166.1 [Arms]
[15:58:30] Starting Magnetic module
[15:58:31] Computing Airgap Flux in FEMM
[15:58:36] Variable Load Results: I0=176.8 [Arms], Phi0=1.222 [rad], N0=2000 [rpm], Id=60.46 [Arms], Iq=166.1 [Arms], Tem_av_ref=125 [N.m], Tem_av=134.1 [N.m], Tem_rip_pp=47.59 [N.m], Tem_rip_norm=0.3549 [-], Pem_av=2.808e+04 [W]
[==========                                        ]  21%
[15:58:36] Running simulation 3/13 with N0=2000 [rpm], Id=30.7 [Arms], Iq=174.1 [Arms]
[15:58:36] Starting Magnetic module
[15:58:36] Computing Airgap Flux in FEMM
[15:58:42] Variable Load Results: I0=176.8 [Arms], Phi0=1.396 [rad], N0=2000 [rpm], Id=30.7 [Arms], Iq=174.1 [Arms], Tem_av_ref=160 [N.m], Tem_av=179 [N.m], Tem_rip_pp=57.22 [N.m], Tem_rip_norm=0.3197 [-], Pem_av=3.748e+04 [W]
[==============                                    ]  28%
[15:58:42] Running simulation 4/13 with N0=2000 [rpm], Id=1.082e-14 [Arms], Iq=176.8 [Arms]
[15:58:42] Starting Magnetic module
[15:58:42] Computing Airgap Flux in FEMM
[15:58:48] Variable Load Results: I0=176.8 [Arms], Phi0=1.571 [rad], N0=2000 [rpm], Id=1.082e-14 [Arms], Iq=176.8 [Arms], Tem_av_ref=192 [N.m], Tem_av=223 [N.m], Tem_rip_pp=55.5 [N.m], Tem_rip_norm=0.2489 [-], Pem_av=4.669e+04 [W]
[=================                                 ]  35%
[15:58:48] Running simulation 5/13 with N0=2000 [rpm], Id=-30.7 [Arms], Iq=174.1 [Arms]
[15:58:48] Starting Magnetic module
[15:58:48] Computing Airgap Flux in FEMM
[15:58:54] Variable Load Results: I0=176.8 [Arms], Phi0=1.745 [rad], N0=2000 [rpm], Id=-30.7 [Arms], Iq=174.1 [Arms], Tem_av_ref=237 [N.m], Tem_av=264.2 [N.m], Tem_rip_pp=58.8 [N.m], Tem_rip_norm=0.2226 [-], Pem_av=5.533e+04 [W]
[=====================                             ]  42%
[15:58:54] Running simulation 6/13 with N0=2000 [rpm], Id=-60.46 [Arms], Iq=166.1 [Arms]
[15:58:54] Starting Magnetic module
[15:58:55] Computing Airgap Flux in FEMM
[15:59:00] Variable Load Results: I0=176.8 [Arms], Phi0=1.92 [rad], N0=2000 [rpm], Id=-60.46 [Arms], Iq=166.1 [Arms], Tem_av_ref=281 [N.m], Tem_av=301 [N.m], Tem_rip_pp=73.05 [N.m], Tem_rip_norm=0.2427 [-], Pem_av=6.303e+04 [W]
[=========================                         ]  50%
[15:59:00] Running simulation 7/13 with N0=2000 [rpm], Id=-88.39 [Arms], Iq=153.1 [Arms]
[15:59:00] Starting Magnetic module
[15:59:00] Computing Airgap Flux in FEMM
[15:59:06] Variable Load Results: I0=176.8 [Arms], Phi0=2.094 [rad], N0=2000 [rpm], Id=-88.39 [Arms], Iq=153.1 [Arms], Tem_av_ref=319 [N.m], Tem_av=331.3 [N.m], Tem_rip_pp=82.32 [N.m], Tem_rip_norm=0.2484 [-], Pem_av=6.939e+04 [W]
[============================                      ]  57%
[15:59:06] Running simulation 8/13 with N0=2000 [rpm], Id=-113.6 [Arms], Iq=135.4 [Arms]
[15:59:06] Starting Magnetic module
[15:59:07] Computing Airgap Flux in FEMM
[15:59:12] Variable Load Results: I0=176.8 [Arms], Phi0=2.269 [rad], N0=2000 [rpm], Id=-113.6 [Arms], Iq=135.4 [Arms], Tem_av_ref=343 [N.m], Tem_av=352.7 [N.m], Tem_rip_pp=85.22 [N.m], Tem_rip_norm=0.2416 [-], Pem_av=7.388e+04 [W]
[================================                  ]  64%
[15:59:12] Running simulation 9/13 with N0=2000 [rpm], Id=-135.4 [Arms], Iq=113.6 [Arms]
[15:59:12] Starting Magnetic module
[15:59:12] Computing Airgap Flux in FEMM
[15:59:19] Variable Load Results: I0=176.8 [Arms], Phi0=2.443 [rad], N0=2000 [rpm], Id=-135.4 [Arms], Iq=113.6 [Arms], Tem_av_ref=353 [N.m], Tem_av=360.6 [N.m], Tem_rip_pp=85.15 [N.m], Tem_rip_norm=0.2361 [-], Pem_av=7.552e+04 [W]
[===================================               ]  71%
[15:59:19] Running simulation 10/13 with N0=2000 [rpm], Id=-153.1 [Arms], Iq=88.39 [Arms]
[15:59:19] Starting Magnetic module
[15:59:19] Computing Airgap Flux in FEMM
[15:59:25] Variable Load Results: I0=176.8 [Arms], Phi0=2.618 [rad], N0=2000 [rpm], Id=-153.1 [Arms], Iq=88.39 [Arms], Tem_av_ref=332 [N.m], Tem_av=343.5 [N.m], Tem_rip_pp=67.67 [N.m], Tem_rip_norm=0.197 [-], Pem_av=7.193e+04 [W]
[=======================================           ]  78%
[15:59:25] Running simulation 11/13 with N0=2000 [rpm], Id=-166.1 [Arms], Iq=60.46 [Arms]
[15:59:25] Starting Magnetic module
[15:59:26] Computing Airgap Flux in FEMM
[15:59:32] Variable Load Results: I0=176.8 [Arms], Phi0=2.793 [rad], N0=2000 [rpm], Id=-166.1 [Arms], Iq=60.46 [Arms], Tem_av_ref=266 [N.m], Tem_av=276.6 [N.m], Tem_rip_pp=38.65 [N.m], Tem_rip_norm=0.1397 [-], Pem_av=5.794e+04 [W]
[==========================================        ]  85%
[15:59:32] Running simulation 12/13 with N0=2000 [rpm], Id=-174.1 [Arms], Iq=30.7 [Arms]
[15:59:32] Starting Magnetic module
[15:59:33] Computing Airgap Flux in FEMM
[15:59:39] Variable Load Results: I0=176.8 [Arms], Phi0=2.967 [rad], N0=2000 [rpm], Id=-174.1 [Arms], Iq=30.7 [Arms], Tem_av_ref=164 [N.m], Tem_av=157 [N.m], Tem_rip_pp=47.61 [N.m], Tem_rip_norm=0.3032 [-], Pem_av=3.289e+04 [W]
[==============================================    ]  92%
[15:59:39] Running simulation 13/13 with N0=2000 [rpm], Id=-176.8 [Arms], Iq=2.165e-14 [Arms]
[15:59:39] Starting Magnetic module
[15:59:39] Computing Airgap Flux in FEMM
[15:59:45] Variable Load Results: I0=176.8 [Arms], Phi0=3.142 [rad], N0=2000 [rpm], Id=-176.8 [Arms], Iq=2.165e-14 [Arms], Tem_av_ref=22 [N.m], Tem_av=-0.4102 [N.m], Tem_rip_pp=39.42 [N.m], Tem_rip_norm=-96.1 [-], Pem_av=-85.92 [W]
[==================================================] 100%
[15:59:45] End of simulation tuto_Id_Iq

Pyleecan will automatically extract some values from each simulation. These values are all gathered in the xoutput_dict:

print("Values available in XOutput:")
print(Xout.xoutput_dict.keys())
print(varload.OP_matrix.col_names)

print("\nI0 for each simulation:")
print(Xout["I0"].result)
print("\nPhi0 for each simulation:")
print(Xout["Phi0"].result)

Values available in XOutput:
dict_keys(['I0', 'Phi0', 'N0', 'Id', 'Iq', 'Tem_av_ref', 'Tem_av', 'Tem_rip_pp', 'Tem_rip_norm', 'Pem_av'])
['N0', 'Id', 'Iq', 'Tem']

I0 for each simulation:
[176.77669529663686, 176.77669529663686, 176.77669529663683, 176.77669529663686, 176.77669529663683, 176.77669529663686, 176.77669529663686, 176.77669529663686, 176.77669529663686, 176.77669529663686, 176.77669529663683, 176.77669529663683, 176.77669529663686]

Phi0 for each simulation:
[1.0471975511965976, 1.2217304763960306, 1.3962634015954636, 1.5707963267948966, 1.7453292519943295, 1.9198621771937625, 2.0943951023931957, 2.2689280275926285, 2.443460952792061, 2.6179938779914944, 2.7925268031909276, 2.9670597283903604, 3.141592653589793]

Any parameter in the XOutput can be plot as a function of any other

fig_1, ax_1 = Xout.plot_multi("Phi0", "Tem_av")
fig_2, ax_2 = Xout.plot_multi("Id", "Iq")

Finally, the computed average torque can be compared to the one in the publication from Yang et al (data has been extracted from their graph using Engauge Digitizer. Note that the generic plot function plot_2D has been used here.

from SciDataTool.Functions.Plot.plot_2D import plot_2D
from pyleecan.definitions import config_dict
from numpy import array

plot_2D(
    array([x*180/pi for x in Xout.xoutput_dict["Phi0"].result]),
    [Xout.xoutput_dict["Tem_av"].result, Xout.xoutput_dict["Tem_av_ref"].result],
    legend_list=["Pyleecan", "Yang et al, 2013"],
    xlabel="Current angle [°]",
    ylabel="Electrical torque [N.m]",
    title="Electrical torque vs current angle",
    **dict_2D
)