diff --git a/docs/source/EM_analyzers/Images/FluxLinkageExampleCircuit.png b/docs/source/EM_analyzers/Images/FluxLinkageExampleCircuit.png new file mode 100644 index 00000000..571e38fe Binary files /dev/null and b/docs/source/EM_analyzers/Images/FluxLinkageExampleCircuit.png differ diff --git a/docs/source/EM_analyzers/Images/a_b_c_inductances.svg b/docs/source/EM_analyzers/Images/a_b_c_inductances.svg new file mode 100644 index 00000000..3770222c --- /dev/null +++ b/docs/source/EM_analyzers/Images/a_b_c_inductances.svg @@ -0,0 +1,1528 @@ + + + + + + + + 2024-06-01T17:13:38.224360 + image/svg+xml + + + Matplotlib v3.5.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/EM_analyzers/Images/alpha_beta_inductances.svg b/docs/source/EM_analyzers/Images/alpha_beta_inductances.svg new file mode 100644 index 00000000..172ac91f --- /dev/null +++ b/docs/source/EM_analyzers/Images/alpha_beta_inductances.svg @@ -0,0 +1,1434 @@ + + + + + + + + 2024-06-01T17:11:38.514690 + image/svg+xml + + + Matplotlib v3.5.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/EM_analyzers/Images/d_q_inductances.svg b/docs/source/EM_analyzers/Images/d_q_inductances.svg new file mode 100644 index 00000000..2a496b0d --- /dev/null +++ b/docs/source/EM_analyzers/Images/d_q_inductances.svg @@ -0,0 +1,1452 @@ + + + + + + + + 2024-06-01T17:11:49.344411 + image/svg+xml + + + Matplotlib v3.5.1, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/source/EM_analyzers/flux_linkage_analyzer.rst b/docs/source/EM_analyzers/flux_linkage_analyzer.rst new file mode 100644 index 00000000..02b2ea8d --- /dev/null +++ b/docs/source/EM_analyzers/flux_linkage_analyzer.rst @@ -0,0 +1,147 @@ +Flux Linkage Analyzer +######################################################################## + +This analyzer enables the flux linkage evaluation of a mutli-phase electric machine after running 2D FEA simulations using JMAG. + +Model Background +**************** + +The flux linkage of a coil is defined as the amount of flux linking together for a multi-coil arrangment with electric current flowing +through them. The flux linkage of a coil within an electric machine comes from all coils present in the machine and has a profound +impact on the machine characteristics. Calculating coil flux linkages over time can lead to inductance calculations for an electric +machine, which are also important for characterizing that machine. The flux linkage is an important parameter for inductance calculations +as can be seen in the following equation: + +.. math:: + + L = \lambda I \\ + +where :math:`\lambda` is the flux linkage, :math:`L` is the inductance, and :math:`I` is the coil current. + +The code is structured such that the ``flux_linkage_analyzer`` contains the code for setting up and running the JMAG simulations based on +1) the machine inputs and conditions of the user and 2) the conditions required of the machine to be able to calculate the +necessary parameters. In the case of each machine, DC excitement of each phase coil occurs with all other coils open. This repeats for +each coil until all of the coils have been excited with DC current and all of the flux linkages have been captured. + +This analyzer calculates the self and mutual flux linkages of each coil using JMAG's transient solver. It models an electric machine +under synchronous operation. The following information document will provide a description of the analyzer inputs and outputs. + +Input from User +********************************* + +Users are utilizing a single Problem class to interface with this analyzer. This class requires the user to provide an instance of the +JMAG FEA application being used, a the model of the current project which needs to be in a file that is already open, a dedicated filepath +for the results, the name of the phases, and the rated current of the machine. The specific requirements are summarized below: + +.. csv-table:: `flux_linkage_analyzer Inputs` + :file: input_flux_linkage_analyzer.csv + :widths: 70, 70, 30 + :header-rows: 1 + +With regards to the model properties, there are several aspects of the JMAG FEA model that must be specific such that this analyzer works +properly. Those properties are as follows + +1. Must be transient model (can be 2D or 3D) +2. The model must be named "Machine_FluxLinkage_Project" +3. The study must be named "Machine_FluxLinkage_Study" +4. The motion condition must be specified and MUST be rotating 1 full mechanincal revolution +5. "FEM Coils" must be applied to the winding and linked to the circuit +6. The study properties must be fully defined, except for the csv output, which is defined by the analyzer +7. The circuit must appear similar to the following image, where the current sources are titled "cs_PhaseName" +8. Under case control, the "cs_PhaseName" amplitudes MUST be the first cases AND listed in order of excitation + +.. figure:: ./Images/FluxLinkageExampleCircuit.png + :alt: Stator Diagram + :align: center + :width: 250 + +Example Code +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The following example demonstrates how to initialize instances of ``Flux_Linkage_Problem`` and ``Flux_Linkage_Analyzer``. An instance of +JMAG FEA is used in this example and is stored in the ``fluxlinkage_inductance_eval`` folder of the ``mach_eval_examples`` in ``eMach``. +An example file containing the example code is also found in that folder. The following code runs the flux linkage analyzer using the +aforementioned example file: + +.. code-block:: python + + import os + import sys + import numpy as np + import matplotlib.pyplot as plt + from time import time as clock_time + + os.chdir(os.path.dirname(__file__)) + sys.path.append("../../../") + + from mach_eval.analyzers.electromagnetic.flux_linkage_analyzer import FluxLinkageJMAG_Problem, FluxLinkageJMAG_Analyzer + + from mach_cad.tools import jmag as JMAG + + filepath = "eMach_location/eMach/examples/mach_eval_examples/fluxlinkage_inductance_eval" + phasenames = ['U', 'V', 'W'] + rated_current = 20 + + #################################################### + # 01 Setting project name and output folder + #################################################### + + toolJmag = JMAG.JmagDesigner() + toolJmag.visible = True + toolJmag.open(filepath + "/Example_FluxLinkage_Machine.jproj") + toolJmag.save() + +This example code does the following: +1. Initializes all of the required libraries and classes +2. Defines the necessary inputs of the ``problem`` and ``analyzer`` classes +3. Opens a fully-defined instance of JMAG +4. Defines the output file location based +5. Names project accordingly + +Output to User +********************************** + +The ``flux_linkage_analyzer`` returns a directory holding the results obtained from the transient analysis of the machine. The elements +of this dictionary and their descriptions are provided below: + +.. csv-table:: `flux_linkage_analyzer Output` + :file: output_flux_linkage_analyzer.csv + :widths: 70, 70 + :header-rows: 1 + +The following code should be used to run the example analysis: + +.. code-block:: python + + ############################ Create Evaluator ##################### + tic = clock_time() + flux_linkage_prob = FluxLinkageJMAG_Problem(toolJmag, phasenames, rated_current) + flux_linkage_analyzer = FluxLinkageJMAG_Analyzer() + flux_linkages, currents = flux_linkage_analyzer.analyze(flux_linkage_prob) + toc = clock_time() + print("Time spent on the flux linkage evaluation is %g min." % ((toc- tic)/60)) + + linkages = flux_linkages["linkages"] + phase_currents = currents + rotor_angle = flux_linkages["rotor_angle"][0] + name_of_phases = flux_linkages["name_of_phases"] + + print("\n************************ FLUX LINKAGE RESULTS ************************") + print("Linkages = ", linkages) + print("phase_currents = ", phase_currents, " A") + print("rotor_angle = ", rotor_angle, " deg") + print("name_of_phases = ", name_of_phases) + print("*************************************************************************\n") + +This example, contained in the aforementioned ``fluxlinkage_inductance_eval`` folder, should produce the following results: + +.. csv-table:: `flux_linkage_analyzer Results` + :file: results_flux_linkage_analyzer.csv + :widths: 70, 70, 30 + :header-rows: 1 + +One should expect the ``/run_data/`` working folder location to differ depending on where the workspace is. Within ``/run_data/`` there should be a +total of 4 csv files that contains the flux linkage calculations for a 3 phase machine, there should be 7 csv files for a 6 phase machine, etc. Each +csv files should contain a total number of columns that equals the phase count of the machine. All of the code shown exists in the +``fluxlinkage_inductance_evaluator.py`` file in the ``eMach/examples/mach_eval_examples/fluxlinkage_inductance_eval`` folder. This analyzer serves +as a first step in conjunction with the `Inductance Analyzer `_. \ No newline at end of file diff --git a/docs/source/EM_analyzers/index.rst b/docs/source/EM_analyzers/index.rst index 3d319c8c..3c49b2db 100644 --- a/docs/source/EM_analyzers/index.rst +++ b/docs/source/EM_analyzers/index.rst @@ -15,4 +15,5 @@ This section provides documentation for the electromagnetic analyzers supported BSPM JMAG 2D FEA Winding Factors SynR JMAG 2D FEA - Inductance/Saliency \ No newline at end of file + Flux Linkage + Inductance \ No newline at end of file diff --git a/docs/source/EM_analyzers/inductance_analyzer.rst b/docs/source/EM_analyzers/inductance_analyzer.rst index d3bb7582..6e2192b6 100644 --- a/docs/source/EM_analyzers/inductance_analyzer.rst +++ b/docs/source/EM_analyzers/inductance_analyzer.rst @@ -1,7 +1,7 @@ -Inductance/Saliency Analyzer +Inductance Analyzer ######################################################################## -This analyzer enables the inductance evaluation of an electric machine after running 2D FEA simulations using JMAG. +This analyzer enables the inductance evaluation of a multi-phase electric machine using raw data from 2D FEA simulations using JMAG. Model Background **************** @@ -11,7 +11,7 @@ within an electric machine can come from multiple sources, including its own ele flowing through other phases of the machine. Understanding the inductance characteristics of an electric machine leads to finding any saliency that an electric machine rotor may have. In some electric machines, such as reluctance or induction machines, saliency exists and aids in producing electromagnetic torque between the rotor and stator. This can be seen in the torque equation -for a 3-phase machine: +for an electric machine: .. math:: @@ -25,243 +25,154 @@ q-axis currents :math:`I_\text{d,q}`, and the d- and q-axis inductances :math:`L :math:`\frac{3p}{2} (L_\text{d} - L_\text{q}) I_\text{d} I_\text{q}` represents the reluctance torque, which is generated by different d- and q-axis inductances. -The code is structured such that the ``inductance_analyzer`` contains the code for setting up and running the JMAG simulations based on -1) the machine inputs and conditions of the user and 2) the conditions required of the machine to be able to calculate the -necessary parameters. In the case of this machine, DC excitement of the U-phase is required with both the V- and W-phases being open. -The post-analyzer script post-processes the .csv file generated to take the self-inductance of the U-phase coil and mutual-inductance of -the V-phase coil to calcualte :math:`L_\text{d}` and :math:`L_\text{q}`. This is done using the following equations: +The code is structured such that the ``inductance_analyzer`` contains the code for taking flux linkage data in the form of .csv files +and processing it into alpha-beta and d-q inductance arrays. The example below utilitzes the existing +`flux linkage analyzer `_. The analyzer uses the +following equations to find :math:`L_{\alpha \beta \gamma}` and :math:`L_\text{dq0}`. These are the standard equations for the Clarke +and Park transforms, in this case applied to each phase of the machine calculated at each individual rotor angle: .. math:: - L_\text{d} &= L_\text{ls} + \frac{3}{2}(L_\text{0} - L_\text{g}) \\ - L_\text{q} &= L_\text{ls} + \frac{3}{2}(L_\text{0} + L_\text{g}) \\ + [L_{\alpha \beta \gamma}] &= [T_\text{C}] [L_\text{abc}] \\ + [L_\text{dq0}] &= [T_\text{P}] [L_{\alpha \beta \gamma}] \\ -where :math:`L_\text{ls}` is the average value of the self leakage inductance, :math:`L_\text{0}` is the inductance component caused -by the air-gap magnetic field, and :math:`L_\text{g}` is the the amplitude of self/mutual inductance variation due to saliency. More -information and images depicting the relationships between these variables can be found using the reference at the conclusion of this -paragraph. Together, both :math:`L_\text{d}` and :math:`L_\text{q}` can be used to find the saliency ratio, which is defined as: - -.. math:: - - \xi &= \frac{L_\text{d}}{L_\text{q}} \\ - -* X. Qiu, and W. Wang, and J. Yang, and J. Jiang, and J. Yang, "Phase-Inductance-Based Position Estimation Method for Interior Permanent Magnet Synchronous Motors," Energies 2017, 10, 2002. - -This analyzer calculates :math:`L_\text{d}`, :math:`L_\text{q}`, and :math:`\xi` using JMAG's transient solver. It models a synchronous -reluctance machine under synchronous operation. The following document will provide a description of the analyzer inputs and outputs: +where :math:`[L_{\alpha \beta \gamma}]` is the alpha-beta inductance matrix, :math:`[T_\text{c}]` is the Clarke transformation matrix, +:math:`[L_\text{abc}]` is the synchronous frame inductance matrix, :math:`[L_\text{dq0}]` is the d-q inductance matrix, and :math:`[T_\text{c}]` +is the Park transformation matrix. This analyzer is designed to handle any number of phases, given the proper :math:`[T_\text{c}]` is +provided as an input. From these equations, post-processing of the data can find, for example, cross coupling of inductances and saliency +ratios. Input from User ********************************* -This analyzer is used in the same way as the ``SynR_JMAG_2D_FEA_Analyzer``. The inputs and initialization are the exact same and are shown -in the tables below: +This analyzer requires several input parameters, all of which are displayed in the table below: -.. csv-table:: `MachineDesign Input` - :file: input_SynR_jmag2d_analyzer.csv - :widths: 70, 70 - :header-rows: 1 - -.. csv-table:: `SynR_induction_analyzer Initialization` - :file: init_SynR_jmag2d_analyzer.csv - :widths: 70, 70 +.. csv-table:: `inductance_analyzer Input` + :file: input_inductance_analyzer.csv + :widths: 70, 70, 30 :header-rows: 1 -The example code initializing the machine design for the optimized SynR design provided in the example found in the ``SynR Design`` section of -``MACHINE DESIGNS`` is also identical, along with the ``example_SynR_machine`` file in the ``SynR_eval`` folder within the ``mach_eval_examples`` -folder in the ``examples`` folder of ``eMach``. This example code is used exactly the same, where the step within the ``SynR_evaluator`` file -should be changed to the ``inductance_step``, which is also contained in the aforementinoed folder. The objective of this file is to call the -example machine (in this case the ``example_SynR_machine.py`` that was just created in the ``SynR_eval`` folder) and create a machine design -object. - -.. code-block:: python - - import os - import sys - from time import time as clock_time +It should be noted that all of the inputs remain requirements regardless of the machine type. However, the Clarke transformation matrix +will be different depending on the number of phases. An adaptation to this analyzer will be required for any machine type with more than one +current type (i.e. BSPM machine with torque and suspension currents). - os.chdir(os.path.dirname(__file__)) +Example Code +~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - from eMach.mach_eval import (MachineEvaluator, MachineDesign) - from eMach.examples.mach_eval_examples.SynR_eval.inductance_step import inductance_step - from eMach.examples.mach_eval_examples.SynR_eval.example_SynR_machine import Example_SynR_Machine, Machine_Op_Pt - - ############################ Create Evaluator ######################## - SynR_evaluator = MachineEvaluator( - [ - inductance_step - ] - ) - - design_variant = MachineDesign(Example_SynR_Machine, Machine_Op_Pt) - - results = SynR_evaluator.evaluate(design_variant) - -Example code defining the inductance step is provided below. This code defines the analyzer problem class (input to the analyzer), -initializes the analyzer class with an explanation of the required configurations, and calls the post-analyzer class. +The following example demonstrates how to initialize instances of ``Inductance_Problem`` and ``Inductance_Analyzer``. An example file containing +the example code is stored in the ``fluxlinkage_inductance_eval`` folder of the ``mach_eval_examples`` in ``eMach``. The following code runs the +inductance analyzer aforementioned example file, which requires the use of the `Flux Linkage Analyzer `_: .. code-block:: python - import os - import sys - import copy - - from mach_eval import AnalysisStep, ProblemDefinition - from mach_eval.analyzers.electromagnetic.SynR import SynR_inductance_analyzer as SynR_inductance - from mach_eval.analyzers.electromagnetic.SynR.SynR_inductance_config import SynR_Inductance_Config - from examples.mach_eval_examples.SynR_eval.SynR_inductance_post_analyzer import SynR_Inductance_PostAnalyzer - - ############################ Define Inductance Step ########################### - class SynR_Inductance_ProblemDefinition(ProblemDefinition): - """Converts a State into a problem""" - - def __init__(self): - pass - - def get_problem(state): - - problem = SynR_inductance.SynR_Inductance_Problem( - state.design.machine, state.design.settings) - return problem - - # initialize inductance analyzer class with FEA configuration - configuration = SynR_Inductance_Config( - no_of_rev = 1, - no_of_steps = 72, - - mesh_size=3, # mm - mesh_size_rotor=1.5, # mm - airgap_mesh_radial_div=4, - airgap_mesh_circum_div=720, - mesh_air_region_scale=1.05, - - only_table_results=False, - csv_results=("FEMCoilFlux"), - del_results_after_calc=False, - run_folder=os.path.dirname(__file__) + "/run_data/", - jmag_csv_folder=os.path.dirname(__file__) + "/run_data/jmag_csv/", + from mach_eval.analyzers.electromagnetic.inductance_analyzer import Inductance_Problem, Inductance_Analyzer - max_nonlinear_iterations=50, - multiple_cpus=True, - num_cpus=4, - jmag_scheduler=False, - jmag_visible=True, - scale_axial_length = True, - jmag_version=None, - ) + linkages = flux_linkages["linkages"] + phase_currents = currents + rotor_angle = flux_linkages["rotor_angle"][0] + name_of_phases = flux_linkages["name_of_phases"] - SynR_inductance_analysis = SynR_inductance.SynR_Inductance_Analyzer(configuration) + clarke_transformation_matrix = 2/3*np.array([[1, -1/2, -1/2], [0, np.sqrt(3)/2, -np.sqrt(3)/2], [1/2, 1/2, 1/2]]) - inductance_step = AnalysisStep(SynR_Inductance_ProblemDefinition, SynR_inductance_analysis, SynR_Inductance_PostAnalyzer) +It should be noted that this code should be contained as an analysis step in the main folder of the eMach repository. It must be contained +within the same folder as the code below in order for the code below to run. Output to User ********************************** -The ``SynR_inductance_analyzer`` returns a dictionary holding the results obtained from the transient analysis of the machine. The elements +The ``inductance_analyzer`` returns a directory holding the results obtained from the transient analysis of the machine. The elements of this dictionary and their descriptions are provided below: -.. csv-table:: `SynR_inductance_analyzer Output` - :file: output_SynR_inductance_analyzer.csv +.. csv-table:: `inductance_analyzer Output` + :file: output_inductance_analyzer.csv :widths: 70, 70 :header-rows: 1 -As mentioned, the post analyzer is necessary to extract and compute the analyzer's computations and to interpret the results. The post analyzer -contains the following code and lies also in the ``eMach\examples\mach_eval_examples\SynR_eval`` folder. The code contained in the post analyzer, -in this case to find inductance quantities the saliency ratio, can be seen here: +The following code should be used to run the example analysis: .. code-block:: python - import copy - import numpy as np - import matplotlib.pyplot as plt - import scipy.optimize - - class SynR_Inductance_PostAnalyzer: - - def get_next_state(results, in_state): - state_out = copy.deepcopy(in_state) - - ############################ Extract required info ########################### - inductances = results["coil_inductances"] - I_hat = results["current_peak"] - - ############################ post processing ########################### - data = inductances.to_numpy() # change csv format to readable array - - t = data[:,0] # define x axis data as time - Uu = data[:,1] # define y axis data as self inductance - Uv = data[:,2] # define y axis data as mutual inductance - - # curve fit inductance values and calculate curve - def fit_sin(t, y): - fft_func = np.fft.fftfreq(len(t), (t[1]-t[0])) # define fft function with assumed uniform spacing - fft_y = abs(np.fft.fft(y)) # carry out fft function for inductance values - guess_freq = abs(fft_func[np.argmax(fft_y[1:])+1]) # excluding the zero frequency "peak", which can cause problematic fits - guess_amp = np.std(y) # guess amplitude based on one standard deviation - guess_offset = np.mean(y) # guess y offset based on average of magnitude - guess = np.array([guess_amp, 2.*np.pi*guess_freq, 0, guess_offset]) # arrage in array - - # define sin function - def sinfunc(t, A, w, p, c): - return A * np.sin(w*t + p) + c - - popt, pcov = scipy.optimize.curve_fit(sinfunc, t, y, p0=guess) # calculate sin function fit - A, w, p, c = popt # assign appropriate variables - fitfunc = lambda t: A * np.sin(w*t + p) + c # define fit function for curve fit - - # define function used to calculate least square - def sumfunc(x): - return sum((sinfunc(t, x[0], x[1], x[2], x[3]) - y)**2) - - sUx = scipy.optimize.minimize(fun=sumfunc, x0=np.array([guess_amp, 2.*np.pi*guess_freq, 0, guess_offset])) # calculate matching curve fit values with minimum error - return [{"amp": A, "omega": w, "phase": p, "offset": c, "fitfunc": fitfunc}, sUx] - - [Uu_fit, sUu] = fit_sin(t, Uu) # carry out calculations on self inductance - [Uv_fit, sUv] = fit_sin(t, Uv) # carry out calculations on mutual inductance - - fig1, ax1 = plt.subplots() - ax1.plot(t, Uu, "-k", label="y", linewidth=2) - ax1.plot(t, Uu_fit["fitfunc"](t), "r-", label="y fit curve", linewidth=2) - ax1.legend(loc="best") - plt.savefig("temp1.svg") - - fig2, ax2 = plt.subplots() - ax2.plot(t, Uv, "-k", label="y", linewidth=2) - ax2.plot(t, Uv_fit["fitfunc"](t), "r-", label="y fit curve", linewidth=2) - ax2.legend(loc="best") - plt.savefig("temp2.svg") - - Lzero = -2*sUv.x[3]/I_hat; # calculate L0 based on equations in publication - Lg = sUv.x[0]/I_hat # calculate Lg based on equations in publication - Lls = (sUu.x[3] + 2*sUv.x[3])/I_hat # calculate Lls based on equations in publication - Ld = Lls + 3/2*(Lzero - Lg) # calculate Ld based on equations in publication - Lq = Lls + 3/2*(Lzero + Lg) # calculate Lq based on equations in publication - saliency_ratio = Ld/Lq # calculate saliency ratio - - ############################ Output ################################# - post_processing = {} - post_processing["Ld"] = Ld - post_processing["Lq"] = Lq - post_processing["saliency_ratio"] = saliency_ratio - - state_out.conditions.inductance = post_processing - - print("\n************************ INDUCTANCE RESULTS ************************") - print("Ld = ", Ld, " H") - print("Lq = ", Lq, " H") - print("Saliency Ratio = ", saliency_ratio) - print("*************************************************************************\n") - - return state_out - -All example SynR evaluation scripts, including the one used for this analyzer, can be found in ``eMach\examples\mach_eval_examples\SynR_eval``, -where the post-analyzer script uses FEA results and calculates machine performance metrics, including torque density, power density, efficiency, -and torque ripple. This analyzer can be run by simply running the ``SynR_evaluator`` file in the aforementioned folder using the ``inductance_step``. -This example should produce the following results: - -.. csv-table:: `SynR_inductance_analyzer Results` - :file: results_SynR_inductance_analyzer.csv + inductance_prob = InductanceProblem(rated_current, linkages, rotor_angle, phasenames) + inductance_analyzer = Inductance_Analyzer(clarke_transformation_matrix) + data = inductance_analyzer.analyze(inductance_prob) + + rotor_angle = data["rotor_angle"] + Labc = data["Labc"] + Lalphabeta = data["Lalphabeta"] + Ldq = data["Ldq"] + L_d = np.mean(Ldq[:,0,0]) + L_q = np.mean(Ldq[:,1,1]) + saliency_ratio = L_d/L_q + + fig1 = plt.figure() + ax1 = plt.axes() + fig1.add_axes(ax1) + ax1.plot(rotor_angle, Labc[0,0,:]*1000) + ax1.plot(rotor_angle, Labc[1,1,:]*1000) + ax1.plot(rotor_angle, Labc[2,2,:]*1000) + ax1.set_xlabel("Rotor Angle [deg]") + ax1.set_ylabel("Inductance [mH]") + ax1.set_title("abc Inductances") + plt.legend(["$L_a$", "$L_b$", "$L_c$"], fontsize=12, loc='center right') + plt.grid(True, linewidth=0.5, color="#A9A9A9", linestyle="-.") + plt.show() + + fig2 = plt.figure() + ax2 = plt.axes() + fig2.add_axes(ax2) + ax2.plot(rotor_angle, Lalphabeta[:,0,0]*1000) + ax2.plot(rotor_angle, Lalphabeta[:,1,1]*1000) + ax2.plot(rotor_angle, Lalphabeta[:,2,2]*1000) + ax2.set_xlabel("Rotor Angle [deg]") + ax2.set_ylabel("Inductance [mH]") + ax2.set_title(r"$\alpha \beta \gamma$ Inductances") + plt.legend([r"$L_{\alpha \alpha}$", r"$L_{\beta \beta}$", r"$L_{\gamma \gamma}$"], fontsize=12, loc='center right') + plt.grid(True, linewidth=0.5, color="#A9A9A9", linestyle="-.") + plt.show() + + fig3 = plt.figure() + ax3 = plt.axes() + fig3.add_axes(ax3) + ax3.plot(rotor_angle, Ldq[:,0,0]*1000) + ax3.plot(rotor_angle, Ldq[:,1,1]*1000) + ax3.plot(rotor_angle, Ldq[:,2,2]*1000) + ax3.set_xlabel("Rotor Angle [deg]") + ax3.set_ylabel("Inductance [mH]") + ax3.set_title("dq0 Inductances") + plt.legend(["$L_d$", "$L_q$", "$L_0$"], fontsize=12, loc='center right') + plt.grid(True, linewidth=0.5, color="#A9A9A9", linestyle="-.") + plt.show() + + print("\n************************ INDUCTANCE RESULTS ************************") + print("Ld = ", L_d*1000, " mH") + print("Lq = ", L_q*1000, " mH") + print("Saliency Ratio = ", saliency_ratio) + print("*************************************************************************\n") + +This example, contained in the aforementioned ``fluxlinkage_inductance_eval`` folder, should produce the following results: + +.. figure:: ./Images/a_b_c_inductances.svg + :alt: a_b_c_inductances + :align: center + :width: 500 + +.. figure:: ./Images/alpha_beta_inductances.svg + :alt: alpha_beta_inductances + :align: center + :width: 500 + +.. figure:: ./Images/d_q_inductances.svg + :alt: d_q_inductances + :align: center + :width: 500 + +.. csv-table:: `inductance_analyzer Results` + :file: results_inductance_analyzer.csv :widths: 70, 70, 30 :header-rows: 1 + :align: center It should be noted that the inductance values calculated will be dependent on the number of turns in the stator. The saliency ratio however will -remain independent of this. \ No newline at end of file +remain independent of this. All of the code shown exists in the ``fluxlinkage_inductance_evaluator.py`` file in the +``eMach/examples/mach_eval_examples/fluxlinkage_inductance_eval`` folder. This analyzer serves as a second step in conjunction with the +`Flux Linkage Analyzer `_. \ No newline at end of file diff --git a/docs/source/EM_analyzers/init_SynR_inductance_analyzer.csv b/docs/source/EM_analyzers/init_SynR_inductance_analyzer.csv deleted file mode 100644 index 3d8a8757..00000000 --- a/docs/source/EM_analyzers/init_SynR_inductance_analyzer.csv +++ /dev/null @@ -1,2 +0,0 @@ -Arguments,Description -configuration,"object of type SynR_Inductance_Config describing time step, JMAG configurations, etc." diff --git a/docs/source/EM_analyzers/input_flux_linkage_analyzer.csv b/docs/source/EM_analyzers/input_flux_linkage_analyzer.csv new file mode 100644 index 00000000..0ebf67ab --- /dev/null +++ b/docs/source/EM_analyzers/input_flux_linkage_analyzer.csv @@ -0,0 +1,4 @@ +Input,Description,Units +tool,FEA application tool, +phase_names,list of names of phases, +rated_current,rated machine current,A diff --git a/docs/source/EM_analyzers/input_inductance_analyzer.csv b/docs/source/EM_analyzers/input_inductance_analyzer.csv new file mode 100644 index 00000000..065c0e90 --- /dev/null +++ b/docs/source/EM_analyzers/input_inductance_analyzer.csv @@ -0,0 +1,6 @@ +Arguments,Description,Units +I_hat,peak current,A +linkages,array of flux linkages, +rotor_angle,rotor angle array,deg +name_of_phases,name of phases dictionary, +clarke_trans_matrix,Clarke transformation matrix, diff --git a/docs/source/EM_analyzers/output_SynR_inductance_analyzer.csv b/docs/source/EM_analyzers/output_SynR_inductance_analyzer.csv deleted file mode 100644 index e664de85..00000000 --- a/docs/source/EM_analyzers/output_SynR_inductance_analyzer.csv +++ /dev/null @@ -1,4 +0,0 @@ -Output,Description -Ld,Scalar value of d-axis inductance -Lq,Scalar value of q-axis inductance -saliency_ratio,Ratio between Ld and Lq diff --git a/docs/source/EM_analyzers/output_flux_linkage_analyzer.csv b/docs/source/EM_analyzers/output_flux_linkage_analyzer.csv new file mode 100644 index 00000000..ebb115b2 --- /dev/null +++ b/docs/source/EM_analyzers/output_flux_linkage_analyzer.csv @@ -0,0 +1,5 @@ +Output,Description +linkages,group of flux linkage arrays +phase_currents,vector of phase currents in A +rotor_angle,rotor angle array in deg +name_of_phases,name of each machine phase diff --git a/docs/source/EM_analyzers/output_inductance_analyzer.csv b/docs/source/EM_analyzers/output_inductance_analyzer.csv new file mode 100644 index 00000000..4d350952 --- /dev/null +++ b/docs/source/EM_analyzers/output_inductance_analyzer.csv @@ -0,0 +1,5 @@ +Output,Description +rotor_angle,array of rotor angles +L_a_b_c,3D array of a-b-c inductances +L_alpha_beta,3D array of alpha-beta inductances +L_d_q,3D array of d-q inductances diff --git a/docs/source/EM_analyzers/results_SynR_inductance_analyzer.csv b/docs/source/EM_analyzers/results_SynR_inductance_analyzer.csv deleted file mode 100644 index 0c97ddf1..00000000 --- a/docs/source/EM_analyzers/results_SynR_inductance_analyzer.csv +++ /dev/null @@ -1,4 +0,0 @@ -Result,Value,Unit -Ld,0.0144,H -Lq,0.004,H -saliency_ratio,3.56, diff --git a/docs/source/EM_analyzers/results_flux_linkage_analyzer.csv b/docs/source/EM_analyzers/results_flux_linkage_analyzer.csv new file mode 100644 index 00000000..f313e82f --- /dev/null +++ b/docs/source/EM_analyzers/results_flux_linkage_analyzer.csv @@ -0,0 +1,5 @@ +Result,Value,Unit +linkages,Four 73x4 flux linkage arrays, +phase_currents,"[[20, 0, 0], [0, 20, 0], [0, 0, 20]]",A +rotor_angle,0 to 360 at 5 deg intervals,deg +name_of_phases,"['U','V','W']", diff --git a/docs/source/EM_analyzers/results_inductance_analyzer.csv b/docs/source/EM_analyzers/results_inductance_analyzer.csv new file mode 100644 index 00000000..8300a482 --- /dev/null +++ b/docs/source/EM_analyzers/results_inductance_analyzer.csv @@ -0,0 +1,4 @@ +Result,Value,Unit +Ld,14.4,mH +Lq,4.03,mH +saliency_ratio,3.57, diff --git a/examples/mach_eval_examples/SynR_eval/SynR_evaluator.py b/examples/mach_eval_examples/SynR_eval/SynR_evaluator.py index 6ce5cf36..217c636a 100644 --- a/examples/mach_eval_examples/SynR_eval/SynR_evaluator.py +++ b/examples/mach_eval_examples/SynR_eval/SynR_evaluator.py @@ -6,14 +6,15 @@ sys.path.append("../../../") from mach_eval import (MachineEvaluator, MachineDesign) -from electromagnetic_step import electromagnetic_step -from inductance_step import inductance_step +from SynR_flux_linkage_step import SynR_flux_linkage_step +from SynR_inductance_step import SynR_inductance_step from example_SynR_machine import Example_SynR_Machine, Machine_Op_Pt ############################ Create Evaluator ##################### SynR_evaluator = MachineEvaluator( [ - inductance_step + SynR_flux_linkage_step, + SynR_inductance_step ] ) diff --git a/examples/mach_eval_examples/SynR_eval/SynR_inductance_post_analyzer.py b/examples/mach_eval_examples/SynR_eval/SynR_inductance_post_analyzer.py deleted file mode 100644 index a24ef7de..00000000 --- a/examples/mach_eval_examples/SynR_eval/SynR_inductance_post_analyzer.py +++ /dev/null @@ -1,82 +0,0 @@ -import copy -import numpy as np -import matplotlib.pyplot as plt -import scipy.optimize - -class SynR_Inductance_PostAnalyzer: - - def get_next_state(results, in_state): - state_out = copy.deepcopy(in_state) - - ############################ Extract required info ########################### - flux_linkages = results["coil_flux_linkages"] - I_hat = results["current_peak"] - - ############################ post processing ########################### - data = flux_linkages.to_numpy() # change csv format to readable array - - t = data[:,0] # define x axis data as time - Uu = data[:,1] # define y axis data as self inductance - Uv = data[:,2] # define y axis data as mutual inductance - - # curve fit inductance values and calculate curve - def fit_sin(t, y): - fft_func = np.fft.fftfreq(len(t), (t[1]-t[0])) # define fft function with assumed uniform spacing - fft_y = abs(np.fft.fft(y)) # carry out fft function for inductance values - guess_freq = abs(fft_func[np.argmax(fft_y[1:])+1]) # excluding the zero frequency "peak", which can cause problematic fits - guess_amp = np.std(y) # guess amplitude based on one standard deviation - guess_offset = np.mean(y) # guess y offset based on average of magnitude - guess = np.array([guess_amp, 2.*np.pi*guess_freq, 0, guess_offset]) # arrage in array - - # define sin function - def sinfunc(t, A, w, p, c): - return A * np.sin(w*t + p) + c - - popt, pcov = scipy.optimize.curve_fit(sinfunc, t, y, p0=guess) # calculate sin function fit - A, w, p, c = popt # assign appropriate variables - fitfunc = lambda t: A * np.sin(w*t + p) + c # define fit function for curve fit - - # define function used to calculate least square - def sumfunc(x): - return sum((sinfunc(t, x[0], x[1], x[2], x[3]) - y)**2) - - sUx = scipy.optimize.minimize(fun=sumfunc, x0=np.array([guess_amp, 2.*np.pi*guess_freq, 0, guess_offset])) # calculate matching curve fit values with minimum error - return [{"amp": A, "omega": w, "phase": p, "offset": c, "fitfunc": fitfunc}, sUx] - - [Uu_fit, sUu] = fit_sin(t, Uu) # carry out calculations on self inductance - [Uv_fit, sUv] = fit_sin(t, Uv) # carry out calculations on mutual inductance - - fig1, ax1 = plt.subplots() - ax1.plot(t, Uu, "-k", label="y", linewidth=2) - ax1.plot(t, Uu_fit["fitfunc"](t), "r-", label="y fit curve", linewidth=2) - ax1.legend(loc="best") - #plt.savefig("temp1.svg") - - fig2, ax2 = plt.subplots() - ax2.plot(t, Uv, "-k", label="y", linewidth=2) - ax2.plot(t, Uv_fit["fitfunc"](t), "r-", label="y fit curve", linewidth=2) - ax2.legend(loc="best") - #plt.savefig("temp2.svg") - - Lzero = -2*sUv.x[3]/I_hat; # calculate L0 based on equations in publication - Lg = sUv.x[0]/I_hat # calculate Lg based on equations in publication - Lls = (sUu.x[3] + 2*sUv.x[3])/I_hat # calculate Lls based on equations in publication - Ld = Lls + 3/2*(Lzero - Lg) # calculate Ld based on equations in publication - Lq = Lls + 3/2*(Lzero + Lg) # calculate Lq based on equations in publication - saliency_ratio = Ld/Lq # calculate saliency ratio - - ############################ Output ################################# - post_processing = {} - post_processing["Ld"] = Ld - post_processing["Lq"] = Lq - post_processing["saliency_ratio"] = saliency_ratio - - state_out.conditions.inductance = post_processing - - print("\n************************ INDUCTANCE RESULTS ************************") - print("Ld = ", Ld, " H") - print("Lq = ", Lq, " H") - print("Saliency Ratio = ", saliency_ratio) - print("*************************************************************************\n") - - return state_out \ No newline at end of file diff --git a/examples/mach_eval_examples/SynR_eval/example_SynR_machine.py b/examples/mach_eval_examples/SynR_eval/example_SynR_machine.py index ad284940..79c17aeb 100644 --- a/examples/mach_eval_examples/SynR_eval/example_SynR_machine.py +++ b/examples/mach_eval_examples/SynR_eval/example_SynR_machine.py @@ -62,6 +62,7 @@ SynR_winding = { "no_of_layers": 2, + "name_of_phases": ['U', 'V', 'W'], "layer_phases": [ ['U', 'W', 'V', 'U', 'W', 'V', 'U', 'W', 'V', 'U', 'W', 'V'], ['W', 'V', 'U', 'W', 'V', 'U', 'W', 'V', 'U', 'W', 'V', 'U'] ], "layer_polarity": [ ['+', '-', '+', '-', '+', '-', '+', '-', '+', '-', '+', '-'], diff --git a/examples/mach_eval_examples/SynR_eval/inductance_step.py b/examples/mach_eval_examples/SynR_eval/inductance_step.py deleted file mode 100644 index a6d33c7a..00000000 --- a/examples/mach_eval_examples/SynR_eval/inductance_step.py +++ /dev/null @@ -1,52 +0,0 @@ -import os -import sys -import copy - -from mach_eval import AnalysisStep, ProblemDefinition -from mach_eval.analyzers.electromagnetic.SynR import SynR_inductance_analyzer as SynR_inductance -from mach_eval.analyzers.electromagnetic.SynR.SynR_em_config import SynR_EM_Config -from examples.mach_eval_examples.SynR_eval.SynR_inductance_post_analyzer import SynR_Inductance_PostAnalyzer - -############################ Define Electromagnetic Step ########################### -class SynR_EM_ProblemDefinition(ProblemDefinition): - """Converts a State into a problem""" - - def __init__(self): - pass - - def get_problem(state): - - problem = SynR_inductance.SynR_Inductance_Problem( - state.design.machine, state.design.settings) - return problem - -# initialize em analyzer class with FEA configuration -configuration = SynR_EM_Config( - no_of_rev = 1, - no_of_steps = 72, - - mesh_size=3, # mm - mesh_size_rotor=1.5, # mm - airgap_mesh_radial_div=4, - airgap_mesh_circum_div=720, - mesh_air_region_scale=1.05, - - only_table_results=False, - csv_results=("FEMCoilFlux"), - del_results_after_calc=False, - run_folder=os.path.dirname(__file__) + "/run_data/", - jmag_csv_folder=os.path.dirname(__file__) + "/run_data/jmag_csv/", - - max_nonlinear_iterations=50, - multiple_cpus=True, - num_cpus=4, - jmag_scheduler=False, - jmag_visible=True, - non_zero_end_ring_res = False, - scale_axial_length = True, - jmag_version="21.1", -) - -SynR_inductance_analysis = SynR_inductance.SynR_Inductance_Analyzer(configuration) - -inductance_step = AnalysisStep(SynR_EM_ProblemDefinition, SynR_inductance_analysis, SynR_Inductance_PostAnalyzer) \ No newline at end of file diff --git a/examples/mach_eval_examples/fluxlinkage_inductance_eval/Example_FluxLinkage_Machine.jproj b/examples/mach_eval_examples/fluxlinkage_inductance_eval/Example_FluxLinkage_Machine.jproj new file mode 100644 index 00000000..e6e123da Binary files /dev/null and b/examples/mach_eval_examples/fluxlinkage_inductance_eval/Example_FluxLinkage_Machine.jproj differ diff --git a/examples/mach_eval_examples/fluxlinkage_inductance_eval/fluxlinkage_inductance_evaluator.py b/examples/mach_eval_examples/fluxlinkage_inductance_eval/fluxlinkage_inductance_evaluator.py new file mode 100644 index 00000000..aad16098 --- /dev/null +++ b/examples/mach_eval_examples/fluxlinkage_inductance_eval/fluxlinkage_inductance_evaluator.py @@ -0,0 +1,105 @@ +import os +import sys +import numpy as np +import matplotlib.pyplot as plt +from time import time as clock_time + +os.chdir(os.path.dirname(__file__)) +sys.path.append("../../../") + +from mach_eval.analyzers.electromagnetic.flux_linkage_analyzer import FluxLinkageJMAG_Problem, FluxLinkageJMAG_Analyzer +from mach_eval.analyzers.electromagnetic.inductance_analyzer import InductanceProblem, InductanceAnalyzer + +from mach_cad.tools import jmag as JMAG + +script_dir = os.path.dirname(__file__) +filename = os.path.join(script_dir, "Example_FluxLinkage_Machine.jproj") +phasenames = ['U', 'V', 'W'] +rated_current = 20 + +#################################################### +# 01 Setting project name and output folder +#################################################### + +toolJmag = JMAG.JmagDesigner() +toolJmag.visible = True +toolJmag.open(filename) +toolJmag.save() + +############################ Create Evaluator ##################### +tic = clock_time() +flux_linkage_prob = FluxLinkageJMAG_Problem(toolJmag, phasenames, rated_current) +flux_linkage_analyzer = FluxLinkageJMAG_Analyzer() +flux_linkages, currents = flux_linkage_analyzer.analyze(flux_linkage_prob) +toc = clock_time() +print("Time spent on the flux linkage evaluation is %g min." % ((toc- tic)/60)) + +linkages = flux_linkages["linkages"] +phase_currents = currents +rotor_angle = flux_linkages["rotor_angle"][0] +name_of_phases = flux_linkages["name_of_phases"] + +print("\n************************ FLUX LINKAGE RESULTS ************************") +print("linkages = ", linkages) +print("phase_currents = ", phase_currents, " A") +print("rotor_angle = ", rotor_angle, " deg") +print("name_of_phases = ", name_of_phases) +print("*************************************************************************\n") + +clarke_transformation_matrix = 2/3*np.array([[1, -1/2, -1/2], [0, np.sqrt(3)/2, -np.sqrt(3)/2], [1/2, 1/2, 1/2]]) +inductance_prob = InductanceProblem(rated_current, linkages, rotor_angle, phasenames) +inductance_analyzer = InductanceAnalyzer(clarke_transformation_matrix) +data = inductance_analyzer.analyze(inductance_prob) + +rotor_angle = data["rotor_angle"] +Labc = data["Labc"] +Lalphabeta = data["Lalphabeta"] +Ldq = data["Ldq"] +L_d = np.mean(Ldq[:,0,0]) +L_q = np.mean(Ldq[:,1,1]) +saliency_ratio = L_d/L_q + +fig1 = plt.figure() +ax1 = plt.axes() +fig1.add_axes(ax1) +ax1.plot(rotor_angle, Labc[0,0,:]*1000) +ax1.plot(rotor_angle, Labc[1,1,:]*1000) +ax1.plot(rotor_angle, Labc[2,2,:]*1000) +ax1.set_xlabel("Rotor Angle [deg]") +ax1.set_ylabel("Inductance [mH]") +ax1.set_title("abc Inductances") +plt.legend(["$L_a$", "$L_b$", "$L_c$"], fontsize=12, loc='center right') +plt.grid(True, linewidth=0.5, color="#A9A9A9", linestyle="-.") +plt.show() + +fig2 = plt.figure() +ax2 = plt.axes() +fig2.add_axes(ax2) +ax2.plot(rotor_angle, Lalphabeta[:,0,0]*1000) +ax2.plot(rotor_angle, Lalphabeta[:,1,1]*1000) +ax2.plot(rotor_angle, Lalphabeta[:,2,2]*1000) +ax2.set_xlabel("Rotor Angle [deg]") +ax2.set_ylabel("Inductance [mH]") +ax2.set_title(r"$\alpha \beta \gamma$ Inductances") +plt.legend([r"$L_{\alpha \alpha}$", r"$L_{\beta \beta}$", r"$L_{\gamma \gamma}$"], fontsize=12, loc='center right') +plt.grid(True, linewidth=0.5, color="#A9A9A9", linestyle="-.") +plt.show() + +fig3 = plt.figure() +ax3 = plt.axes() +fig3.add_axes(ax3) +ax3.plot(rotor_angle, Ldq[:,0,0]*1000) +ax3.plot(rotor_angle, Ldq[:,1,1]*1000) +ax3.plot(rotor_angle, Ldq[:,2,2]*1000) +ax3.set_xlabel("Rotor Angle [deg]") +ax3.set_ylabel("Inductance [mH]") +ax3.set_title("dq0 Inductances") +plt.legend(["$L_d$", "$L_q$", "$L_0$"], fontsize=12, loc='center right') +plt.grid(True, linewidth=0.5, color="#A9A9A9", linestyle="-.") +plt.show() + +print("\n************************ INDUCTANCE RESULTS ************************") +print("Ld = ", L_d*1000, " mH") +print("Lq = ", L_q*1000, " mH") +print("Saliency Ratio = ", saliency_ratio) +print("*************************************************************************\n") \ No newline at end of file diff --git a/mach_eval/analyzers/electromagnetic/SynR/SynR_em_analyzer.py b/mach_eval/analyzers/electromagnetic/SynR/SynR_em_analyzer.py index 8774c2a5..d29e7ed1 100644 --- a/mach_eval/analyzers/electromagnetic/SynR/SynR_em_analyzer.py +++ b/mach_eval/analyzers/electromagnetic/SynR/SynR_em_analyzer.py @@ -118,7 +118,7 @@ def analyze(self, problem): model.SetName(self.project_name) model.SetDescription(self.show(self.project_name, toString=True)) - valid_design = self.pre_process(model) + valid_design = self.pre_process(app, model) if not valid_design: raise InvalidDesign @@ -403,7 +403,7 @@ def show(self, name, toString=False): "%s = %s" % item for item in tuple_list ) - def pre_process(self, model): + def pre_process(self, app, model): # pre-process : you can select part by coordinate! """Group""" @@ -427,6 +427,13 @@ def group(name, id_list): group("Coils", partIDRange_Coil) + # """ Set Parts names """ + + app.GetModel(0).SetPartName(id_rotorCore, u"RotorCore") + app.GetModel(0).SetPartName(id_shaft, u"Shaft") + app.GetModel(0).SetPartName(id_statorCore, u"StatorCore") + + """ Add Part to Set for later references """ def add_part_to_set(name, x, y, ID=None): diff --git a/mach_eval/analyzers/electromagnetic/SynR/SynR_inductance_analyzer.py b/mach_eval/analyzers/electromagnetic/SynR/SynR_inductance_analyzer.py deleted file mode 100644 index 241c9670..00000000 --- a/mach_eval/analyzers/electromagnetic/SynR/SynR_inductance_analyzer.py +++ /dev/null @@ -1,1056 +0,0 @@ -import os -import numpy as np -import pandas as pd -import sys -from time import time as clock_time - -from mach_cad import model_obj as mo -from mach_opt import InvalidDesign -from mach_eval.analyzers.electromagnetic.stator_wdg_res import( - StatorWindingResistanceProblem, StatorWindingResistanceAnalyzer -) -from mach_cad.tools import jmag as JMAG - -class SynR_Inductance_Problem: - def __init__(self, machine, operating_point): - self.machine = machine - self.operating_point = operating_point - self._validate_attr() - - def _validate_attr(self): - if 'SynR_Machine' in str(type(self.machine)): - pass - else: - raise TypeError("Invalid machine type") - - if 'SynR_Machine_Oper_Pt' in str(type(self.operating_point)): - pass - else: - raise TypeError("Invalid settings type") - - -class SynR_Inductance_Analyzer: - def __init__(self, configuration): - self.config = configuration - - def analyze(self, problem): - self.machine_variant = problem.machine - self.operating_point = problem.operating_point - - #################################################### - # 01 Setting project name and output folder - #################################################### - - self.project_name = self.machine_variant.name - expected_project_file = self.config.run_folder + "_Ind" + "%s.jproj" % self.project_name - - # Create output folder - if not os.path.isdir(self.config.jmag_csv_folder): - os.makedirs(self.config.jmag_csv_folder) - - attempts = 1 - if os.path.exists(expected_project_file): - print( - "JMAG project exists already, I will not delete it but create a new one with a different name instead." - ) - attempts = 2 - temp_path = expected_project_file[ - : -len(".jproj") - ] + "_attempts_%d.jproj" % (attempts) - while os.path.exists(temp_path): - attempts += 1 - temp_path = expected_project_file[ - : -len(".jproj") - ] + "_attempts_%d.jproj" % (attempts) - - expected_project_file = temp_path - - if attempts > 1: - self.project_name = self.project_name + "_attempts_%d" % (attempts) - - JMAG_exist = os.path.exists("C:/Program Files/JMAG-Designer%s" % self.config.jmag_version) - if JMAG_exist is True: - toolJmag = JMAG.JmagDesigner(self.config.jmag_version) - else: - print('WARNING: The JMAG version does not exist, defaulting to newest JMAG version available!') - toolJmag = JMAG.JmagDesigner() - - toolJmag.visible = self.config.jmag_visible - toolJmag.open(comp_filepath=expected_project_file, length_unit="DimMillimeter", study_type="Transient2D") - toolJmag.save() - - self.study_name = self.project_name + "_Ind_SynR" - self.design_results_folder = ( - self.config.run_folder + "%s_results/" % self.project_name - ) - if not os.path.isdir(self.design_results_folder): - os.makedirs(self.design_results_folder) - - ################################################################ - # 02 Run Electromagnetic analysis - ################################################################ - - # Draw cross_section - draw_success = self.draw_machine(toolJmag) - - if not draw_success: - raise InvalidDesign - - toolJmag.doc.SaveModel(False) - app = toolJmag.jd - model = app.GetCurrentModel() - - # Pre-processing - model.SetName(self.project_name) - model.SetDescription(self.show(self.project_name, toString=True)) - - valid_design = self.pre_process(model) - - if not valid_design: - raise InvalidDesign - - # Create transient study with two time step sections - study = self.add_em_study(app, model, self.config.jmag_csv_folder, self.study_name) - self.create_custom_material( - app, self.machine_variant.stator_iron_mat["core_material"] - ) - app.SetCurrentStudy(self.study_name) - - # Mesh study - self.mesh_study(app, model, study) - - self.breakdown_torque = None - self.stator_slot_area = None - - study.GetDesignTable().GetEquation("freq").SetExpression("%g" % self.drive_freq) - - # Set current excitation - I = self.I_hat - phi_0 = self.operating_point.phi_0 - self.set_currents_sequence(I, 0.0001, - phi_0, app, study) - - # Add time step settings - no_of_steps = self.config.no_of_steps - no_of_rev = self.config.no_of_rev - time_interval = no_of_rev / (self.drive_freq) - self.add_time_step_settings(time_interval, no_of_steps, app, study) - - self.run_study(app, study, clock_time()) - - toolJmag.save() - app.Quit() - - #################################################### - # 03 Load FEA output - #################################################### - - fea_rated_output = self.extract_JMAG_results( - self.config.jmag_csv_folder, self.study_name - ) - - return fea_rated_output - - @property - def drive_freq(self): - speed_in_elec_ang = (self.operating_point.speed * self.operating_point.speed_ratio) / 60 * self.machine_variant.p - drive_freq = speed_in_elec_ang - return drive_freq - - @property - def speed(self): - return (self.operating_point.speed * self.operating_point.speed_ratio) - - @property - def elec_omega(self): - return 2 * np.pi * self.drive_freq - - @property - def I_hat(self): - I_hat = self.machine_variant.rated_current - return I_hat - - @property - def phi_0(self): - return self.operating_point.phi_0 - - @property - def z_C(self): - if len(self.machine_variant.layer_phases) == 1: - z_C = self.machine_variant.Q / 3 - elif len(self.machine_variant.layer_phases) == 2: - z_C = self.machine_variant.Q / 3 - - return z_C - - @property - def stator_resistance(self): - res_prob = StatorWindingResistanceProblem( - r_si=self.machine_variant.r_si/1000, - d_sp=self.machine_variant.d_sp/1000, - d_st=self.machine_variant.d_st/1000, - w_st=self.machine_variant.w_st/1000, - l_st=self.machine_variant.l_st/1000, - Q=self.machine_variant.Q, - y=self.machine_variant.pitch, - z_Q=self.machine_variant.Z_q, - z_C=self.z_C, - Kcu=self.machine_variant.Kcu, - Kov=self.machine_variant.Kov, - sigma_cond=self.machine_variant.coil_mat["copper_elec_conductivity"], - slot_area=self.machine_variant.s_slot*1e-6, - n_layers=self.machine_variant.no_of_layers, - ) - res_analyzer = StatorWindingResistanceAnalyzer() - stator_resistance = res_analyzer.analyze(res_prob) - return stator_resistance - - @property - def R_wdg(self): - return self.stator_resistance["R_wdg"] - - @property - def R_wdg_coil_ends(self): - return 2 * self.stator_resistance["R_ew"] * self.z_C - - @property - def R_wdg_coil_sides(self): - return self.R_wdg - self. R_wdg_coil_ends - - - def draw_machine(self, tool): - #################################################### - # Adding parts objects - #################################################### - self.stator_core = mo.CrossSectInnerRotorStatorPartial( - name="StatorCore", - dim_alpha_st=mo.DimDegree(self.machine_variant.alpha_st), - dim_alpha_so=mo.DimDegree(self.machine_variant.alpha_so), - dim_r_si=mo.DimMillimeter(self.machine_variant.r_si), - dim_d_so=mo.DimMillimeter(self.machine_variant.d_so), - dim_d_sp=mo.DimMillimeter(self.machine_variant.d_sp), - dim_d_st=mo.DimMillimeter(self.machine_variant.d_st), - dim_d_sy=mo.DimMillimeter(self.machine_variant.d_sy), - dim_w_st=mo.DimMillimeter(self.machine_variant.w_st), - dim_r_st=mo.DimMillimeter(0), - dim_r_sf=mo.DimMillimeter(0), - dim_r_sb=mo.DimMillimeter(0), - Q=self.machine_variant.Q, - location=mo.Location2D(anchor_xy=[mo.DimMillimeter(0), mo.DimMillimeter(0)], - theta=mo.DimDegree(-180 / self.machine_variant.Q)), - ) - - self.winding_layer1 = mo.CrossSectInnerRotorStatorRightSlot( - name="WindingLayer1", - stator_core=self.stator_core, - location=mo.Location2D(anchor_xy=[mo.DimMillimeter(0), mo.DimMillimeter(0)], - theta=mo.DimDegree(-180 / self.machine_variant.Q)), - ) - - self.winding_layer2 = mo.CrossSectInnerRotorStatorLeftSlot( - name="WindingLayer2", - stator_core=self.stator_core, - location=mo.Location2D(anchor_xy=[mo.DimMillimeter(0), mo.DimMillimeter(0)], - theta=mo.DimDegree(-180 / self.machine_variant.Q)), - ) - - self.rotor_core = mo.CrossSectFluxBarrierRotorPartial( - name="RotorCore", - dim_alpha_b=mo.DimDegree(self.machine_variant.alpha_b), - dim_r_ri=mo.DimMillimeter(self.machine_variant.r_ri), - dim_r_ro=mo.DimMillimeter(self.machine_variant.r_ro), - dim_r_f1=mo.DimMillimeter(self.machine_variant.r_f1), - dim_r_f2=mo.DimMillimeter(self.machine_variant.r_f2), - dim_r_f3=mo.DimMillimeter(self.machine_variant.r_f3), - dim_d_r1=mo.DimMillimeter(self.machine_variant.d_r1), - dim_d_r2=mo.DimMillimeter(self.machine_variant.d_r2), - dim_d_r3=mo.DimMillimeter(self.machine_variant.d_r3), - dim_w_b1=mo.DimMillimeter(self.machine_variant.w_b1), - dim_w_b2=mo.DimMillimeter(self.machine_variant.w_b2), - dim_w_b3=mo.DimMillimeter(self.machine_variant.w_b3), - dim_l_b1=mo.DimMillimeter(self.machine_variant.l_b1), - dim_l_b2=mo.DimMillimeter(self.machine_variant.l_b2), - dim_l_b3=mo.DimMillimeter(self.machine_variant.l_b3), - dim_l_b4=mo.DimMillimeter(self.machine_variant.l_b4), - dim_l_b5=mo.DimMillimeter(self.machine_variant.l_b5), - dim_l_b6=mo.DimMillimeter(self.machine_variant.l_b6), - p=2, - location=mo.Location2D(anchor_xy=[mo.DimMillimeter(0), mo.DimMillimeter(0)], theta=mo.DimDegree(-180 / (2 * self.machine_variant.p))), - ) - - self.shaft = mo.CrossSectHollowCylinder( - name="Shaft", - dim_t=mo.DimMillimeter(self.machine_variant.r_ri), - dim_r_o=mo.DimMillimeter(self.machine_variant.r_ri), - location=mo.Location2D(anchor_xy=[mo.DimMillimeter(0), mo.DimMillimeter(0)]), - ) - - - self.comp_stator_core = mo.Component( - name="StatorCore", - cross_sections=[self.stator_core], - material=mo.MaterialGeneric(name=self.machine_variant.stator_iron_mat["core_material"], color=r"#808080"), - make_solid=mo.MakeExtrude(location=mo.Location3D(), - dim_depth=mo.DimMillimeter(self.machine_variant.l_st)), - ) - - self.comp_winding_layer1 = mo.Component( - name="WindingLayer1", - cross_sections=[self.winding_layer1], - material=mo.MaterialGeneric(name=self.machine_variant.coil_mat["coil_material"]), - make_solid=mo.MakeExtrude(location=mo.Location3D(), - dim_depth=mo.DimMillimeter(self.machine_variant.l_st)), - ) - - self.comp_winding_layer2 = mo.Component( - name="WindingLayer2", - cross_sections=[self.winding_layer2], - material=mo.MaterialGeneric(name=self.machine_variant.coil_mat["coil_material"]), - make_solid=mo.MakeExtrude(location=mo.Location3D(), - dim_depth=mo.DimMillimeter(self.machine_variant.l_st)), - ) - - self.comp_rotor_core = mo.Component( - name="RotorCore", - cross_sections=[self.rotor_core], - material=mo.MaterialGeneric(name=self.machine_variant.rotor_iron_mat["core_material"], color=r"#808080"), - make_solid=mo.MakeExtrude(location=mo.Location3D(), - dim_depth=mo.DimMillimeter(self.machine_variant.l_st)), - ) - - self.comp_shaft = mo.Component( - name="Shaft", - cross_sections=[self.shaft], - material=mo.MaterialGeneric(name=self.machine_variant.shaft_mat["shaft_material"], color=r"#71797E"), - make_solid=mo.MakeExtrude(location=mo.Location3D(), - dim_depth=mo.DimMillimeter(self.machine_variant.l_st)), - ) - - tool.bMirror = False - - tool.sketch = tool.create_sketch() - tool.sketch.SetProperty("Name", self.stator_core.name) - tool.sketch.SetProperty("Color", r"#808080") - cs_stator = self.stator_core.draw(tool) - stator_tool = tool.prepare_section(cs_stator, self.machine_variant.Q) - - tool.sketch = tool.create_sketch() - tool.sketch.SetProperty("Name", self.winding_layer1.name) - tool.sketch.SetProperty("Color", r"#B87333") - cs_winding_layer1 = self.winding_layer1.draw(tool) - winding_tool1 = tool.prepare_section(cs_winding_layer1, self.machine_variant.Q) - self.winding_layer1_inner_coord = cs_winding_layer1.inner_coord - - tool.sketch = tool.create_sketch() - tool.sketch.SetProperty("Name", self.winding_layer2.name) - tool.sketch.SetProperty("Color", r"#B87333") - cs_winding_layer2 = self.winding_layer2.draw(tool) - winding_tool2 = tool.prepare_section(cs_winding_layer2, self.machine_variant.Q) - self.winding_layer2_inner_coord = cs_winding_layer2.inner_coord - - tool.sketch = tool.create_sketch() - tool.sketch.SetProperty("Name", self.rotor_core.name) - tool.sketch.SetProperty("Color", r"#808080") - cs_rotor_core = self.rotor_core.draw(tool) - rotor_tool = tool.prepare_section(cs_rotor_core, 2*self.machine_variant.p) - - tool.sketch = tool.create_sketch() - tool.sketch.SetProperty("Name", self.shaft.name) - tool.sketch.SetProperty("Color", r"#71797E") - cs_shaft = self.shaft.draw(tool) - shaft_tool = tool.prepare_section(cs_shaft) - - return True - - - def show(self, name, toString=False): - attrs = list(vars(self).items()) - key_list = [el[0] for el in attrs] - val_list = [el[1] for el in attrs] - the_dict = dict(list(zip(key_list, val_list))) - sorted_key = sorted( - key_list, - key=lambda item: ( - int(item.partition(" ")[0]) if item[0].isdigit() else float("inf"), - item, - ), - ) # this is also useful for string beginning with digiterations '15 Steel'. - tuple_list = [(key, the_dict[key]) for key in sorted_key] - if not toString: - print("- SynR Individual #%s\n\t" % name, end=" ") - print(", \n\t".join("%s = %s" % item for item in tuple_list)) - return "" - else: - return "\n- SynR Individual #%s\n\t" % name + ", \n\t".join( - "%s = %s" % item for item in tuple_list - ) - - def pre_process(self, model): - # pre-process : you can select part by coordinate! - """Group""" - - def group(name, id_list): - model.GetGroupList().CreateGroup(name) - for the_id in id_list: - model.GetGroupList().AddPartToGroup(name, the_id) - - part_ID_list = model.GetPartIDs() - - if len(part_ID_list) != int( - 1 + 1 + self.machine_variant.Q * 2 + 1 - ): - print("Parts are missing in this machine") - return False - - self.id_statorCore = id_statorCore = part_ID_list[0] - partIDRange_Coil = part_ID_list[1 : int(2 * self.machine_variant.Q + 1)] - self.id_rotorCore = id_rotorCore = part_ID_list[int(2 * self.machine_variant.Q + 1)] - id_shaft = part_ID_list[-1] - - group("Coils", partIDRange_Coil) - - """ Add Part to Set for later references """ - - def add_part_to_set(name, x, y, ID=None): - model.GetSetList().CreatePartSet(name) - model.GetSetList().GetSet(name).SetMatcherType("Selection") - model.GetSetList().GetSet(name).ClearParts() - sel = model.GetSetList().GetSet(name).GetSelection() - if ID is None: - # print x,y - sel.SelectPartByPosition(x, y, 0) # z=0 for 2D - else: - sel.SelectPart(ID) - model.GetSetList().GetSet(name).AddSelected(sel) - - # RotorSet - add_part_to_set("RotorSet", 0.0, 0.0, ID=id_shaft) - - # Create Set for right layer - Angle_StatorSlotSpan = 360 / self.machine_variant.Q - R = np.sqrt(self.winding_layer1_inner_coord[0] ** 2 + self.winding_layer1_inner_coord[1] ** 2) - THETA = np.arctan(self.winding_layer2_inner_coord[1] / self.winding_layer2_inner_coord[0]) - X = R * np.cos(THETA) - Y = R * np.sin(THETA) - count = 0 - for UVW, UpDown in zip( - self.machine_variant.layer_phases[0], self.machine_variant.layer_polarity[0] - ): - count += 1 - add_part_to_set("coil_right_%s%s %d" % (UVW, UpDown, count), X, Y) - - THETA += Angle_StatorSlotSpan / 180.0 * np.pi - X = R * np.cos(THETA) - Y = R * np.sin(THETA) - - # Create Set for left layer - THETA = np.arctan(self.winding_layer1_inner_coord[1] / self.winding_layer1_inner_coord[0]) - X = R * np.cos(THETA) - Y = R * np.sin(THETA) - count = 0 - for UVW, UpDown in zip( - self.machine_variant.layer_phases[1], self.machine_variant.layer_polarity[1] - ): - count += 1 - add_part_to_set("coil_left_%s%s %d" % (UVW, UpDown, count), X, Y) - - THETA += Angle_StatorSlotSpan / 180.0 * np.pi - X = R * np.cos(THETA) - Y = R * np.sin(THETA) - - # Create Set for Motion Region - def part_list_set(name, list_part_id=None, prefix=None): - model.GetSetList().CreatePartSet(name) - model.GetSetList().GetSet(name).SetMatcherType("Selection") - model.GetSetList().GetSet(name).ClearParts() - sel = model.GetSetList().GetSet(name).GetSelection() - - if list_part_id is not None: - for ID in list_part_id: - sel.SelectPart(ID) - model.GetSetList().GetSet(name).AddSelected(sel) - - part_list_set( - "Motion_Region", list_part_id=[id_rotorCore, id_shaft] - ) - - return True - - - def create_custom_material(self, app, steel_name): - - core_mat_obj = app.GetMaterialLibrary().GetCustomMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ) - app.GetMaterialLibrary().DeleteCustomMaterialByObject(core_mat_obj) - - app.GetMaterialLibrary().CreateCustomMaterial( - self.machine_variant.stator_iron_mat["core_material"], "Custom Materials" - ) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue( - "Density", self.machine_variant.stator_iron_mat["core_material_density"] / 1000 - ) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("MagneticSteelPermeabilityType", 2) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("CoerciveForce", 0) - BH = np.loadtxt( - self.machine_variant.stator_iron_mat["core_bh_file"], - unpack=True, - usecols=(0, 1), - ) # values from Nishanth Magnet BH curve - refarray = BH.T.tolist() - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).GetTable("BhTable").SetTable(refarray) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("DemagnetizationCoerciveForce", 0) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("MagnetizationSaturated", 0) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("MagnetizationSaturated2", 0) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("ExtrapolationMethod", 0) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue( - "YoungModulus", self.machine_variant.stator_iron_mat["core_youngs_modulus"] / 1000000 - ) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue("Loss_Type", 1) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue( - "LossConstantKhX", self.machine_variant.stator_iron_mat["core_ironloss_Kh"] - ) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue( - "LossConstantKeX", self.machine_variant.stator_iron_mat["core_ironloss_Ke"] - ) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue( - "LossConstantAlphaX", - self.machine_variant.stator_iron_mat["core_ironloss_a"], - ) - app.GetMaterialLibrary().GetUserMaterial( - self.machine_variant.stator_iron_mat["core_material"] - ).SetValue( - "LossConstantBetaX", self.machine_variant.stator_iron_mat["core_ironloss_b"] - ) - - - - def add_em_study( - self, app, model, dir_csv_output_folder, study_name - ): - - model.CreateStudy("Transient2D", study_name) - app.SetCurrentStudy(self.study_name) - study = model.GetStudy(study_name) - - # Study properties - study.GetStudyProperties().SetValue("ApproximateTransientAnalysis", 1) # psuedo steady state freq is for PWM drive to use - study.GetStudyProperties().SetValue("OutputSteadyResultAs1stStep", 0) - study.GetStudyProperties().SetValue("ConversionType", 0) - study.GetStudyProperties().SetValue( - "NonlinearMaxIteration", self.config.max_nonlinear_iterations - ) - study.GetStudyProperties().SetValue( - "ModelThickness", self.machine_variant.l_st - ) # [mm] Stack Length - - # Material - self.add_materials(study) - - # Conditions - Motion - study.CreateCondition("RotationMotion", - "RotCon") - study.GetCondition("RotCon").SetValue("AngularVelocity", - int(self.speed)) - study.GetCondition("RotCon").ClearParts() - study.GetCondition("RotCon").AddSet( - model.GetSetList().GetSet("Motion_Region"), 0 - ) - - study.CreateCondition( - "Torque", "TorCon" - ) - study.GetCondition("TorCon").SetValue("TargetType", 1) - study.GetCondition("TorCon").SetLinkWithType("LinkedMotion", "RotCon") - study.GetCondition("TorCon").ClearParts() - - study.CreateCondition("Force", "ForCon") - study.GetCondition("ForCon").SetValue("TargetType", 1) - study.GetCondition("ForCon").SetLinkWithType("LinkedMotion", "RotCon") - study.GetCondition("ForCon").ClearParts() - - # Conditions - FEM Coils & Conductors (i.e. stator/rotor winding) - self.add_circuit(app, model, study, bool_3PhaseCurrentSource=True) - - # True: no mesh or field results are needed - study.GetStudyProperties().SetValue( - "OnlyTableResults", self.config.only_table_results - ) - - study.GetStudyProperties().SetValue("DirectSolverType", 1) - - if self.config.multiple_cpus: - # This SMP(shared memory process) is effective only if there are tons of elements. e.g., over 100,000. - # too many threads will in turn make them compete with each other and slow down the solve. 2 is good enough - # for eddy current solve. 6~8 is enough for transient solve. - study.GetStudyProperties().SetValue("UseMultiCPU", True) - study.GetStudyProperties().SetValue("MultiCPU", self.config.num_cpus) - - # two sections of different time step - no_of_rev = self.config.no_of_rev - no_of_steps = self.config.no_of_steps - number_of_total_steps = ( - 1 + no_of_steps - ) - # add equations - study.GetDesignTable().AddEquation("freq") - study.GetDesignTable().AddEquation("speed") - study.GetDesignTable().GetEquation("freq").SetType(0) - study.GetDesignTable().GetEquation("freq").SetDescription( - "Excitation Frequency" - ) - study.GetDesignTable().GetEquation("speed").SetType(1) - study.GetDesignTable().GetEquation("speed").SetExpression("freq * %d"%(60/(self.machine_variant.p))) - study.GetDesignTable().GetEquation("speed").SetDescription( - "mechanical speed of rotor" - ) - - # speed, freq - study.GetCondition("RotCon").SetValue("AngularVelocity", "speed") - study.GetStudyProperties().SetValue("ApproximateTransientAnalysis", 1) - study.GetStudyProperties().SetValue("OutputSteadyResultAs1stStep", 0) - - # Iron Loss Calculation Condition - # Stator - if True: - cond = study.CreateCondition("Ironloss", "IronLossConStator") - cond.SetValue("RevolutionSpeed", "freq*60/%d" % self.machine_variant.p) - cond.ClearParts() - sel = cond.GetSelection() - sel.SelectPart(self.id_statorCore) - cond.AddSelected(sel) - # Use FFT for hysteresis to be consistent with FEMM's results and to have a FFT plot - cond.SetValue("HysteresisLossCalcType", 1) - cond.SetValue("PresetType", 3) # 3:Custom - # Specify the reference steps yourself because you don't really know what JMAG is doing behind you - cond.SetValue( - "StartReferenceStep", - number_of_total_steps + 1 - no_of_steps / no_of_rev * 0.25, - ) # 1/4 period in no of steps = no_of_steps / no_of_rev * 0.25 - cond.SetValue("EndReferenceStep", number_of_total_steps) - cond.SetValue("UseStartReferenceStep", 1) - cond.SetValue("UseEndReferenceStep", 1) - cond.SetValue( - "Cyclicity", 4 - ) # specify reference steps for 1/4 period and extend it to whole period - cond.SetValue("UseFrequencyOrder", 1) - cond.SetValue("FrequencyOrder", "1-50") # Harmonics up to 50th orders - # Check CSV results for iron loss (You cannot check this for Freq study) # CSV and save space - study.GetStudyProperties().SetValue( - "CsvOutputPath", dir_csv_output_folder - ) # it's folder rather than file! - study.GetStudyProperties().SetValue("CsvResultTypes", self.config.csv_results) - study.GetStudyProperties().SetValue( - "DeleteResultFiles", self.config.del_results_after_calc - ) - study.GetMaterial("Shaft").SetValue("OutputResult", 0) - study.GetMaterial("Coils").SetValue("OutputResult", 0) - - - # Rotor - if True: - cond = study.CreateCondition("Ironloss", "IronLossConRotor") - cond.SetValue("RevolutionSpeed", "freq*60/%d" % self.machine_variant.p) - cond.ClearParts() - sel = cond.GetSelection() - sel.SelectPartByPosition(self.machine_variant.r_sh + 0.1 * self.machine_variant.d_r1, 1, 0) - - cond.AddSelected(sel) - # Use FFT for hysteresis to be consistent with JMAG's results - cond.SetValue("HysteresisLossCalcType", 1) - cond.SetValue("PresetType", 3) - # Specify the reference steps yourself because you don't really know what JMAG is doing behind you - cond.SetValue( - "StartReferenceStep", - number_of_total_steps + 1 - no_of_steps / no_of_rev * 0.25, - ) # 1/4 period in no of steps = no_of_steps / no_of_rev * 0.25 - cond.SetValue("EndReferenceStep", number_of_total_steps) - cond.SetValue("UseStartReferenceStep", 1) - cond.SetValue("UseEndReferenceStep", 1) - cond.SetValue( - "Cyclicity", 4 - ) # specify reference steps for 1/4 period and extend it to whole period - cond.SetValue("UseFrequencyOrder", 1) - cond.SetValue("FrequencyOrder", "1-50") # Harmonics up to 50th orders - self.study_name = study_name - - return study - - - def add_materials(self, study): - # if 'M19' in self.machine_variant.stator_iron_mat["core_material"]: - # study.SetMaterialByName(self.comp_stator_core.name, "M-19 Steel Gauge-29") - # study.GetMaterial(self.comp_stator_core.name).SetValue("Laminated", 1) - # study.GetMaterial(self.comp_stator_core.name).SetValue("LaminationFactor", - # self.machine_variant.stator_iron_mat["core_stacking_factor"]) - - # study.SetMaterialByName(self.comp_rotor_core.name, "M-19 Steel Gauge-29") - # study.GetMaterial(self.comp_rotor_core.name).SetValue("Laminated", 1) - # study.GetMaterial(self.comp_rotor_core.name).SetValue("LaminationFactor", - # self.machine_variant.rotor_iron_mat["core_stacking_factor"]) - - study.SetMaterialByName(self.comp_stator_core.name, - self.machine_variant.stator_iron_mat["core_material"]) - study.GetMaterial(self.comp_stator_core.name).SetValue("Laminated", 1) - study.GetMaterial(self.comp_stator_core.name).SetValue("LaminationFactor", - self.machine_variant.stator_iron_mat["core_stacking_factor"]) - - study.SetMaterialByName(self.comp_rotor_core.name, - self.machine_variant.rotor_iron_mat["core_material"]) - study.GetMaterial(self.comp_rotor_core.name).SetValue("Laminated", 1) - study.GetMaterial(self.comp_rotor_core.name).SetValue("LaminationFactor", - self.machine_variant.rotor_iron_mat["core_stacking_factor"]) - - study.SetMaterialByName(self.comp_shaft.name, - self.machine_variant.shaft_mat["shaft_material"]) - study.GetMaterial(self.comp_shaft.name).SetValue("Laminated", 0) - study.GetMaterial(self.comp_shaft.name).SetValue("EddyCurrentCalculation", 1) - - study.SetMaterialByName(self.comp_winding_layer1.name, - self.machine_variant.coil_mat["coil_material"]) - study.GetMaterial(self.comp_shaft.name).SetValue("UserConductivityType", 1) - - - def add_circuit(self, app, model, study, bool_3PhaseCurrentSource=False): - def add_mp_circuit(study, turns, Rs, x=10, y=10): - # Placing coils/phase windings - coil_name = [] - for i in range(0, 3): - coil_name.append("coil_" + ['U', 'V', 'W'][i]) - study.GetCircuit().CreateComponent("Coil", - coil_name[i]) - study.GetCircuit().CreateInstance(coil_name[i], - x + 4 * i, y) - study.GetCircuit().GetComponent(coil_name[i]).SetValue("Turn", turns) - study.GetCircuit().GetComponent(coil_name[i]).SetValue("Resistance", Rs) - study.GetCircuit().GetInstance(coil_name[i], 0).RotateTo(90) - - self.coil_name = coil_name - - # Connecting all phase windings to a neutral point - for i in range(0, 2): - study.GetCircuit().CreateWire(x + 4 * i, y - 2, x + 4 * (i + 1), y - 2) - - study.GetCircuit().CreateComponent("Ground", "Ground") - study.GetCircuit().CreateInstance("Ground", x + 8, y - 4) - - # Placing current sources - cs_name = [] - for i in range(0, 1): - cs_name.append("cs_" + ['U', 'V', 'W'][i]) - study.GetCircuit().CreateComponent("CurrentSource", cs_name[i]) - study.GetCircuit().CreateInstance(cs_name[i], x + 4 * i, y + 4) - study.GetCircuit().GetInstance(cs_name[i], 0).RotateTo(90) - - self.cs_name = cs_name - - # Terminal Voltage/Circuit Voltage: Check for outputting CSV results - terminal_name = [] - for i in range(0, 1): - terminal_name.append("vp_" + ['U', 'V', 'W'][i]) - study.GetCircuit().CreateTerminalLabel(terminal_name[i], x + 4 * i, y + 2) - study.GetCircuit().CreateComponent("VoltageProbe", terminal_name[i]) - study.GetCircuit().CreateInstance(terminal_name[i], x + 2 + 4 * i, y + 2) - study.GetCircuit().GetInstance(terminal_name[i], 0).RotateTo(90) - - self.terminal_name = terminal_name - - app.ShowCircuitGrid(True) - study.CreateCircuit() - - add_mp_circuit(study, self.machine_variant.Z_q, Rs=self.R_wdg) - - for phase_name in ['U', 'V', 'W']: - study.CreateCondition("FEMCoil", phase_name) - # link between FEM Coil Condition and Circuit FEM Coil - condition = study.GetCondition(phase_name) - condition.SetLink("coil_%s" % (phase_name)) - condition.GetSubCondition("untitled").SetName("delete") - - count = 0 # count indicates which slot the current rightlayer is in. - index = 0 - dict_dir = {"+": 1, "-": 0} - coil_pitch = self.machine_variant.pitch # self.dict_coil_connection[0] - # select the part (via `Set') to assign the FEM Coil condition - for UVW, UpDown in zip( - self.machine_variant.layer_phases[0], self.machine_variant.layer_polarity[0] - ): - - count += 1 - condition = study.GetCondition(UVW) - - # right layer - condition.CreateSubCondition("FEMCoilData", "Coil Set Right %d" % count) - subcondition = condition.GetSubCondition("Coil Set Right %d" % count) - subcondition.ClearParts() - subcondition.AddSet( - model.GetSetList().GetSet( - "coil_%s%s%s %d" % ("right_", UVW, UpDown, count) - ), - 0, - ) # right layer - subcondition.SetValue("Direction2D", dict_dir[UpDown]) - - # left layer - if coil_pitch > 0: - if count + coil_pitch <= self.machine_variant.Q: - count_leftlayer = count + coil_pitch - index_leftlayer = index + coil_pitch - else: - count_leftlayer = int(count + coil_pitch - self.machine_variant.Q) - index_leftlayer = int(index + coil_pitch - self.machine_variant.Q) - else: - if count + coil_pitch > 0: - count_leftlayer = count + coil_pitch - index_leftlayer = index + coil_pitch - else: - count_leftlayer = int(count + coil_pitch + self.machine_variant.Q) - index_leftlayer = int(index + coil_pitch + self.machine_variant.Q) - - # Check if it is a distributed windg??? - if self.machine_variant.pitch == 1: - UVW = self.machine_variant.layer_phases[1][index_leftlayer] - UpDown = self.machine_variant.layer_polarity[1][index_leftlayer] - else: - if self.machine_variant.layer_phases[1][index_leftlayer] != UVW: - print("[Warn] Potential bug in your winding layout detected.") - raise Exception("Bug in winding layout detected.") - if UpDown == "+": - UpDown = "-" - else: - UpDown = "+" - - condition.CreateSubCondition( - "FEMCoilData", "Coil Set Left %d" % count_leftlayer - ) - subcondition = condition.GetSubCondition( - "Coil Set Left %d" % count_leftlayer - ) - subcondition.ClearParts() - subcondition.AddSet( - model.GetSetList().GetSet( - "coil_%s%s%s %d" % ("left_", UVW, UpDown, count_leftlayer) - ), - 0, - ) # left layer - subcondition.SetValue("Direction2D", dict_dir[UpDown]) - index += 1 - # clean up - for phase_name in ['U', 'V', 'W']: - condition = study.GetCondition(phase_name) - condition.RemoveSubCondition("delete") - - - def set_currents_sequence(self, I, freq, phi_0, app, study): - # Setting current values after creating a circuit using "add_mp_circuit" method - # "freq" variable cannot be used here. So pay extra attention when you - # create new case of a different freq. - for i in range(0, 1): - func = app.FunctionFactory().Composite() - f1 = app.FunctionFactory().Sin(I, freq, - - 360 / 3 * i + phi_0 + 90) - func.AddFunction(f1) - study.GetCircuit().GetComponent(self.cs_name[i]).SetFunction(func) - - - def add_time_step_settings(self, time_interval, no_of_steps, app, study): - - DM = app.GetDataManager() - DM.CreatePointArray("point_array/timevsdivision", "SectionStepTable") - refarray = [[0 for i in range(3)] for j in range(3)] - refarray[0][0] = 0 - refarray[0][1] = 1 - refarray[0][2] = 50 - refarray[1][0] = time_interval - refarray[1][1] = no_of_steps - refarray[1][2] = 50 - DM.GetDataSet("SectionStepTable").SetTable(refarray) - number_of_total_steps = ( - 1 + no_of_steps - ) # don't forget to modify here! - study.GetStep().SetValue("Step", number_of_total_steps) - study.GetStep().SetValue("StepType", 3) - study.GetStep().SetTableProperty("Division", DM.GetDataSet("SectionStepTable")) - - - def mesh_study(self, app, model, study): - - # mesh - print("------------------Adding mesh") - self.add_mesh(study, model) - - # Export Image - app.View().ShowAllAirRegions() - app.View().ShowMesh() - app.View().Zoom(3) - app.View().Pan(-self.machine_variant.r_si / 1000, 0) - app.ExportImageWithSize( - self.design_results_folder + self.project_name + "mesh.png", 2000, 2000 - ) - app.View().ShowModel() - - - def add_mesh(self, study, model): - # this is for multi slide planes, which we will not be usin - refarray = [[0 for i in range(2)] for j in range(1)] - refarray[0][0] = 3 - refarray[0][1] = 1 - study.GetMeshControl().GetTable("SlideTable2D").SetTable(refarray) - - study.GetMeshControl().SetValue("MeshType", 1) # make sure this has been exe'd: - study.GetMeshControl().SetValue( - "RadialDivision", self.config.airgap_mesh_radial_div - ) # for air region near which motion occurs - study.GetMeshControl().SetValue( - "CircumferentialDivision", self.config.airgap_mesh_circum_div - ) # 1440) # for air region near which motion occurs - study.GetMeshControl().SetValue( - "AirRegionScale", self.config.mesh_air_region_scale - ) # [Model Length]: Specify a value within (1.05 <= value < 1000) - study.GetMeshControl().SetValue("MeshSize", self.config.mesh_size) - study.GetMeshControl().SetValue("AutoAirMeshSize", 0) - study.GetMeshControl().SetValue( - "AirMeshSize", self.config.mesh_size_rotor - ) # mm - study.GetMeshControl().SetValue("Adaptive", 0) - - study.GetMeshControl().CreateCondition("RotationPeriodicMeshAutomatic", - "autoRotMesh") # with this you can choose to set CircumferentialDivision automatically - - study.GetMeshControl().CreateCondition("Part", "ShaftMeshCtrl") - study.GetMeshControl().GetCondition("ShaftMeshCtrl").SetValue("Size", 1) # 10 mm - study.GetMeshControl().GetCondition("ShaftMeshCtrl").ClearParts() - study.GetMeshControl().GetCondition("ShaftMeshCtrl").AddSet(model.GetSetList().GetSet("ShaftSet"), 0) - - def mesh_all_cases(study): - numCase = study.GetDesignTable().NumCases() - for case in range(0, numCase): - study.SetCurrentCase(case) - if not study.HasMesh(): - study.CreateMesh() - - mesh_all_cases(study) - - - def run_study(self, app, study, toc): - if not self.config.jmag_scheduler: - print("-----------------------Running JMAG...") - study.RunAllCases() - msg = "Time spent on %s is %g s." % (study.GetName(), clock_time() - toc) - print(msg) - else: - print("Submit to JMAG_Scheduler...") - job = study.CreateJob() - job.SetValue("Title", study.GetName()) - job.SetValue("Queued", True) - job.Submit(False) # False:CurrentCase, True:AllCases - # wait and check - - app.Save() - - - def prepare_section( - self, list_regions, tool, bMirrorMerge=True, bRotateMerge=True - ): # csToken is a list of cross section's token - - def regionCircularPattern360Origin(region, tool, bMerge=True): - # index is used to define name of region - - Q_float = float(tool.iRotateCopy) # don't ask me, ask JSOL - circular_pattern = tool.sketch.CreateRegionCircularPattern() - circular_pattern.SetProperty("Merge", bMerge) - - ref2 = tool.doc.CreateReferenceFromItem(region) - circular_pattern.SetPropertyByReference("Region", ref2) - face_region_string = circular_pattern.GetProperty("Region") - - circular_pattern.SetProperty("CenterType", 2) # origin I guess - - circular_pattern.SetProperty("Angle", "360/%d" % Q_float) - circular_pattern.SetProperty("Instance", str(Q_float)) - - list_region_objects = [] - for idx, list_segments in enumerate(list_regions): - # Region - tool.doc.GetSelection().Clear() - for segment in list_segments: - tool.doc.GetSelection().Add(tool.sketch.GetItem(segment.draw_token.GetName())) - - tool.sketch.CreateRegions() - - if idx == 0: - region_object = tool.sketch.GetItem( - "Region" - ) # This is how you get access to the region you create. - else: - region_object = tool.sketch.GetItem( - "Region.%d" % (idx + 1) - ) # This is how you get access to the region you create. - list_region_objects.append(region_object) - # raise - - for idx, region_object in enumerate(list_region_objects): - # Mirror - if tool.bMirror == True: - if tool.edge4Ref is None: - tool.regionMirrorCopy( - region_object, - edge4Ref=None, - symmetryType=2, - bMerge=bMirrorMerge, - ) # symmetryType=2 means x-axis as ref - else: - tool.regionMirrorCopy( - region_object, - edge4Ref=tool.edge4ref, - symmetryType=None, - bMerge=bMirrorMerge, - ) # symmetryType=2 means x-axis as ref - - # RotateCopy - if tool.iRotateCopy != 0: - # print('Copy', self.iRotateCopy) - regionCircularPattern360Origin( - region_object, tool, bMerge=bRotateMerge - ) - - tool.sketch.CloseSketch() - return list_region_objects - - def extract_JMAG_results(self, path, study_name): - fem_coil_flux_path = path + study_name + "_flux_of_fem_coil.csv" - - flux_df = pd.read_csv(fem_coil_flux_path, skiprows=7) - - fea_data = { - "coil_flux_linkages": flux_df, - "current_peak": self.I_hat, - } - - return fea_data \ No newline at end of file diff --git a/mach_eval/analyzers/electromagnetic/SynR/SynR_inductance_config.py b/mach_eval/analyzers/electromagnetic/SynR/SynR_inductance_config.py deleted file mode 100644 index ab9dac7c..00000000 --- a/mach_eval/analyzers/electromagnetic/SynR/SynR_inductance_config.py +++ /dev/null @@ -1,28 +0,0 @@ -class SynR_Inductance_Config: - def __init__(self, **kwargs) -> None: - # attributes for the number of rev and steps - self.no_of_rev = kwargs["no_of_rev"] # number of revolutions - self.no_of_steps = kwargs["no_of_steps"] # number steps - - self.mesh_size = kwargs["mesh_size"] # generic mesh size for overall model [mm] - self.mesh_size_rotor = kwargs["mesh_size_rotor"] # mesh size for rotor [mm] - self.airgap_mesh_radial_div = kwargs["airgap_mesh_radial_div"] # radial divisions of airgap mesh - self.airgap_mesh_circum_div = kwargs["airgap_mesh_circum_div"] # circumferential divisions of airgap mesh - self.mesh_air_region_scale = kwargs["mesh_air_region_scale"] # factor by which air region is scaled beyond machine model - - # results and paths - self.only_table_results = kwargs["only_table_results"] # if True: no mesh or field results are extracted - self.csv_results = kwargs["csv_results"] # data to be extracted from JMAG in csv format - self.del_results_after_calc = kwargs["del_results_after_calc"] # Flag to delete result plot files after calculation - self.run_folder = kwargs["run_folder"] # folder in which JMAG files will reside - self.jmag_csv_folder = kwargs["jmag_csv_folder"] # folder in which csv files from solve are extracted to - - # other settings - # set maximum iteration number of nonlinear calculation required until solution converges - self.max_nonlinear_iterations = kwargs["max_nonlinear_iterations"] - self.multiple_cpus = kwargs["multiple_cpus"] # True if multiple cores are required - self.num_cpus = kwargs["num_cpus"] # number of cpus or cores used. Only value if multiple_cpus = True - self.jmag_scheduler = kwargs["jmag_scheduler"] # True if it is desired to schedule jobs instead of solving immediately - self.jmag_visible = kwargs["jmag_visible"] # JMAG application visible if true - self.scale_axial_length = kwargs["scale_axial_length"] # True: scale axial length to get the required rated torque - self.jmag_version = kwargs["jmag_version"] # JMAG application version \ No newline at end of file diff --git a/mach_eval/analyzers/electromagnetic/bspm/jmag_2d.py b/mach_eval/analyzers/electromagnetic/bspm/jmag_2d.py index 490c1545..a9f6eb7b 100644 --- a/mach_eval/analyzers/electromagnetic/bspm/jmag_2d.py +++ b/mach_eval/analyzers/electromagnetic/bspm/jmag_2d.py @@ -305,7 +305,7 @@ def group(name, id_list): id_shaft = part_ID_list[1] partIDRange_Magnet = part_ID_list[2 : int(2 + self.machine_variant.p * 2)] # id_sleeve = part_ID_list[int(2 + self.machine_variant.p * 2)] - id_statorCore = part_ID_list[int(2 + self.machine_variant.p * 2) + 1] + id_statorCore = part_ID_list[int(2 + self.machine_variant.p * 2)] partIDRange_Coil = part_ID_list[ int(1 + self.machine_variant.p * 2) + 2 : int(2 + self.machine_variant.p * 2) @@ -318,6 +318,12 @@ def group(name, id_list): group("Magnet", partIDRange_Magnet) group("Coils", partIDRange_Coil) + # """ Set Parts names """ + + app.GetModel(0).SetPartName(id_backiron, u"NotchedRotor") + app.GetModel(0).SetPartName(id_shaft, u"Shaft") + app.GetModel(0).SetPartName(id_statorCore, u"StatorCore") + """ Add Part to Set for later references """ def add_part_to_set(name, x, y, ID=None): diff --git a/mach_eval/analyzers/electromagnetic/flux_linkage_analyzer.py b/mach_eval/analyzers/electromagnetic/flux_linkage_analyzer.py new file mode 100644 index 00000000..6f950af5 --- /dev/null +++ b/mach_eval/analyzers/electromagnetic/flux_linkage_analyzer.py @@ -0,0 +1,164 @@ +import os +import pandas as pd +import numpy as np +from time import time as clock_time + +class FluxLinkageJMAG_Problem: + """Problem class for processing FEA flux linkages + Attributes: + toolJmag: JMAG tool required for FEA scripting + phase_names: names of phases in machine + rated_current: rated current of machine + """ + + def __init__(self, toolJmag, phase_names, rated_current): + self.toolJmag = toolJmag + self.phase_names = phase_names + self.rated_current = rated_current + + def add_study_parameters(self, study_name): + + # Pre-processing + self.model.SetName("Machine_FluxLinkage_Project") + + # Initialize JMAG application + self.app.SetCurrentStudy(study_name) + study = self.model.GetStudy(study_name) + + # Set csv folder output + self.results_filepath = os.path.dirname(__file__) + "/run_data/" + + # Create output folder + if not os.path.isdir(self.results_filepath): + os.makedirs(self.results_filepath) + + # Create csv output file loction + study.GetStudyProperties().SetValue("CsvOutputPath", self.results_filepath) + study.GetStudyProperties().SetValue("CsvResultTypes", "FEMCoilFlux") + + self.app.SetCurrentStudy(study_name) + + return study + + + def delete_operating_points(self): + + self.app = self.toolJmag.jd + self.model = self.app.GetCurrentModel() + + self.app.SetCurrentStudy("Machine_FluxLinkage_Study") + study = self.model.GetStudy("Machine_FluxLinkage_Study") + # Creating circuit function and zeroing all values + cs_name = [] + for i in range(0, len(self.phase_names)): + cs_name.append("cs_" + self.phase_names[i]) + f1 = self.app.FunctionFactory().Sin(0, 0.000001, 90) + func = self.app.FunctionFactory().Composite() + func.AddFunction(f1) + study.GetCircuit().GetComponent(cs_name[i]).SetFunction(func) + self.app.GetModel(0).GetStudy(0).GetDesignTable().AddParameterVariableName("cs_%s (CurrentSource): 1 (Composite Function): Amplitude" % self.phase_names[i]) + + def new_operating_point(self): + + # Creating new operating points + self.app.GetModel(0).GetStudy(0).GetDesignTable().AddCase() + op_pt = self.app.GetModel(0).GetStudy(0).GetDesignTable() + + return op_pt + + + def set_phase_currents(self, op_pt, m, i, I_hat): + + # Setting current values of single current source depending on phase excitation + op_pt.SetValue(m+1, i, I_hat) + + return I_hat + + + def run_all_studies(self): + + study = self.add_study_parameters("Machine_FluxLinkage_Study") + print("-----------------------Running JMAG...") + study.RunAllCases() + + + def get_flux_linkages(self, study_name): + + # Read data from output file + linkage_files = {} + linkages = [] + for i in range(len(self.phase_names)+1): + linkage_files[i] = pd.read_csv(self.results_filepath + study_name + "_flux_of_fem_coil.csv", skiprows=4+len(self.phase_names), usecols=[0]+np.arange(len(self.phase_names)*i+1,len(self.phase_names)*i+len(self.phase_names)+1).tolist()) + linkages.append(linkage_files[i]) + + # Find rotor angle array from time data + zero_current_linkage = linkages[0].to_numpy() + time = zero_current_linkage[:,0] + rotor_angle = 360*time/(max(time) - min(time)), + + # Create current vector + current_vector = np.full((1, len(self.phase_names)), self.rated_current, dtype=int) + + # Write output data of analyzer + fea_data = { + "phase_currents": current_vector, + "rotor_angle": rotor_angle, + "linkages": linkages, + "name_of_phases": self.phase_names, + } + + return fea_data + + +class FluxLinkageJMAG_Analyzer: + """Calcuates generates flux linkages from FEA + + Args: + problem: object of type Flux_Linkage_Problem holding flux linkage simulation instructions + Returns: + fea_data: Data dictionary containing information for post-processing + """ + + def analyze(self, problem): + + ################################################################ + # 02. Create all operating points! One per phase + 0 current case + ################################################################ + + # Zero all Currents + problem.delete_operating_points() + + currents = [] + op_pt = {} + m = 0 + # Create cases for each operating point + #for m in range(len(problem.phase_names)+1): + # op_pt[m] = problem.new_operating_point(m) + # problem.set_phase_currents(op_pt[m], problem.rated_current, m) + while m < len(problem.phase_names): + op_pt[m] = problem.new_operating_point() + i = 0 + currents.append([]) + while i < len(problem.phase_names): + if i == m: + currents[m].append([i]) + currents[m][i] = problem.set_phase_currents(op_pt[m], m, i, problem.rated_current) + else: + currents[m].append([i]) + currents[m][i] = problem.set_phase_currents(op_pt[m], m, i, 0) + i = i + 1 + m = m + 1 + + ################################################################ + # 03. Run electromagnetic studies + ################################################################ + + problem.run_all_studies() + + #################################################### + # 04. Extract Results + #################################################### + + flux_linkages = problem.get_flux_linkages("Machine_FluxLinkage_Study") + + return flux_linkages, currents \ No newline at end of file diff --git a/mach_eval/analyzers/electromagnetic/inductance_analyzer.py b/mach_eval/analyzers/electromagnetic/inductance_analyzer.py new file mode 100644 index 00000000..8fa54421 --- /dev/null +++ b/mach_eval/analyzers/electromagnetic/inductance_analyzer.py @@ -0,0 +1,93 @@ +import numpy as np +import pandas as pd +import numpy as np + +class InductanceProblem: + """Problem class for inductance data processing + Attributes: + I_hat: peak current value for inductance calculations + linkages: arrays of flux linkages, should be (m+1) x m for m-phase machine + rotor_angle: array of rotor angle values for inductance calculations + name_of_phases: array of name of m-phases for machine + """ + + def __init__(self, I_hat, linkages, rotor_angle, name_of_phases): + self.I_hat = I_hat + self.linkages = linkages + self.rotor_angle = rotor_angle + self.name_of_phases = name_of_phases + + +class InductanceAnalyzer: + def __init__(self, clarke_trans_matrix): + self.clarke_trans_matrix = clarke_trans_matrix + + def analyze(self, problem): + """Calcuates abc, alpha-beta, and dq inductances + + Args: + problem: object of type Inductance_Problem holding torque data + Returns: + data: dictionary of inductances and rotor angles + """ + + L_abc = [] + flux_linkages = [] # make mxm list + flux_linkages_files = {} + L_zero = [] + flux_linkages_zero = [] # make 1xm list + flux_linkages_files_zero = {} + + for col in range(len(problem.name_of_phases)+1): + if col == 0: + flux_linkages_files_zero[col] = problem.linkages[0] + flux_linkages_zero.append([]) + L_zero.append([]) + for row in range(len(problem.name_of_phases)): + flux_linkages_zero[col].append(row) + L_zero[col].append(row) + flux_linkages_zero[col][row] = np.array(flux_linkages_files_zero[col]["coil_%s" % problem.name_of_phases[row]]) + else: + flux_linkages_files[col-1] = problem.linkages[col] + flux_linkages.append([]) + L_abc.append([]) + for row in range(len(problem.name_of_phases)): + flux_linkages[col-1].append(row) + L_abc[col-1].append(row) + flux_linkages[col-1][row] = np.array(flux_linkages_files[col-1]["coil_%s.%s" % (problem.name_of_phases[row], str(col))]) + L_abc[col-1][row] = (flux_linkages[col-1][row] - flux_linkages_zero[0][row])/problem.I_hat + + L_abc = np.array(L_abc) + + L_alpha_beta = [] + for i in range(len(L_abc[0][0])): + L_alpha_beta.append([]) + L_alpha_beta[i] = np.dot(self.clarke_trans_matrix,L_abc[:,:,i]) + L_alpha_beta[i] = np.dot(L_alpha_beta[i],np.linalg.inv(self.clarke_trans_matrix)) + + L_alpha_beta = np.array(L_alpha_beta) + + L_dq = [] + for i in range(len(L_alpha_beta)): + L_dq.append([]) + theta = -problem.rotor_angle[i]*np.pi/180 + park_trans_matrix = np.array([[np.cos(theta), -np.sin(theta), 0], [np.sin(theta), np.cos(theta), 0], [0, 0, 1]]) + L_dq[i] = np.dot(park_trans_matrix,L_alpha_beta[i,:,:]) + L_dq[i] = np.dot(L_dq[i],np.linalg.inv(park_trans_matrix)) + + L_dq = np.array(L_dq) + + data = self.extract_results(problem, L_abc, L_alpha_beta, L_dq) + + return data + + def extract_results(self, problem, L_abc, L_alpha_beta, L_dq): + + data = { + "rotor_angle": problem.rotor_angle, + "Labc": L_abc, + "Lalphabeta": L_alpha_beta, + "Ldq": L_dq, + } + + return data \ No newline at end of file diff --git a/mach_eval/machines/machine.py b/mach_eval/machines/machine.py index d53deb0d..48bf1a04 100644 --- a/mach_eval/machines/machine.py +++ b/mach_eval/machines/machine.py @@ -62,6 +62,10 @@ def required_materials(): @property def no_of_layers(self): return self._winding_dict["no_of_layers"] + + @property + def name_of_phases(self): + return self._winding_dict["name_of_phases"] @property def layer_phases(self):