diff --git a/machine_design/design.py b/machine_design/design.py index 0489031..3324a50 100644 --- a/machine_design/design.py +++ b/machine_design/design.py @@ -27,10 +27,24 @@ def load(cls, file_name: str, **kwargs) -> "Design": return cls(m2d) def set_parameters(self) -> None: - # materials + self.setup_name = "Setup1" + self.set_iron() + self.set_geom_params() + self.set_slot_params() + self.set_winds_params() + self.set_mod_params() + self.set_oper_params() + self.set_rot_points() + self.set_derived_params() + self.set_solution_expressions() + self.set_udp_par_list_stator() + self.set_output_vars() + self.set_post_params() + + def set_iron(self): self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" - # main definitions + def set_geom_params(self): self.geom_params = { "DiaStatorGap": "79mm", "DiaStatorYoke": "125mm", @@ -40,7 +54,8 @@ def set_parameters(self) -> None: "DiaShaft": "25mm", "StackLength": "85mm", } - # stator slot + + def set_slot_params(self): self.slot_params = { "Hs0": "0.95mm", "Hs1": "0.31mm", @@ -51,7 +66,8 @@ def set_parameters(self) -> None: "Rs": "1.5mm", "SetAngle": "10deg", } - # winding + + def set_winds_params(self): self.wind_params = { "Layers": "1", "ParallelPaths": "1", @@ -60,17 +76,19 @@ def set_parameters(self) -> None: "SpaceLayers": "0.2mm", "Nc": "68", # turns per coil } - # model parameters - PolePairs = 2 - f = 50 # [Hz] - RotSpeed = 60 * f / PolePairs # [rpm] + + def set_mod_params(self): + self.PolePairs = 2 self.mod_params = { - "Poles": f"2*{PolePairs}", + "Poles": f"2*{self.PolePairs}", "ModelLength": "85mm", "SymmetryFactor": "Poles", "StatorSkewAngle": "0deg", } - # operation parameters + + def set_oper_params(self): + f = 50 # [Hz] + RotSpeed = 60 * f / self.PolePairs # [rpm] self.oper_params = { "Im": "1.5*sqrt(2)A", "epsI": "pi/4", # current angle @@ -80,26 +98,89 @@ def set_parameters(self) -> None: "Nper": "1/6", # number of included periods "PointPer": "101", # number of time points per period } + + def set_rot_points(self): self.rot_points = [ ["DiaShaft/2*cos(360deg/SymmetryFactor)", "DiaShaft/2*sin(360deg/SymmetryFactor)", "0mm"], ["DiaShaft/2*cos(360deg/(2*SymmetryFactor))", "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", "0mm"], ["DiaShaft/2", "0mm", "0mm"], ["DiaStatorGap/2-Airgap", "0mm", "0mm"], - [ - "(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", - "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", - "0mm", - ], - [ - "(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", - "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", - "0mm", - ], + ["(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", "0mm"], + ["(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", "0mm"], ] - self.setup_name = "Setup1" + + def set_derived_params(self): self.rotor_r_min = self.mm_to_str("geom_params", "DiaShaft") / 2 self.rotor_r_max = self.mm_to_str("geom_params", "DiaStatorGap") / 2 - self.mm_to_str("geom_params", "Airgap") + def set_solution_expressions(self): + self.solution_expressions = "Moving1.Torque" + + def set_udp_par_list_stator(self): + self.udp_par_list_stator = [ + ["DiaGap", "DiaStatorGap"], + ["DiaYoke", "DiaStatorYoke"], + ["Length", "0mm"], + ["Skew", "0deg"], + ["Slots", "SlotNumber"], + ["SlotType", "SlotType"], + ["Hs0", "Hs0"], + ["Hs01", "0mm"], + ["Hs1", "Hs1"], + ["Hs2", "Hs2"], + ["Bs0", "Bs0"], + ["Bs1", "Bs1"], + ["Bs2", "Bs2"], + ["Rs", "Rs"], + ["FilletType", "0"], + ["HalfSlot", "0"], + ["SegAngle", "0deg"], + ["LenRegion", "0mm"], + ["InfoCore", "0"], + ] + + def set_output_vars(self): + self.output_vars = { + "pos": "(Moving1.Position -InitPos) * Poles/2", + "cos0": "cos(pos)", + "cos1": "cos(pos-2*PI/3)", + "cos2": "cos(pos-4*PI/3)", + "sin0": "sin(pos)", + "sin1": "sin(pos-2*PI/3)", + "sin2": "sin(pos-4*PI/3)", + "Lad": "L(PhaseA,PhaseA)*cos0 + L(PhaseA,PhaseB)*cos1 + L(PhaseA,PhaseC)*cos2", + "Laq": "L(PhaseA,PhaseA)*sin0 + L(PhaseA,PhaseB)*sin1 + L(PhaseA,PhaseC)*sin2", + "Lbd": "L(PhaseB,PhaseA)*cos0 + L(PhaseB,PhaseB)*cos1 + L(PhaseB,PhaseC)*cos2", + "Lbq": "L(PhaseB,PhaseA)*sin0 + L(PhaseB,PhaseB)*sin1 + L(PhaseB,PhaseC)*sin2", + "Lcd": "L(PhaseC,PhaseA)*cos0 + L(PhaseC,PhaseB)*cos1 + L(PhaseC,PhaseC)*cos2", + "Lcq": "L(PhaseC,PhaseA)*sin0 + L(PhaseC,PhaseB)*sin1 + L(PhaseC,PhaseC)*sin2", + "L_d": "(Lad*cos0 + Lbd*cos1 + Lcd*cos2) * 2/3", + "L_q": "(Laq*sin0 + Lbq*sin1 + Lcq*sin2) * 2/3", + "Flux_d": "(FluxLinkage(PhaseA)*cos0+FluxLinkage(PhaseB)*cos1+FluxLinkage(PhaseC)*cos2)*2/3", + "Flux_q": "-(FluxLinkage(PhaseA)*sin0+FluxLinkage(PhaseB)*sin1+FluxLinkage(PhaseC)*sin2)*2/3", + "Ui_d": "(InducedVoltage(PhaseA)*cos0+InducedVoltage(PhaseB)*cos1+InducedVoltage(PhaseC)*cos2)*2/3", + "Ui_q": "-(InducedVoltage(PhaseA)*sin0+InducedVoltage(PhaseB)*sin1+InducedVoltage(PhaseC)*sin2)*2/3", + "I_d": "(InputCurrent(PhaseA)*cos0 + InputCurrent(PhaseB)*cos1 + InputCurrent(PhaseC)*cos2)*2/3", + "I_q": "-(InputCurrent(PhaseA)*sin0 + InputCurrent(PhaseB)*sin1 + InputCurrent(PhaseC)*sin2)*2/3", + "Irms": "sqrt(I_d^2+I_q^2)/sqrt(2)", + } + + def set_post_params(self): + self.post_params = { # reports + ("InducedVoltage(PhaseA)", "InducedVoltage(PhaseB)", "InducedVoltage(PhaseC)"): "InducedVoltage", + ("Moving1.Torque"): "Torque", + ("InputCurrent(PhaseA)", "InputCurrent(PhaseB)", "InputCurrent(PhaseC)"): "Current", + ( + "FluxLinkage(PhaseA)", + "FluxLinkage(PhaseB)", + "FluxLinkage(PhaseC)", + ): "FluxLinkage", + ("I_d", "I_q"): "Current_dq", + ("Flux_d", "Flux_q"): "FluxLinkage_dq", + ("Ui_d", "Ui_q"): "InducedVoltage_dq", + ("L_d", "L_q"): "Inductance_dq", + } + def create_stator(self) -> None: m2d = self.m2d modeler = m2d.modeler @@ -168,30 +249,9 @@ def create_stator(self) -> None: modeler.fit_all() # Stator geometry - udp_par_list_stator = [ - ["DiaGap", "DiaStatorGap"], - ["DiaYoke", "DiaStatorYoke"], - ["Length", "0mm"], - ["Skew", "0deg"], - ["Slots", "SlotNumber"], - ["SlotType", "SlotType"], - ["Hs0", "Hs0"], - ["Hs01", "0mm"], - ["Hs1", "Hs1"], - ["Hs2", "Hs2"], - ["Bs0", "Bs0"], - ["Bs1", "Bs1"], - ["Bs2", "Bs2"], - ["Rs", "Rs"], - ["FilletType", "0"], - ["HalfSlot", "0"], - ["SegAngle", "0deg"], - ["LenRegion", "0mm"], - ["InfoCore", "0"], - ] stator_id = modeler.create_udp( dll="RMxprt/SlotCore.dll", - parameters=udp_par_list_stator, + parameters=self.udp_par_list_stator, library="syslib", name="Stator", # SolveInside="True", @@ -211,7 +271,7 @@ def create_stator(self) -> None: coil_id.color = (255, 128, 0) coil_id.transparency = 0.0 modeler.rotate(assignment=coil_id, axis="Z", angle="360deg/SlotNumber/2") - coil_id.duplicate_around_axis(axis="Z", angle="360deg/SlotNumber", clones="CoilPitch", create_new_objects=True) + coil_id.duplicate_around_axis(axis="Z", angle="360deg/SlotNumber", clones="SlotNumber/Poles", create_new_objects=True) id_coils = modeler.get_objects_w_string(string_name="Coil", case_sensitive=True) # Create section of machine @@ -260,6 +320,68 @@ def create_stator(self) -> None: ) m2d.assign_vector_potential(assignment=id_bc_az, vector_value=0, boundary="A0") + self.assign_stator_coils() + + # Mesh operation + m2d.mesh.assign_length_mesh( + assignment=id_coils, + inside_selection=True, + maximum_length=3, + maximum_elements=None, + name="coils", + ) + m2d.mesh.assign_length_mesh( + assignment=stator_id, + inside_selection=True, + maximum_length=3, + maximum_elements=None, + name="stator", + ) + + # core loss + m2d.set_core_losses("Stator", core_loss_on_field=False) + + # inductance calculation + self.inductance_computation() + + # model depth + m2d.model_depth = "StackLength" + # symmetry + m2d.change_symmetry_multiplier("SymmetryFactor") + # Calculation setup + setup = m2d.create_setup(name=self.setup_name) + setup.props["StopTime"] = "Nper/f" + setup.props["TimeStep"] = "1/(f*(PointPer-1))" + setup.props["SaveFieldsType"] = "None" + setup.props["OutputPerObjectCoreLoss"] = False + setup.props["OutputPerObjectSolidLoss"] = True + setup.props["OutputError"] = True + setup.update() + m2d.validate_simple() + + for k, v in self.output_vars.items(): + m2d.create_output_variable(k, v) + + for k, v in self.post_params.items(): + expressions = list(k) if isinstance(k, tuple) else [k] # if multiple report, use list(k). Else, use k + m2d.post.create_report( + expressions=expressions, + setup_sweep_name="", + domain="Sweep", + variations=None, + primary_sweep_variable="Time", + secondary_sweep_variable=None, + report_category=None, + plot_type="Rectangular Plot", + context=None, + subdesign_id=None, + polyline_points=1001, + plot_name=v, + ) + + def assign_stator_coils(self): + m2d = self.m2d + # Excitations I_A = "Im * cos(2*pi*f*time+epsI)" I_B = "Im * cos(2*pi*f*time-120deg+epsI)" @@ -348,110 +470,14 @@ def create_stator(self) -> None: ) m2d.add_winding_coils(assignment="PhaseC", coils=["CS4", "CS5", "CS6"]) - # Mesh operation - m2d.mesh.assign_length_mesh( - assignment=id_coils, - inside_selection=True, - maximum_length=3, - maximum_elements=None, - name="coils", - ) - m2d.mesh.assign_length_mesh( - assignment=stator_id, - inside_selection=True, - maximum_length=3, - maximum_elements=None, - name="stator", - ) - - # core loss - m2d.set_core_losses("Stator", core_loss_on_field=False) - - # inductance calculation - m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=False) - # model depth - m2d.model_depth = "StackLength" - # symmetry - m2d.change_symmetry_multiplier("SymmetryFactor") - # Calculation setup - setup = m2d.create_setup(name=self.setup_name) - setup.props["StopTime"] = "Nper/f" - setup.props["TimeStep"] = "1/(f*(PointPer-1))" - setup.props["SaveFieldsType"] = "None" - setup.props["OutputPerObjectCoreLoss"] = False - setup.props["OutputPerObjectSolidLoss"] = True - setup.props["OutputError"] = True - setup.update() - m2d.validate_simple() - - # ooutput variables - output_vars = { - "pos": "(Moving1.Position -InitPos) * Poles/2", - "cos0": "cos(pos)", - "cos1": "cos(pos-2*PI/3)", - "cos2": "cos(pos-4*PI/3)", - "sin0": "sin(pos)", - "sin1": "sin(pos-2*PI/3)", - "sin2": "sin(pos-4*PI/3)", - "Lad": "L(PhaseA,PhaseA)*cos0 + L(PhaseA,PhaseB)*cos1 + L(PhaseA,PhaseC)*cos2", - "Laq": "L(PhaseA,PhaseA)*sin0 + L(PhaseA,PhaseB)*sin1 + L(PhaseA,PhaseC)*sin2", - "Lbd": "L(PhaseB,PhaseA)*cos0 + L(PhaseB,PhaseB)*cos1 + L(PhaseB,PhaseC)*cos2", - "Lbq": "L(PhaseB,PhaseA)*sin0 + L(PhaseB,PhaseB)*sin1 + L(PhaseB,PhaseC)*sin2", - "Lcd": "L(PhaseC,PhaseA)*cos0 + L(PhaseC,PhaseB)*cos1 + L(PhaseC,PhaseC)*cos2", - "Lcq": "L(PhaseC,PhaseA)*sin0 + L(PhaseC,PhaseB)*sin1 + L(PhaseC,PhaseC)*sin2", - "L_d": "(Lad*cos0 + Lbd*cos1 + Lcd*cos2) * 2/3", - "L_q": "(Laq*sin0 + Lbq*sin1 + Lcq*sin2) * 2/3", - "Flux_d": "(FluxLinkage(PhaseA)*cos0+FluxLinkage(PhaseB)*cos1+FluxLinkage(PhaseC)*cos2)*2/3", - "Flux_q": "-(FluxLinkage(PhaseA)*sin0+FluxLinkage(PhaseB)*sin1+FluxLinkage(PhaseC)*sin2)*2/3", - "Ui_d": "(InducedVoltage(PhaseA)*cos0+InducedVoltage(PhaseB)*cos1+InducedVoltage(PhaseC)*cos2)*2/3", - "Ui_q": "-(InducedVoltage(PhaseA)*sin0+InducedVoltage(PhaseB)*sin1+InducedVoltage(PhaseC)*sin2)*2/3", - "I_d": "(InputCurrent(PhaseA)*cos0 + InputCurrent(PhaseB)*cos1 + InputCurrent(PhaseC)*cos2)*2/3", - "I_q": "-(InputCurrent(PhaseA)*sin0 + InputCurrent(PhaseB)*sin1 + InputCurrent(PhaseC)*sin2)*2/3", - "Irms": "sqrt(I_d^2+I_q^2)/sqrt(2)", - } - for k, v in output_vars.items(): - m2d.create_output_variable(k, v) - - # Definitions for plots - post_params = { # reports - ("InducedVoltage(PhaseA)", "InducedVoltage(PhaseB)", "InducedVoltage(PhaseC)"): "InducedVoltage", - ("Moving1.Torque"): "Torque", - ("InputCurrent(PhaseA)", "InputCurrent(PhaseB)", "InputCurrent(PhaseC)"): "Current", - ( - "FluxLinkage(PhaseA)", - "FluxLinkage(PhaseB)", - "FluxLinkage(PhaseC)", - ): "FluxLinkage", - ("I_d", "I_q"): "Current_dq", - ("Flux_d", "Flux_q"): "FluxLinkage_dq", - ("Ui_d", "Ui_q"): "InducedVoltage_dq", - ("L_d", "L_q"): "Inductance_dq", - } - # Create Report - for k, v in post_params.items(): - expressions = list(k) if isinstance(k, tuple) else [k] # if multiple report, use list(k). Else, use k - m2d.post.create_report( - expressions=expressions, - setup_sweep_name="", - domain="Sweep", - variations=None, - primary_sweep_variable="Time", - secondary_sweep_variable=None, - report_category=None, - plot_type="Rectangular Plot", - context=None, - subdesign_id=None, - polyline_points=1001, - plot_name=v, - ) + def inductance_computation(self): + self.m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=False) def add_rotor(self) -> None: modeler = self.m2d.modeler assert isinstance(modeler, Modeler2D) - rotor_id = modeler.create_polyline( - points=self.rot_points, segment_type=["Arc", "Line", "Arc"], cover_surface=True, name="Rotor" - ) + rotor_id = modeler.create_polyline(points=self.rot_points, segment_type=["Arc", "Line", "Arc"], cover_surface=True, name="Rotor") self.rotor_id = rotor_id rotor_id.material_name = self.Fe rotor_id.color = (192, 192, 192) # rgb @@ -470,15 +496,14 @@ def add_rotor_barrier(self, barrier_points, segment_type=None) -> None: # Convert them into a string format and interpolate points_str = [[str(y) for y in x] for x in barrier_points] - barrier_id = modeler.create_polyline( - points=points_str, segment_type=segment_type, cover_surface=True, name="Barrier" - ) + barrier_id = modeler.create_polyline(points=points_str, segment_type=segment_type, cover_surface=True, name="Barrier") # Remove the barrier self.rotor_id.subtract(barrier_id) modeler.delete(barrier_id) - def compute(self, NUM_CORES: int = 1): + # TODO: change the other arguments to kwargs + def compute(self, *args, NUM_CORES: int = 1): m2d = self.m2d assert m2d.mesh is not None assert m2d.post is not None @@ -493,12 +518,14 @@ def compute(self, NUM_CORES: int = 1): # core loss rotor m2d.set_core_losses("Rotor", core_loss_on_field=False) + self.set_variables(*args) + # Analyze m2d.analyze_setup(self.setup_name, use_auto_settings=False, cores=NUM_CORES) - solutions = m2d.post.get_solution_data(expressions="Moving1.Torque", primary_sweep_variable="Time") + solutions = m2d.post.get_solution_data(expressions=self.solution_expressions, primary_sweep_variable="Time") try: - result = solutions.data_magnitude() + result = self.extract_results(solutions) except AttributeError: result = None @@ -507,6 +534,12 @@ def compute(self, NUM_CORES: int = 1): return result + def set_variables(self, *args): + pass + + def extract_results(self, solutions): + return solutions.data_magnitude() + def delete_rotor(self) -> None: assert isinstance(self.m2d.modeler, Modeler2D) self.m2d.modeler.delete(self.rotor_id) diff --git a/machine_design/design2.py b/machine_design/design2.py new file mode 100644 index 0000000..8b982db --- /dev/null +++ b/machine_design/design2.py @@ -0,0 +1,325 @@ +import numpy as np + +from .design import Design + + +class Design2(Design): + def set_geom_params(self): + super().set_geom_params() + self.geom_params["SlotNumber"] = "40" + + def set_slot_params(self): + super().set_slot_params() + self.slot_params["Bs1"] = "3.0mm" + self.slot_params["Bs2"] = "4.3mm" + self.slot_params["SetAngle"] = "9deg" + + def set_winds_params(self): + super().set_winds_params() + self.wind_params["Nc"] = "113" + + def set_oper_params(self): + f = 50 # [Hz] + RotSpeed = 60 * f / self.PolePairs # [rpm] + self.oper_params = { + "Id1": "0.0A", + "Iq1": "0.0A", + "Id3": "0.0A", + "Iq3": "0.0A", + "epsI1": "atan2(Iq1,Id1)", # current angle, 1st harmonic + "epsI3": "atan2(Iq3,Id3)", # current angle, 1st harmonic + "Im1": "sqrt(Id1^2+Iq1^2)", + "Im3": "sqrt(Id3^2+Iq3^2)", + "InitPos": "-45deg", + "f": f"{f}Hz", + "RotSpeed": f"{RotSpeed}rpm", + "Nper": "1/10", # number of included periods + "PointPer": "101", # number of time points per period + } + + def set_derived_params(self): + pass + + def set_solution_expressions(self): + self.solution_expressions = [ + "V_d1", + "V_q1", + "V_d3", + "V_q3", + "Flux_e_d1", + "Flux_e_q1", + "Flux_e_d3", + "Flux_e_q3", + "I_d1", + "I_q1", + "I_d3", + "I_q3", + "Ld1", + "Ld1q1", + "Ld1d3", + "Ld1q3", + "Lq1", + "Lq1d3", + "Lq1q3", + "Ld3", + "Ld3q3", + "Lq3", + "Moving1.Torque", + ] + + def set_output_vars(self): + self.output_vars = { + "PolePairs": "2", + "RotSign": "1", + "Rstat": "19", + "Lew": "0", + "theta_el": "RotSign*(Moving1.Position - InitPos) * PolePairs - pi", + "cos0_1": "cos(1*(theta_el - 2*PI*0/5))", + "sin0_1": "sin(-1*(theta_el - 2*PI*0/5))", + "cos1_1": "cos(1*(theta_el - 2*PI*1/5))", + "sin1_1": "sin(-1*(theta_el - 2*PI*1/5))", + "cos2_1": "cos(1*(theta_el - 2*PI*2/5))", + "sin2_1": "sin(-1*(theta_el - 2*PI*2/5))", + "cos3_1": "cos(1*(theta_el - 2*PI*3/5))", + "sin3_1": "sin(-1*(theta_el - 2*PI*3/5))", + "cos4_1": "cos(1*(theta_el - 2*PI*4/5))", + "sin4_1": "sin(-1*(theta_el - 2*PI*4/5))", + "cos0_3": "cos(3*(theta_el - 2*PI*0/5))", + "sin0_3": "sin(-3*(theta_el - 2*PI*0/5))", + "cos1_3": "cos(3*(theta_el - 2*PI*1/5))", + "sin1_3": "sin(-3*(theta_el - 2*PI*1/5))", + "cos2_3": "cos(3*(theta_el - 2*PI*2/5))", + "sin2_3": "sin(-3*(theta_el - 2*PI*2/5))", + "cos3_3": "cos(3*(theta_el - 2*PI*3/5))", + "sin3_3": "sin(-3*(theta_el - 2*PI*3/5))", + "cos4_3": "cos(3*(theta_el - 2*PI*4/5))", + "sin4_3": "sin(-3*(theta_el - 2*PI*4/5))", + "Flux_d1": "(FluxLinkage(PhaseA)*cos0_1 + FluxLinkage(PhaseB)*cos1_1 + FluxLinkage(PhaseC)*cos2_1 + FluxLinkage(PhaseD)*cos3_1 + FluxLinkage(PhaseE)*cos4_1) * 2/5", + "Flux_q1": "(FluxLinkage(PhaseA)*sin0_1 + FluxLinkage(PhaseB)*sin1_1 + FluxLinkage(PhaseC)*sin2_1 + FluxLinkage(PhaseD)*sin3_1 + FluxLinkage(PhaseE)*sin4_1) * 2/5", + "Flux_d3": "(FluxLinkage(PhaseA)*cos0_3 + FluxLinkage(PhaseB)*cos1_3 + FluxLinkage(PhaseC)*cos2_3 + FluxLinkage(PhaseD)*cos3_3 + FluxLinkage(PhaseE)*cos4_3) * 2/5", + "Flux_q3": "(FluxLinkage(PhaseA)*sin0_3 + FluxLinkage(PhaseB)*sin1_3 + FluxLinkage(PhaseC)*sin2_3 + FluxLinkage(PhaseD)*sin3_3 + FluxLinkage(PhaseE)*sin4_3) * 2/5", + "Vind_d1": "(InducedVoltage(PhaseA)*cos0_1 + InducedVoltage(PhaseB)*cos1_1 + InducedVoltage(PhaseC)*cos2_1 + InducedVoltage(PhaseD)*cos3_1 + InducedVoltage(PhaseE)*cos4_1) * 2/5", + "Vind_q1": "(InducedVoltage(PhaseA)*sin0_1 + InducedVoltage(PhaseB)*sin1_1 + InducedVoltage(PhaseC)*sin2_1 + InducedVoltage(PhaseD)*sin3_1 + InducedVoltage(PhaseE)*sin4_1) * 2/5", + "Vind_d3": "(InducedVoltage(PhaseA)*cos0_3 + InducedVoltage(PhaseB)*cos1_3 + InducedVoltage(PhaseC)*cos2_3 + InducedVoltage(PhaseD)*cos3_3 + InducedVoltage(PhaseE)*cos4_3) * 2/5", + "Vind_q3": "(InducedVoltage(PhaseA)*sin0_3 + InducedVoltage(PhaseB)*sin1_3 + InducedVoltage(PhaseC)*sin2_3 + InducedVoltage(PhaseD)*sin3_3 + InducedVoltage(PhaseE)*sin4_3) * 2/5", + "V_A": "InducedVoltage(PhaseA) + Rstat*InputCurrent(PhaseA) + Lew*ddt(InputCurrent(PhaseA))", + "V_B": "InducedVoltage(PhaseB) + Rstat*InputCurrent(PhaseB) + Lew*ddt(InputCurrent(PhaseB))", + "V_C": "InducedVoltage(PhaseC) + Rstat*InputCurrent(PhaseC) + Lew*ddt(InputCurrent(PhaseC))", + "V_D": "InducedVoltage(PhaseD) + Rstat*InputCurrent(PhaseD) + Lew*ddt(InputCurrent(PhaseD))", + "V_E": "InducedVoltage(PhaseE) + Rstat*InputCurrent(PhaseE) + Lew*ddt(InputCurrent(PhaseE))", + "V_AC": "V_A - V_C", + "V_BD": "V_B - V_D", + "V_CE": "V_C - V_E", + "V_DA": "V_D - V_A", + "V_EB": "V_E - V_B", + "Vterm_A": "1/5*(2*V_AC + -1*V_BD + 1*V_CE + -2*V_DA)", + "Vterm_B": "1/5*(2*V_AC + 4*V_BD + 1*V_CE + 3*V_DA)", + "Vterm_C": "1/5*(-3*V_AC + -1*V_BD + 1*V_CE + -2*V_DA)", + "Vterm_D": "1/5*(2*V_AC + -1*V_BD + 1*V_CE + 3*V_DA)", + "Vterm_E": "1/5*(-3*V_AC + -1*V_BD + -4*V_CE + -2*V_DA)", + "I_d1": "(InputCurrent(PhaseA)*cos0_1 + InputCurrent(PhaseB)*cos1_1 + InputCurrent(PhaseC)*cos2_1 + InputCurrent(PhaseD)*cos3_1 + InputCurrent(PhaseE)*cos4_1) * 2/5", + "I_q1": "(InputCurrent(PhaseA)*sin0_1 + InputCurrent(PhaseB)*sin1_1 + InputCurrent(PhaseC)*sin2_1 + InputCurrent(PhaseD)*sin3_1 + InputCurrent(PhaseE)*sin4_1) * 2/5", + "I_d3": "(InputCurrent(PhaseA)*cos0_3 + InputCurrent(PhaseB)*cos1_3 + InputCurrent(PhaseC)*cos2_3 + InputCurrent(PhaseD)*cos3_3 + InputCurrent(PhaseE)*cos4_3) * 2/5", + "I_q3": "(InputCurrent(PhaseA)*sin0_3 + InputCurrent(PhaseB)*sin1_3 + InputCurrent(PhaseC)*sin2_3 + InputCurrent(PhaseD)*sin3_3 + InputCurrent(PhaseE)*sin4_3) * 2/5", + "V_d1": "(V_A*cos0_1 + V_B*cos1_1 + V_C*cos2_1 + V_D*cos3_1 + V_E*cos4_1) * 2/5", + "V_q1": "(V_A*sin0_1 + V_B*sin1_1 + V_C*sin2_1 + V_D*sin3_1 + V_E*sin4_1) * 2/5", + "V_d3": "(V_A*cos0_3 + V_B*cos1_3 + V_C*cos2_3 + V_D*cos3_3 + V_E*cos4_3) * 2/5", + "V_q3": "(V_A*sin0_3 + V_B*sin1_3 + V_C*sin2_3 + V_D*sin3_3 + V_E*sin4_3) * 2/5", + "L0d_1": "L(PhaseA,PhaseA)*cos0_1 + L(PhaseA,PhaseB)*cos1_1 + L(PhaseA,PhaseC)*cos2_1 + L(PhaseA,PhaseD)*cos3_1 + L(PhaseA,PhaseE)*cos4_1", + "L0q_1": "L(PhaseA,PhaseA)*sin0_1 + L(PhaseA,PhaseB)*sin1_1 + L(PhaseA,PhaseC)*sin2_1 + L(PhaseA,PhaseD)*sin3_1 + L(PhaseA,PhaseE)*sin4_1", + "L1d_1": "L(PhaseB,PhaseA)*cos0_1 + L(PhaseB,PhaseB)*cos1_1 + L(PhaseB,PhaseC)*cos2_1 + L(PhaseB,PhaseD)*cos3_1 + L(PhaseB,PhaseE)*cos4_1", + "L1q_1": "L(PhaseB,PhaseA)*sin0_1 + L(PhaseB,PhaseB)*sin1_1 + L(PhaseB,PhaseC)*sin2_1 + L(PhaseB,PhaseD)*sin3_1 + L(PhaseB,PhaseE)*sin4_1", + "L2d_1": "L(PhaseC,PhaseA)*cos0_1 + L(PhaseC,PhaseB)*cos1_1 + L(PhaseC,PhaseC)*cos2_1 + L(PhaseC,PhaseD)*cos3_1 + L(PhaseC,PhaseE)*cos4_1", + "L2q_1": "L(PhaseC,PhaseA)*sin0_1 + L(PhaseC,PhaseB)*sin1_1 + L(PhaseC,PhaseC)*sin2_1 + L(PhaseC,PhaseD)*sin3_1 + L(PhaseC,PhaseE)*sin4_1", + "L3d_1": "L(PhaseD,PhaseA)*cos0_1 + L(PhaseD,PhaseB)*cos1_1 + L(PhaseD,PhaseC)*cos2_1 + L(PhaseD,PhaseD)*cos3_1 + L(PhaseD,PhaseE)*cos4_1", + "L3q_1": "L(PhaseD,PhaseA)*sin0_1 + L(PhaseD,PhaseB)*sin1_1 + L(PhaseD,PhaseC)*sin2_1 + L(PhaseD,PhaseD)*sin3_1 + L(PhaseD,PhaseE)*sin4_1", + "L4d_1": "L(PhaseE,PhaseA)*cos0_1 + L(PhaseE,PhaseB)*cos1_1 + L(PhaseE,PhaseC)*cos2_1 + L(PhaseE,PhaseD)*cos3_1 + L(PhaseE,PhaseE)*cos4_1", + "L4q_1": "L(PhaseE,PhaseA)*sin0_1 + L(PhaseE,PhaseB)*sin1_1 + L(PhaseE,PhaseC)*sin2_1 + L(PhaseE,PhaseD)*sin3_1 + L(PhaseE,PhaseE)*sin4_1", + "L0d_3": "L(PhaseA,PhaseA)*cos0_3 + L(PhaseA,PhaseB)*cos1_3 + L(PhaseA,PhaseC)*cos2_3 + L(PhaseA,PhaseD)*cos3_3 + L(PhaseA,PhaseE)*cos4_3", + "L0q_3": "L(PhaseA,PhaseA)*sin0_3 + L(PhaseA,PhaseB)*sin1_3 + L(PhaseA,PhaseC)*sin2_3 + L(PhaseA,PhaseD)*sin3_3 + L(PhaseA,PhaseE)*sin4_3", + "L1d_3": "L(PhaseB,PhaseA)*cos0_3 + L(PhaseB,PhaseB)*cos1_3 + L(PhaseB,PhaseC)*cos2_3 + L(PhaseB,PhaseD)*cos3_3 + L(PhaseB,PhaseE)*cos4_3", + "L1q_3": "L(PhaseB,PhaseA)*sin0_3 + L(PhaseB,PhaseB)*sin1_3 + L(PhaseB,PhaseC)*sin2_3 + L(PhaseB,PhaseD)*sin3_3 + L(PhaseB,PhaseE)*sin4_3", + "L2d_3": "L(PhaseC,PhaseA)*cos0_3 + L(PhaseC,PhaseB)*cos1_3 + L(PhaseC,PhaseC)*cos2_3 + L(PhaseC,PhaseD)*cos3_3 + L(PhaseC,PhaseE)*cos4_3", + "L2q_3": "L(PhaseC,PhaseA)*sin0_3 + L(PhaseC,PhaseB)*sin1_3 + L(PhaseC,PhaseC)*sin2_3 + L(PhaseC,PhaseD)*sin3_3 + L(PhaseC,PhaseE)*sin4_3", + "L3d_3": "L(PhaseD,PhaseA)*cos0_3 + L(PhaseD,PhaseB)*cos1_3 + L(PhaseD,PhaseC)*cos2_3 + L(PhaseD,PhaseD)*cos3_3 + L(PhaseD,PhaseE)*cos4_3", + "L3q_3": "L(PhaseD,PhaseA)*sin0_3 + L(PhaseD,PhaseB)*sin1_3 + L(PhaseD,PhaseC)*sin2_3 + L(PhaseD,PhaseD)*sin3_3 + L(PhaseD,PhaseE)*sin4_3", + "L4d_3": "L(PhaseE,PhaseA)*cos0_3 + L(PhaseE,PhaseB)*cos1_3 + L(PhaseE,PhaseC)*cos2_3 + L(PhaseE,PhaseD)*cos3_3 + L(PhaseE,PhaseE)*cos4_3", + "L4q_3": "L(PhaseE,PhaseA)*sin0_3 + L(PhaseE,PhaseB)*sin1_3 + L(PhaseE,PhaseC)*sin2_3 + L(PhaseE,PhaseD)*sin3_3 + L(PhaseE,PhaseE)*sin4_3", + "Ld1": "(L0d_1*cos0_1 + L1d_1*cos1_1 + L2d_1*cos2_1 + L3d_1*cos3_1 + L4d_1*cos4_1) * 2/5", + "Ld1q1": "(L0d_1*sin0_1 + L1d_1*sin1_1 + L2d_1*sin2_1 + L3d_1*sin3_1 + L4d_1*sin4_1) * 2/5", + "Lq1d1": "(L0q_1*cos0_1 + L1q_1*cos1_1 + L2q_1*cos2_1 + L3q_1*cos3_1 + L4q_1*cos4_1) * 2/5", + "Lq1": "(L0q_1*sin0_1 + L1q_1*sin1_1 + L2q_1*sin2_1 + L3q_1*sin3_1 + L4q_1*sin4_1) * 2/5", + "Ld1d3": "(L0d_1*cos0_3 + L1d_1*cos1_3 + L2d_1*cos2_3 + L3d_1*cos3_3 + L4d_1*cos4_3) * 2/5", + "Ld1q3": "(L0d_1*sin0_3 + L1d_1*sin1_3 + L2d_1*sin2_3 + L3d_1*sin3_3 + L4d_1*sin4_3) * 2/5", + "Lq1d3": "(L0q_1*cos0_3 + L1q_1*cos1_3 + L2q_1*cos2_3 + L3q_1*cos3_3 + L4q_1*cos4_3) * 2/5", + "Lq1q3": "(L0q_1*sin0_3 + L1q_1*sin1_3 + L2q_1*sin2_3 + L3q_1*sin3_3 + L4q_1*sin4_3) * 2/5", + "Ld3d1": "(L0d_3*cos0_1 + L1d_3*cos1_1 + L2d_3*cos2_1 + L3d_3*cos3_1 + L4d_3*cos4_1) * 2/5", + "Ld3q1": "(L0d_3*sin0_1 + L1d_3*sin1_1 + L2d_3*sin2_1 + L3d_3*sin3_1 + L4d_3*sin4_1) * 2/5", + "Lq3d1": "(L0q_3*cos0_1 + L1q_3*cos1_1 + L2q_3*cos2_1 + L3q_3*cos3_1 + L4q_3*cos4_1) * 2/5", + "Lq3q1": "(L0q_3*sin0_1 + L1q_3*sin1_1 + L2q_3*sin2_1 + L3q_3*sin3_1 + L4q_3*sin4_1) * 2/5", + "Ld3": "(L0d_3*cos0_3 + L1d_3*cos1_3 + L2d_3*cos2_3 + L3d_3*cos3_3 + L4d_3*cos4_3) * 2/5", + "Ld3q3": "(L0d_3*sin0_3 + L1d_3*sin1_3 + L2d_3*sin2_3 + L3d_3*sin3_3 + L4d_3*sin4_3) * 2/5", + "Lq3d3": "(L0q_3*cos0_3 + L1q_3*cos1_3 + L2q_3*cos2_3 + L3q_3*cos3_3 + L4q_3*cos4_3) * 2/5", + "Lq3": "(L0q_3*sin0_3 + L1q_3*sin1_3 + L2q_3*sin2_3 + L3q_3*sin3_3 + L4q_3*sin4_3) * 2/5", + "Flux_e_d1": "Flux_d1 - (Ld1*I_d1 + Ld1q1*I_q1 + Ld1d3*I_d3 + Ld1q3*I_q3)", + "Flux_e_q1": "Flux_q1 - (Lq1d1*I_d1 + Lq1*I_q1 + Lq1d3*I_d3 + Lq1q3*I_q3)", + "Flux_e_d3": "Flux_d3 - (Ld3d1*I_d1 + Ld3q1*I_q1 + Ld3*I_d3 + Ld3q3*I_q3)", + "Flux_e_q3": "Flux_q3 - (Lq3d1*I_d1 + Lq3q1*I_q1 + Lq3d3*I_d3 + Lq3*I_q3)", + "Torque_dq": "5/2*PolePairs*(1*(Flux_d1*I_q1 - Flux_q1*I_d1) + 3*(Flux_d3*I_q3 - Flux_q3*I_d3))", + } + + def set_post_params(self): + self.post_params = { # reports + ("InducedVoltage(PhaseA)", "InducedVoltage(PhaseB)", "InducedVoltage(PhaseC)", "InducedVoltage(PhaseD)", "InducedVoltage(PhaseE)"): "InducedVoltage", + ("Moving1.Torque"): "Torque", + ("InputCurrent(PhaseA)", "InputCurrent(PhaseB)", "InputCurrent(PhaseC)", "InputCurrent(PhaseD)", "InputCurrent(PhaseE)"): "Current", + ("FluxLinkage(PhaseA)", "FluxLinkage(PhaseB)", "FluxLinkage(PhaseC)", "FluxLinkage(PhaseD)", "FluxLinkage(PhaseE)"): "FluxLinkage", + ("I_d1", "I_q1", "I_d3", "I_q3"): "Current_dq", + ("Flux_d1", "Flux_q1", "Flux_d3", "Flux_q3"): "FluxLinkage_dq", + ("Flux_e_d1", "Flux_e_q1", "Flux_e_d3", "Flux_e_q3"): "FluxLinkage excitation_dq", + ("Vind_d1", "Vind_q1", "Vind_d3", "Vind_q3"): "InducedVoltage_dq", + ("V_d1", "V_q1", "V_d3", "V_q3"): "TerminalVoltage_dq", + ("Ld1", "Lq1", "Ld3", "Lq3"): "Inductance_dq main", + ("Ld1q1", "Ld1d3", "Ld1q3", "Lq1d3", "Lq1q3", "Ld3q3"): "Inductance_dq cross-coupling", + } + + def assign_stator_coils(self): + m2d = self.m2d + + # Excitations + I_A = "Im1*cos(2*pi*f*time+epsI1-pi) + Im3*cos(3*(2*pi*f*time)+epsI3-pi)" + I_B = "Im1*cos(2*pi*f*time-72deg+epsI1-pi) + Im3*cos(3*(2*pi*f*time-72deg)+epsI3-pi)" + I_C = "Im1*cos(2*pi*f*time-144deg+epsI1-pi) + Im3*cos(3*(2*pi*f*time-144deg)+epsI3-pi)" + I_D = "Im1*cos(2*pi*f*time-216deg+epsI1-pi) + Im3*cos(3*(2*pi*f*time-216deg)+epsI3-pi)" + I_E = "Im1*cos(2*pi*f*time-288deg+epsI1-pi) + Im3*cos(3*(2*pi*f*time-288deg)+epsI3-pi)" + m2d.assign_coil + # Define phase windings + m2d.assign_coil( + assignment=["Coil"], + conductors_number="Nc", + polarity="Positive", + name="CS1", + ) + m2d.assign_coil( + assignment=["Coil_1"], + conductors_number="Nc", + polarity="Negative", + name="CS2", + ) + m2d.assign_coil( + assignment=["Coil_2"], + conductors_number="Nc", + polarity="Negative", + name="CS3", + ) + m2d.assign_coil( + assignment=["Coil_3"], + conductors_number="Nc", + polarity="Positive", + name="CS4", + ) + m2d.assign_coil( + assignment=["Coil_4"], + conductors_number="Nc", + polarity="Positive", + name="CS5", + ) + m2d.assign_coil( + assignment=["Coil_5"], + conductors_number="Nc", + polarity="Negative", + name="CS6", + ) + m2d.assign_coil( + assignment=["Coil_6"], + conductors_number="Nc", + polarity="Negative", + name="CS7", + ) + m2d.assign_coil( + assignment=["Coil_7"], + conductors_number="Nc", + polarity="Positive", + name="CS8", + ) + m2d.assign_coil( + assignment=["Coil_8"], + conductors_number="Nc", + polarity="Positive", + name="CS9", + ) + m2d.assign_coil( + assignment=["Coil_9"], + conductors_number="Nc", + polarity="Negative", + name="CS10", + ) + + m2d.assign_winding( + assignment=None, + winding_type="Current", + is_solid=False, + current=I_A, + parallel_branches="ParallelPaths", + name="PhaseA", + ) + m2d.assign_winding( + assignment=None, + winding_type="Current", + is_solid=False, + current=I_B, + parallel_branches="ParallelPaths", + name="PhaseB", + ) + m2d.assign_winding( + assignment=None, + winding_type="Current", + is_solid=False, + current=I_C, + parallel_branches="ParallelPaths", + name="PhaseC", + ) + m2d.assign_winding( + assignment=None, + winding_type="Current", + is_solid=False, + current=I_D, + parallel_branches="ParallelPaths", + name="PhaseD", + ) + m2d.assign_winding( + assignment=None, + winding_type="Current", + is_solid=False, + current=I_E, + parallel_branches="ParallelPaths", + name="PhaseE", + ) + + m2d.add_winding_coils(assignment="PhaseA", coils=["CS1", "CS10"]) + m2d.add_winding_coils(assignment="PhaseB", coils=["CS4", "CS5"]) + m2d.add_winding_coils(assignment="PhaseC", coils=["CS8", "CS9"]) + m2d.add_winding_coils(assignment="PhaseD", coils=["CS2", "CS3"]) + m2d.add_winding_coils(assignment="PhaseE", coils=["CS6", "CS7"]) + + def inductance_computation(self): + self.m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=True) + + def set_variables(self, Id1, Iq1, Id3, Iq3): + self.m2d.variable_manager["Id1"] = f"{Id1}A" + self.m2d.variable_manager["Iq1"] = f"{Iq1}A" + self.m2d.variable_manager["Id3"] = f"{Id3}A" + self.m2d.variable_manager["Iq3"] = f"{Iq3}A" + + def extract_results(self, solutions): + # TODO: this works only because torque is the last one in the array + out = np.zeros(len(self.solution_expressions)) + for i, expr in enumerate(self.solution_expressions): + data = solutions.data_real(expr) + val = float(np.mean(data[:-1])) + + if expr.startswith("L_"): + val /= 1e9 + + out[i] = val + return data diff --git a/machine_design/generators.py b/machine_design/generators.py index dc98a6f..0798ccb 100644 --- a/machine_design/generators.py +++ b/machine_design/generators.py @@ -35,9 +35,7 @@ def save_params(params, file_name: str) -> None: class BarrierGenerator(ABC): # TODO: rename offset - def __init__( - self, design: Design, r_stator_end: float, offset: None | float = None, n_curve: int = 500, n_flat: int = 20 - ) -> None: + def __init__(self, design: Design, r_stator_end: float, offset: None | float = None, n_curve: int = 500, n_flat: int = 20) -> None: self.r_max = design.rotor_r_max self.r_min = design.rotor_r_min self.R = self.r_max - r_stator_end @@ -259,12 +257,7 @@ def _get_bezier_curve(self, phi_deg, is_inner, i): # Formula 33: Polynomial expansion to generate continuous points z_vals = np.linspace(0, 1, self.n_curve)[:, None] - return ( - (1 - z_vals) ** 3 * r0 - + 3 * (1 - z_vals) ** 2 * z_vals * r1 - + 3 * (1 - z_vals) * z_vals**2 * r2 - + z_vals**3 * r3 - ) + return (1 - z_vals) ** 3 * r0 + 3 * (1 - z_vals) ** 2 * z_vals * r1 + 3 * (1 - z_vals) * z_vals**2 * r2 + z_vals**3 * r3 class HacklGenerator_OneLambda(AbstractHacklGenerator): diff --git a/machine_design/optimization.py b/machine_design/optimization.py index 8f9aa98..9f18383 100644 --- a/machine_design/optimization.py +++ b/machine_design/optimization.py @@ -40,7 +40,7 @@ def objective_single(X: Tensor, design, generator, bounds, NUM_CORES, **kwargs) # Compute the torque try: - Tor = design.compute(NUM_CORES) + Tor = design.compute(NUM_CORES=NUM_CORES) TorAvg, _, TorRippleRms = analyze_results(Tor) except Exception: TorAvg, TorRippleRms = np.nan, np.nan diff --git a/notebooks/analyze.ipynb b/notebooks/analyze.ipynb index 06b6081..7e8f3fb 100644 --- a/notebooks/analyze.ipynb +++ b/notebooks/analyze.ipynb @@ -184,27 +184,27 @@ " for root in roots:\n", " key = method, use_constraints, root\n", " Y_pareto, hypervolume = get_pareto_hv(Y[key], hypervolume_fun)\n", - " plt.scatter([i - 0.1 + 0.2*np.random.random()], [hypervolume], color=\"blue\")\n", + " plt.scatter([i - 0.1 + 0.2 * np.random.random()], [hypervolume], color=\"blue\")\n", " i += 1\n", " x_labels.append(label_map.get(method, method))\n", "\n", "plt.xticks(range(len(x_labels)), x_labels, rotation=0)\n", "plt.ylabel(\"Hypervolume\")\n", - " \n", + "\n", "plt.axvline(2.5, color=\"gray\", linestyle=\"--\", alpha=0.6)\n", "plt.grid(True, alpha=0.3)\n", - " \n", + "\n", "ax = plt.gca()\n", "ax.set_axisbelow(True)\n", - " \n", + "\n", "ymin, ymax = plt.ylim()\n", "offset = 0.03 * (ymax - ymin)\n", - " \n", + "\n", "plt.text(1, ymax + offset, \"With constraints\", ha=\"center\", va=\"bottom\", fontsize=11)\n", "plt.text(4, ymax + offset, \"Without constraints\", ha=\"center\", va=\"bottom\", fontsize=11)\n", - " \n", + "\n", "plt.ylim(ymin, ymax + 3 * offset)\n", - " \n", + "\n", "plt.savefig(\"hypervolume_runs.png\", dpi=600, bbox_inches=\"tight\")\n", "plt.show()" ] diff --git a/notebooks/run.py b/notebooks/run.py index 09b8c06..1686e05 100644 --- a/notebooks/run.py +++ b/notebooks/run.py @@ -56,7 +56,8 @@ design.add_rotor_barrier(barrier) # Compute the torque - Tor = design.compute(num_cores) + Tor = design.compute(NUM_CORES=num_cores) + # Tor = design.compute(num_cores) if Tor is None: TorAvg, TorRippleRms = np.nan, np.nan else: diff --git a/pyproject.toml b/pyproject.toml index aafe508..9a8a1f5 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -37,7 +37,7 @@ packages = [ ] [tool.ruff] -line-length = 120 +line-length = 250 [tool.ruff.lint] select = ["E", "F", "I", "UP"]