diff --git a/.gitignore b/.gitignore index 0584a77..045eab6 100644 --- a/.gitignore +++ b/.gitignore @@ -129,4 +129,8 @@ dmypy.json .pyre/ *.csv *.fem -*.ans \ No newline at end of file +*.ans + +# Outputs +*.xlsx +*.sql diff --git a/LICENSE b/LICENSE index 261eeb9..4c9ad98 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ - Apache License + Apache License Version 2.0, January 2004 http://www.apache.org/licenses/ diff --git a/README.md b/README.md index f4c3a96..7bd4bdd 100644 --- a/README.md +++ b/README.md @@ -7,4 +7,6 @@ To run the code, follow these steps 2. In an anaconda environment, type pip install pyfemm - conda install openmdao \ No newline at end of file + conda install openmdao + +The technical support of David Meeker, Ph.D., author of the Femm library is gratefully acknowledged \ No newline at end of file diff --git a/examples/convergence_plots.py b/examples/convergence_plots.py new file mode 100755 index 0000000..3fe0a80 --- /dev/null +++ b/examples/convergence_plots.py @@ -0,0 +1,79 @@ +import openmdao.api as om +import os +import matplotlib.pyplot as plt + + +class Convergence_Trends_Opt(om.ExplicitComponent): + """ + Deprecating this for now and using OptView from PyOptSparse instead. + """ + + def initialize(self): + + self.options.declare("opt_options") + + def compute(self, inputs, outputs): + + folder_output = self.options["opt_options"]["general"]["folder_output"] + print("Figures saved here: ", folder_output) + optimization_log = os.path.join(folder_output, self.options["opt_options"]["recorder"]["file_name"]) + + if os.path.exists(optimization_log): + cr = om.CaseReader(optimization_log) + cases = cr.get_cases() + rec_data = {} + design_vars = {} + responses = {} + iterations = [] + for i, it_data in enumerate(cases): + iterations.append(i) + + # Collect DVs and responses separately for DOE + for design_var in [it_data.get_design_vars()]: + for dv in design_var: + if i == 0: + design_vars[dv] = [] + design_vars[dv].append(design_var[dv]) + + for response in [it_data.get_responses()]: + for resp in response: + if i == 0: + responses[resp] = [] + responses[resp].append(response[resp]) + + # parameters = it_data.get_responses() + for parameters in [it_data.get_responses(), it_data.get_design_vars()]: + for j, param in enumerate(parameters.keys()): + if i == 0: + rec_data[param] = [] + rec_data[param].append(parameters[param]) + + for param in rec_data.keys(): + fig, ax = plt.subplots(1, 1, figsize=(5.3, 4)) + ax.plot(iterations, rec_data[param]) + ax.set(xlabel="Number of Iterations", ylabel=param) + fig_name = "Convergence_trend_" + param + ".png" + fig.savefig(os.path.join(folder_output, fig_name)) + plt.close(fig) + + +class PlotRecorder(om.Group): + def initialize(self): + self.options.declare("opt_options") + + def setup(self): + self.add_subsystem("conv_plots", Convergence_Trends_Opt(opt_options=self.options["opt_options"])) + + +mydir = os.path.dirname(os.path.realpath(__file__)) # get path to this file +output_folder_name = "test5" + +opt_options = {} +opt_options["general"] = {} +opt_options["general"]["folder_output"] = os.path.join(mydir, "outputs", output_folder_name) +opt_options["recorder"] = {} +opt_options["recorder"]["file_name"] = "log.sql" + +wt_opt = om.Problem(model=PlotRecorder(opt_options=opt_options)) +wt_opt.setup(derivatives=False) +wt_opt.run_model() diff --git a/examples/run_LTS.py b/examples/run_LTS.py old mode 100644 new mode 100755 index 2e5d05c..7dec9a1 --- a/examples/run_LTS.py +++ b/examples/run_LTS.py @@ -1,43 +1,94 @@ import openmdao.api as om from lts.lts import LTS_Outer_Rotor_Opt +from lts.structural import LTS_Outer_Rotor_Structural import os import pandas as pd +mydir = os.path.dirname(os.path.realpath(__file__)) # get path to this file + def cleanup_femm_files(clean_dir): files = os.listdir(clean_dir) for file in files: if file.endswith(".ans") or file.endswith(".fem") or file.endswith(".csv"): os.remove(os.path.join(clean_dir, file)) -if __name__ == "__main__": +def copy_data(prob_in, prob_out): + + # Get all OpenMDAO inputs and outputs into a dictionary + def create_dict(prob): + dict_omdao = prob.model.list_inputs(val=True, hierarchical=False, prom_name=True, units=False, desc=False, out_stream=None) + temp = prob.model.list_outputs(val=True, hierarchical=False, prom_name=True, units=False, desc=False, out_stream=None) + dict_omdao.extend(temp) + my_dict = {} + for k in range(len(dict_omdao)): + my_dict[ dict_omdao[k][1]["prom_name"] ] = dict_omdao[k][1]["val"] + + return my_dict + + in_dict = create_dict(prob_in) + out_dict = create_dict(prob_out) + + for k in in_dict: + if k in out_dict: + prob_out[k] = in_dict[k] + + return prob_out - mydir = os.path.dirname(os.path.realpath(__file__)) # get path to this file - output_dir = os.path.join(os.path.dirname(os.path.dirname(mydir)), 'outputs', 'test1') +def save_data(fname, prob): + # Remove file extension + froot = os.path.splitext(fname)[0] + + # Get all OpenMDAO inputs and outputs into a dictionary + var_dict = prob.model.list_inputs(val=True, hierarchical=True, prom_name=True, units=True, desc=True, out_stream=None) + out_dict = prob.model.list_outputs(val=True, hierarchical=True, prom_name=True, units=True, desc=True, out_stream=None) + var_dict.extend(out_dict) + + data = {} + data["variables"] = [] + data["units"] = [] + data["values"] = [] + data["description"] = [] + for k in range(len(var_dict)): + unit_str = var_dict[k][1]["units"] + if unit_str is None: + unit_str = "" + + iname = var_dict[k][1]["prom_name"] + if iname in data["variables"]: + continue + + data["variables"].append(iname) + data["units"].append(unit_str) + data["values"].append(var_dict[k][1]["val"]) + data["description"].append(var_dict[k][1]["desc"]) + df = pd.DataFrame(data) + df.to_excel(froot + ".xlsx", index=False) + df.to_csv(froot + ".csv", index=False) + +def optimize_magnetics_design(prob_in=None, output_dir=None, cleanup_flag=True): + if output_dir is None: + output_dir = "outputs" os.makedirs(output_dir, exist_ok=True) modeling_options = {} - modeling_options['output_dir'] = output_dir + modeling_options["output_dir"] = output_dir - cleanup_flag = True # Clean run directory before the run if cleanup_flag: cleanup_femm_files(mydir) prob = om.Problem() - prob.model = LTS_Outer_Rotor_Opt(modeling_options = modeling_options) - - prob.driver = om.ScipyOptimizeDriver() # pyOptSparseDriver() - prob.driver.options['optimizer'] = 'SLSQP' #'COBYLA' #' - prob.driver.options["maxiter"] = 10 - # prob.driver.opt_settings['IPRINT'] = 4 - # prob.driver.opt_settings['ITRM'] = 3 - # prob.driver.opt_settings['ITMAX'] = 10 - # prob.driver.opt_settings['DELFUN'] = 1e-3 - # prob.driver.opt_settings['DABFUN'] = 1e-3 - # prob.driver.opt_settings['IFILE'] = 'CONMIN_LST.out' - # prob.root.deriv_options['type']='fd' - - recorder = om.SqliteRecorder(os.path.join(output_dir,"log.sql")) + prob.model = LTS_Outer_Rotor_Opt(modeling_options=modeling_options) + + #prob.driver = om.ScipyOptimizeDriver() + #prob.driver.options['optimizer'] = 'COBYLA' #'SLSQP' # + #prob.driver.options["maxiter"] = 500 #50 + prob.driver = om.DifferentialEvolutionDriver() + prob.driver.options["max_gen"] = 10 + prob.driver.options["pop_size"] = 40 + prob.driver.options["penalty_exponent"] = 2 + + recorder = om.SqliteRecorder(os.path.join(output_dir, "log.sql")) prob.driver.add_recorder(recorder) prob.add_recorder(recorder) prob.driver.recording_options["excludes"] = ["*_df"] @@ -45,135 +96,268 @@ def cleanup_femm_files(clean_dir): prob.driver.recording_options["record_desvars"] = True prob.driver.recording_options["record_objectives"] = True - prob.model.add_design_var("D_a", lower=7, upper=8, ref=7.5) + prob.model.add_design_var("D_a", lower=5, upper=9, ref=8) prob.model.add_design_var("delta_em", lower=0.060, upper=0.10, ref=0.08) - prob.model.add_design_var("h_sc", lower=0.03, upper=0.25, ref=0.01) + prob.model.add_design_var("h_sc", lower=0.03, upper=0.15, ref=0.06) prob.model.add_design_var("h_s", lower=0.1, upper=0.4, ref=0.1) - prob.model.add_design_var("p", lower=10, upper=30, ref=20) + prob.model.add_design_var("p", lower=20, upper=30, ref=25) prob.model.add_design_var("h_yr", lower=0.01, upper=0.4, ref=0.1) - prob.model.add_design_var("l_s", lower=1, upper=2.5, ref=1.625) - prob.model.add_design_var("alpha", lower=0.5, upper=20, ref=10) - prob.model.add_design_var("dalpha", lower=1, upper=10, ref=10) - # prob.model.add_design_var('beta', lower=0.75, upper=20,ref=10) + prob.model.add_design_var("l_s", lower=1, upper=1.75) + prob.model.add_design_var("alpha", lower=0.1, upper=1) + prob.model.add_design_var("dalpha", lower=1, upper=5) prob.model.add_design_var("I_sc", lower=200, upper=700, ref=450) - prob.model.add_design_var("N_sc", lower=1500, upper=2500, ref=1500) - prob.model.add_design_var("N_c", lower=2, upper=30, ref=16) + prob.model.add_design_var("N_sc", lower=1500, upper=3000, ref=1500) + prob.model.add_design_var("N_c", lower=1, upper=15, ref=8) prob.model.add_design_var("I_s", lower=500, upper=3000, ref=1750) - #prob.model.add_design_var("J_s", lower=1.5, upper=6, ref=3.75) - prob.model.add_design_var("h_yr_s", lower=0.0250, upper=0.5, ref=0.3) - prob.model.add_design_var("h_ys", lower=0.025, upper=0.6, ref=0.35) - prob.model.add_design_var("t_rdisc", lower=0.025, upper=0.5, ref=0.3) - prob.model.add_design_var("t_sdisc", lower=0.025, upper=0.5, ref=0.3) - prob.model.add_objective("mass_total", ref=1e6) - #prob.model.add_objective("l_sc", ref=1e3) - - # prob.model.add_constraint('K_rad', lower=0.15,upper=0.3) #10 + # prob.model.add_design_var("J_s", lower=1.5, upper=6, ref=3.75) + + # prob.model.add_objective("mass_total", ref=1e6) + prob.model.add_objective("cost_total", ref=1e6) + # prob.model.add_constraint("Slot_aspect_ratio", lower=4.0, upper=10.0) # 11 prob.model.add_constraint("con_angle", lower=0.001) - prob.model.add_constraint("con_angle2", lower=0.001) - #prob.model.add_constraint("E_p", lower=3300, upper=3350) - #prob.model.add_constraint("con_N_sc", lower=-5, upper=5) - prob.model.add_constraint("B_rymax", upper=2.1) - prob.model.add_constraint("Torque_actual", lower=23.07e6, ref=1e7) - prob.model.add_constraint("Critical_current_ratio",upper=1.) - prob.model.add_constraint("Coil_max_ratio",upper=1.) + # Differential evolution driver cannot do double-sided constraints, so have to hack it + prob.model.add_constraint("E_p", lower=0.8 * 3300, ref=3000) + prob.model.add_constraint("E_p_ratio", upper=1.20) + prob.model.add_constraint("B_coil_max", lower=6.0) + prob.model.add_constraint("Coil_max_ratio", upper=1.2) # Consider user-defined limit instead of load line - prob.model.add_constraint("U_rotor_radial_constraint", lower=0.01) - prob.model.add_constraint("U_rotor_axial_constraint", lower=0.01) - prob.model.add_constraint("U_stator_radial_constraint", lower=0.01) - prob.model.add_constraint("U_stator_axial_constraint", lower=0.01) + prob.model.add_constraint("B_rymax", upper=2.1) + + prob.model.add_constraint("gen_eff", lower=0.97) + prob.model.add_constraint("torque_ratio", lower=1.0) + prob.model.add_constraint("Torque_actual", upper=24.5e6, ref=20e6) + # prob.model.add_constraint("Critical_current_ratio",upper=1.) prob.model.approx_totals(method="fd") prob.setup() # --- Design Variables --- - # UNUSED - # Assign values to universal constants - #prob["B_r"] = 1.279 # Tesla remnant flux density - #prob["E"] = 2e11 # N/m^2 young's modulus - #prob["ratio"] = 0.8 # ratio of magnet width to pole pitch(bm/self.tau_p) - #prob["mu_0"] = np.pi * 4e-7 # permeability of free space - - #prob["mu_r"] = 1.06 # relative permeability - #prob["cofi"] = 0.85 # power factor - - # Assign values to design constants - #prob["h_0"] = 0.005 # Slot opening height - #prob["h_1"] = 0.004 # Slot wedge height - #prob["k_sfil"] = 0.65 # Slot fill factor - #prob["P_Fe0h"] = 4 # specific hysteresis losses W/kg @ 1.5 T - #prob["P_Fe0e"] = 1 # specific hysteresis losses W/kg @ 1.5 T - #prob["k_fes"] = 0.8 # Iron fill factor - - # Assign values to universal constants - #prob["gravity"] = 9.8106 # m/s**2 acceleration due to gravity - #prob["E"] = 2e11 # Young's modulus - #prob["phi"] = 90 * 2 * np.pi / 360 # tilt angle (rotor tilt -90 degrees during transportation) - #prob["v"] = 0.3 # Poisson's ratio - #prob["G"] = 79.3e9 - prob["m"] = 6 # phases - prob["q"] = 2 # slots per pole - prob["b_s_tau_s"] = 0.45 - prob["conductor_area"] = 1.8 * 1.2e-6 - prob["K_h"] = 2.0 - prob["K_e"] = 0.5 - - # Initial design variables for a PMSG designed for a 15MW turbine - prob["P_rated"] = 17e6 - prob["T_rated"] = 23.07e6 # rev 1 9.94718e6 - prob["N_nom"] = 7.7 # 7.5598598 #8.68 # rpm 9.6 - prob["l_s"] = 1.00390095 # 8.68 # rpm 9.6 - prob["D_a"] = 7.74736313 # rev 1 6.8 - prob["delta_em"] = 0.02 # rev 2.1 - prob["h_s"] = 0.1803019703 # rev 1 0.3 - prob["p"] = 21 # 100.0 # rev 1 160 - prob["h_sc"] = 0.0503409354 # rev 1 0.034 - prob["h_yr"] = 0.1254730934 # rev 1 0.045 - prob["alpha"] = 0.53805442 # rev 1 0.045 - prob["dalpha"] = 1.0 # rev 1 0.045 - # prob['beta'] = 1.75 # rev 1 0.045 - prob["I_sc"] = 400.4427393 # rev 1 0.045 - prob["N_sc"] = 1502 # rev 1 0.045 - prob["N_c"] = 2 - prob["I_s"] = 2995.06090335 - prob["J_s"] = 3.0 - - # Material properties - prob["rho_steel"] = 7850 - prob["rho_Fe"] = 7700.0 # Steel density - prob["rho_Copper"] = 8900.0 # Kg/m3 copper density - prob["rho_NbTi"] = 8442.37093661195 # magnet density - prob["rho_Cu"] = 1.8e-8 * 1.4 # Copper resisitivty - prob["U_b"] = 2 # brush contact voltage - prob["Y"] = 10 #Short pitch - - prob["Tilt_angle"] = 6.0 - prob["R_shaft_outer"] = 1.25 - prob["R_nose_outer"] = 0.95 - prob["u_allow_pcent"] = 50 - prob["y_allow_pcent"] = 20 - #prob["h_yr"] = 0.1254730934 - prob["h_yr_s"] = 0.025 - prob["h_ys"] = 0.050 - prob["t_rdisc"] = 0.05 - prob["t_sdisc"] = 0.100 - prob["y_bd"] = 0.00 - prob["theta_bd"] = 0.00 - prob["y_sh"] = 0.00 - prob["theta_sh"] = 0.00 - - #prob.model.approx_totals(method="fd") - - #prob.run_model() - prob.run_driver() + if prob_in is None: + prob["m"] = 6 # phases + prob["q"] = 2 # slots per pole per phase + prob["b_s_tau_s"] = 0.45 + prob["conductor_area"] = 1.8 * 1.2e-6 + prob["K_h"] = 2 # specific hysteresis losses W/kg @ 1.5 T + prob["K_e"] = 0.5 # specific hysteresis losses W/kg @ 1.5 T + + # ## Initial design variables for a PMSG designed for a 15MW turbine + prob["P_rated"] = 17e6 + prob["T_rated"] = 23.07e6 + prob["E_p_target"] = 3300.0 + prob["N_nom"] = 7.7 + prob["D_a"] = 9.0 #7.86277266 + prob["delta_em"] = 0.07114809 #0.096 + prob["h_s"] = 0.26448025 #0.1 + prob["p"] = 30.0 + prob["h_sc"] = 0.0798385 #0.15 + prob["h_yr"] = 0.4 #0.18651016 + prob["alpha"] = 0.1 #0.73251198 + prob["dalpha"] = 1.0 + prob["I_sc"] = 675.23529314 #660.09589169 + prob["N_sc"] = 1500.0 #1941.1454573 + prob["N_c"] = 1.0 #2.07793977 + prob["I_s"] = 2933.37488172 #3000.0 + prob["J_s"] = 3.0 + prob["l_s"] = 1.75 #1.05499191 + + # Specific costs + prob["C_Cu"] = 10.3 # https://markets.businessinsider.com/commodities/copper-price + prob["C_Fe"] = 0.556 + prob["C_Fes"] = 0.50139 + prob["C_NbTi"] = 30.0 + + # Material properties + prob["rho_steel"] = 7850 + prob["rho_Fe"] = 7700.0 # Steel density + prob["rho_Copper"] = 8900.0 # Kg/m3 copper density + prob["rho_NbTi"] = 8442.37093661195 # magnet density + prob["resisitivty_Cu"] = 1.724e-8 # 1.8e-8 * 1.4 # Copper resisitivty + prob["U_b"] = 1 # brush voltage drop + prob["Y"] = 10 # Short pitch + + prob["Tilt_angle"] = 90.0 + prob["R_shaft_outer"] = 1.25 + prob["R_nose_outer"] = 0.95 + prob["u_allow_pcent"] = 30 + prob["y_allow_pcent"] = 20 + prob["h_yr_s"] = 0.18866371 #0.43394501 + prob["h_ys"] = 0.025 #19642297 + prob["t_rdisc"] = 0.5 #19016171 + prob["t_sdisc"] = 0.5 #42709254 + prob["y_bd"] = 0.00 + prob["theta_bd"] = 0.00 + prob["y_sh"] = 0.00 + prob["theta_sh"] = 0.00 + + else: + prob = copy_data(prob_in, prob) + + prob.model.approx_totals(method="fd") + + prob.run_model() + #prob.run_driver() # Clean run directory after the run if cleanup_flag: cleanup_femm_files(mydir) - prob.model.list_outputs(values = True, hierarchical=True) + prob.model.list_outputs(val=True, hierarchical=True) + print("Final solution:") + print("E_p_ratio", prob["E_p_ratio"]) + # print("con_angle", prob["con_angle"]) + print("gen_eff", prob["gen_eff"]) + print("N_c", prob["N_c"]) + print("N_sc", prob["N_sc"]) + print("B_coil_max", prob["B_coil_max"]) + print("l_s", prob["l_s"]) + print("Torque_actual", prob["Torque_actual"]) + + return prob + +def optimize_structural_design(prob_in=None, output_dir=None): + if output_dir is None: + output_dir = "outputs" + os.makedirs(output_dir, exist_ok=True) + + prob_struct = om.Problem() + prob_struct.model = LTS_Outer_Rotor_Structural() + + prob_struct.driver = om.ScipyOptimizeDriver() + prob_struct.driver.options["optimizer"] = "SLSQP" + prob_struct.driver.options["maxiter"] = 100 + + #recorder = om.SqliteRecorder("log.sql") + #prob_struct.driver.add_recorder(recorder) + #prob_struct.add_recorder(recorder) + #prob_struct.driver.recording_options["excludes"] = ["*_df"] + #prob_struct.driver.recording_options["record_constraints"] = True + #prob_struct.driver.recording_options["record_desvars"] = True + #prob_struct.driver.recording_options["record_objectives"] = True + + prob_struct.model.add_design_var("h_yr_s", lower=0.0250, upper=0.5, ref=0.1) + prob_struct.model.add_design_var("h_ys", lower=0.025, upper=0.6, ref=0.1) + prob_struct.model.add_design_var("t_rdisc", lower=0.025, upper=0.5, ref=0.1) + prob_struct.model.add_design_var("t_sdisc", lower=0.025, upper=0.5, ref=0.1) + prob_struct.model.add_objective("mass_structural", ref=1e6) + + prob_struct.model.add_constraint("U_rotor_radial_constraint", upper=1.0) + prob_struct.model.add_constraint("U_rotor_axial_constraint", upper=1.0) + prob_struct.model.add_constraint("U_stator_radial_constraint", upper=1.0) + prob_struct.model.add_constraint("U_stator_axial_constraint", upper=1.0) + + prob_struct.model.approx_totals(method="fd") + + prob_struct.setup() + # --- Design Variables --- + + if prob_in is None: + # Initial design variables for a PMSG designed for a 15MW turbine + #prob_struct["Sigma_shear"] = 74.99692029e3 + prob_struct["Sigma_normal"] = 378.45123826e3 + prob_struct["T_e"] = 9e06 + prob_struct["l_eff_stator"] = 1.44142189 # rev 1 9.94718e6 + prob_struct["l_eff_rotor"] = 1.2827137 + prob_struct["D_a"] = 7.74736313 + prob_struct["delta_em"] = 0.0199961 + prob_struct["h_s"] = 0.1803019703 + prob_struct["D_sc"] = 7.78735533 + prob_struct["rho_steel"] = 7850 + prob_struct["rho_Fe"] = 7700 + prob_struct["Tilt_angle"] = 90.0 + prob_struct["R_shaft_outer"] = 1.25 + prob_struct["R_nose_outer"] = 0.95 + prob_struct["u_allow_pcent"] = 50 + prob_struct["y_allow_pcent"] = 20 + prob_struct["h_yr"] = 0.1254730934 + prob_struct["h_yr_s"] = 0.025 + prob_struct["h_ys"] = 0.050 + prob_struct["t_rdisc"] = 0.05 + prob_struct["t_sdisc"] = 0.1 + prob_struct["y_bd"] = 0.00 + prob_struct["theta_bd"] = 0.00 + prob_struct["y_sh"] = 0.00 + prob_struct["theta_sh"] = 0.00 + + #prob_struct["mass_copper"] = 60e3 + prob_struct["mass_SC"] = 4000 + else: + prob_struct = copy_data(prob_in, prob_struct) + + prob_struct.model.approx_totals(method="fd") + + #prob_struct.run_model() + prob_struct.run_driver() + + #prob_struct.model.list_outputs(val=True, hierarchical=True) + + raw_data = { + "Parameters": [ + "Rotor disc thickness", + "Rotor yoke thickness", + "Stator disc thickness", + "Stator yoke thickness", + "Rotor radial deflection", + "Rotor axial deflection", + "Stator radial deflection", + "Stator axial deflection", + "Rotor structural mass", + "Stator structural mass", + "Total structural mass", + ], + "Values": [ + prob_struct.get_val("t_rdisc", units="mm"), + prob_struct.get_val("h_yr_s", units="mm"), + prob_struct.get_val("t_sdisc", units="mm"), + prob_struct.get_val("h_ys", units="mm"), + prob_struct.get_val("u_ar", units="mm"), + prob_struct.get_val("y_ar", units="mm"), + prob_struct.get_val("u_as", units="mm"), + prob_struct.get_val("y_as", units="mm"), + prob_struct.get_val("Structural_mass_rotor", units="t"), + prob_struct.get_val("Structural_mass_stator", units="t"), + prob_struct.get_val("mass_structural", units="t"), + ], + "Limit": [ + "", + "", + "", + "", + prob_struct.get_val("u_allowable_r", units="mm"), + prob_struct.get_val("y_allowable_r", units="mm"), + prob_struct.get_val("u_allowable_s", units="mm"), + prob_struct.get_val("y_allowable_s", units="mm"), + "", + "", + "", + ], + "Units": [ + "mm", + "mm", + "mm", + "mm", + "mm", + "mm", + "mm", + "mm", + "tons", + "tons", + "tons", + ], + } + #df = pd.DataFrame(raw_data, columns=["Parameters", "Values", "Limit", "Units"]) + #df.to_excel(os.path.join(output_dir, "Optimized_structure_LTSG_MW.xlsx")) + + return prob_struct + +def write_all_data(prob, output_dir=None): + if output_dir is None: + output_dir = "outputs" + os.makedirs(output_dir, exist_ok=True) + + save_data(os.path.join(output_dir, "LTS_output"), prob) raw_data = { "Parameters": [ @@ -221,6 +405,21 @@ def cleanup_femm_files(clean_dir): "Iron mass", "Total active mass", "Efficiency", + "Rotor disc thickness", + "Rotor yoke thickness", + "Stator disc thickness", + "Stator yoke thickness", + "Rotor radial deflection", + "Rotor axial deflection", + "Rotor torsional twist", + "Stator radial deflection", + "Stator axial deflection", + "Stator torsional twist", + "Rotor structural mass", + "Stator structural mass", + "Total structural mass", + "Total generator mass", + "Total generator cost", ], "Values": [ prob.get_val("P_rated", units="MW"), @@ -246,7 +445,7 @@ def cleanup_femm_files(clean_dir): prob["p1"], prob["E_p"], prob["I_s"], - prob["S"], + prob["Slots"], int(prob["N_c"]), prob["J_s"], prob["R_s"], @@ -258,15 +457,30 @@ def cleanup_femm_files(clean_dir): prob["N_l"], prob["N_sc_layer"], prob.get_val("l_sc", units="km"), - prob["mass_SC"], - prob.get_val("Total_mass_SC",units="t"), + prob["mass_SC_racetrack"], + prob.get_val("mass_SC", units="t"), prob["B_rymax"], prob["B_g"], prob["B_coil_max"], - prob.get_val("Copper", units="t"), - prob.get_val("Iron", units="t"), - prob.get_val("Mass", units="t"), + prob.get_val("mass_copper", units="t"), + prob.get_val("mass_iron", units="t"), + prob.get_val("mass_active", units="t"), prob["gen_eff"] * 100, + prob.get_val("t_rdisc", units="mm"), + prob.get_val("h_yr_s", units="mm"), + prob.get_val("t_sdisc", units="mm"), + prob.get_val("h_ys", units="mm"), + prob.get_val("u_ar", units="mm"), + prob.get_val("y_ar", units="mm"), + prob.get_val("twist_r"), + prob.get_val("u_as", units="mm"), + prob.get_val("y_as", units="mm"), + prob.get_val("twist_s"), + prob.get_val("Structural_mass_rotor", units="t"), + prob.get_val("Structural_mass_stator", units="t"), + prob.get_val("mass_structural", units="t"), + prob.get_val("mass_total", units="t"), + prob.get_val("cost_total", units="USD"), ], "Limit": [ "", @@ -313,6 +527,21 @@ def cleanup_femm_files(clean_dir): "", "", "", + "", + "", + "", + "", + prob.get_val("u_allowable_r", units="mm"), + prob.get_val("y_allowable_r", units="mm"), + "", + prob.get_val("u_allowable_s", units="mm"), + prob.get_val("y_allowable_s", units="mm"), + "", + "", + "", + "", + "", + "", ], "Units": [ "MW", @@ -359,19 +588,41 @@ def cleanup_femm_files(clean_dir): "Tons", "Tons", "%", + "mm", + "mm", + "mm", + "mm", + "mm", + "mm", + "deg", + "mm", + "mm", + "deg", + "tons", + "tons", + "tons", + "tons", + "$", ], } - #print(raw_data) + df = pd.DataFrame(raw_data, columns=["Parameters", "Values", "Limit", "Units"]) + df.to_excel(os.path.join(output_dir, "Optimized_LTSG_" + str(prob["P_rated"][0] / 1e6) + "_MW.xlsx")) - #print(df) - df.to_excel(os.path.join(output_dir,"Optimized_LTSG_" + str(prob["P_rated"][0] / 1e6) + "_MW.xlsx")) - print("Final solution:") - print("Slot", prob["Slot_aspect_ratio"]) - print("con_angle", prob["con_angle"]) - #print("con_I_sc", prob["con_I_sc"]) - #print("con_N_sc", prob["con_N_sc"]) - print("B_g", prob["B_g"]) - print("l_s", prob["l_s"]) - print("Torque_actual", prob["Torque_actual"]) +if __name__ == "__main__": + + output_dir = os.path.join(mydir, "outputs") + + # Optimize just magnetrics with GA and then structural with SLSQP + prob = optimize_magnetics_design(output_dir=output_dir) + prob_struct = optimize_structural_design(prob_in=prob, output_dir=output_dir) + + # Bring all data together + prob = copy_data(prob_struct, prob) + prob.run_model() + + # Write to xlsx and csv files + write_all_data(prob, output_dir=output_dir) + prob.model.list_outputs(val=True, hierarchical=True) + cleanup_femm_files(mydir) diff --git a/lts/femm_fea.py b/lts/femm_fea.py index 904f913..da52c97 100644 --- a/lts/femm_fea.py +++ b/lts/femm_fea.py @@ -9,6 +9,16 @@ import openmdao.api as om +def bad_inputs(outputs): + outputs["B_g"] = 7 + outputs["B_coil_max"] = 12 + outputs["B_rymax"] = 5 + outputs["Torque_actual"] = 1.0e6 # 50e6 + outputs["Sigma_shear"] = 1.0e3 # 300000 + outputs["Sigma_normal"] = 1.0e3 # 200000 + return outputs + + def run_post_process(D_a, radius_sc, h_sc, slot_radius, theta_p_r, alpha_r, beta_r, n): # After looking at the femm results, this function post-processes them, no electrical load condition @@ -25,10 +35,11 @@ def run_post_process(D_a, radius_sc, h_sc, slot_radius, theta_p_r, alpha_r, beta femm.mo_makeplot(1, 1500, "gap_" + str(n) + ".csv", 1) femm.mo_makeplot(2, 100, "B_r_normal" + str(n) + ".csv", 1) femm.mo_makeplot(3, 100, "B_t_normal" + str(n) + ".csv", 1) + # temp_force = femm.mo_lineintegral(3) + # temp_torque = femm.mo_lineintegral(4) femm.mo_clearcontour() # femm.mo_addcontour((radius_sc)*np.cos(beta_r),(radius_sc)*np.sin(beta_r)) - # femm.mo_addcontour((radius_sc)*np.cos(alpha_r),(radius_sc)*np.sin(alpha_r)) # femm.mo_addcontour((radius_sc+h_sc)*np.cos(alpha_r),(radius_sc+h_sc)*np.sin(alpha_r)) # femm.mo_addcontour((radius_sc+h_sc)*np.cos(beta_r),(radius_sc+h_sc)*np.sin(beta_r)) # femm.mo_addcontour((radius_sc)*np.cos(beta_r),(radius_sc)*np.sin(beta_r)) @@ -62,11 +73,11 @@ def run_post_process(D_a, radius_sc, h_sc, slot_radius, theta_p_r, alpha_r, beta B_t_normal = np.loadtxt("B_t_normal" + str(n) + ".csv") circ = B_r_normal[-1, 0] - delta_L = np.diff(B_r_normal[:, 0])[0] - - force = np.sum((B_r_normal[:, 1]) ** 2 - (B_t_normal[:, 1]) ** 2) * delta_L + force = np.trapz(B_r_normal[:, 1] ** 2 - B_t_normal[:, 1] ** 2, B_r_normal[:, 0]) sigma_n = abs(1 / (8 * np.pi * 1e-7) * force) / circ - + # print(force, force/delta_L, sigma_n, sigma_n*circ) + # print(temp_force, np.linalg.norm(temp_force)) + # print(temp_torque) # print ((radius_sc)*np.cos(beta_r),(radius_sc)*np.sin(beta_r),(radius_sc)*np.cos(alpha_r),(radius_sc)*np.sin(alpha_r), # (radius_sc+h_sc)*np.cos(alpha_r),(radius_sc+h_sc)*np.sin(alpha_r),(radius_sc+h_sc)*np.cos(beta_r), # (radius_sc+h_sc)*np.sin(beta_r),(radius_sc)*np.cos(beta_r),(radius_sc)*np.sin(beta_r)) @@ -115,34 +126,40 @@ def B_r_B_t(D_a, l_s, p1, delta_em, theta_p_r, I_s, theta_b_t, theta_b_s, layer_ count = 0 angle_r = theta_b_t * 0.5 + theta_b_s * 0.5 delta_theta = theta_b_t + theta_b_s - for pitch in range(1, Y_q, 2): - femm.mi_selectlabel( - layer_2 * np.cos(angle_r + (pitch - 1) * (delta_theta)), - layer_2 * np.sin(angle_r + (pitch - 1) * (delta_theta)), - ) - femm.mi_selectlabel( - layer_2 * np.cos(angle_r + pitch * delta_theta), layer_2 * np.sin(angle_r + pitch * delta_theta) - ) - femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c) - # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c_a1[count]) - femm.mi_clearselected() - count = count + 1 + for pitch in range(1, int(np.ceil(Y_q)), 2): + try: + femm.mi_selectlabel( + layer_2 * np.cos(angle_r + (pitch - 1) * (delta_theta)), + layer_2 * np.sin(angle_r + (pitch - 1) * (delta_theta)), + ) + femm.mi_selectlabel( + layer_2 * np.cos(angle_r + pitch * delta_theta), layer_2 * np.sin(angle_r + pitch * delta_theta) + ) + femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c) + # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c_a1[count]) + femm.mi_clearselected() + count = count + 1 + except: + continue count = 0 angle_r = theta_b_t * 0.5 + theta_b_s * 0.5 delta_theta = theta_b_t + theta_b_s - for pitch in range(1, Y_q, 2): - femm.mi_selectlabel( - layer_1 * np.cos(angle_r + (pitch - 1) * (delta_theta)), - layer_1 * np.sin(angle_r + (pitch - 1) * (delta_theta)), - ) - femm.mi_selectlabel( - layer_1 * np.cos(angle_r + pitch * delta_theta), layer_1 * np.sin(angle_r + pitch * delta_theta) - ) - femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c) - # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c_a1[count+1]) - femm.mi_clearselected() - count = count + 1 + for pitch in range(1, int(np.ceil(Y_q)), 2): + try: + femm.mi_selectlabel( + layer_1 * np.cos(angle_r + (pitch - 1) * (delta_theta)), + layer_1 * np.sin(angle_r + (pitch - 1) * (delta_theta)), + ) + femm.mi_selectlabel( + layer_1 * np.cos(angle_r + pitch * delta_theta), layer_1 * np.sin(angle_r + pitch * delta_theta) + ) + femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c) + # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c_a1[count+1]) + femm.mi_clearselected() + count = count + 1 + except: + continue femm.mi_modifycircprop("D+", 1, I_s * np.sin(T_elec + np.pi / 6)) femm.mi_modifycircprop("C-", 1, -I_s * np.sin(T_elec - 4 * np.pi / 6)) @@ -160,6 +177,9 @@ def B_r_B_t(D_a, l_s, p1, delta_em, theta_p_r, I_s, theta_b_t, theta_b_s, layer_ femm.mo_bendcontour(theta_p_d, 0.25) femm.mo_makeplot(2, 100, "B_r_2.csv", 1) femm.mo_makeplot(3, 100, "B_t_2.csv", 1) + # temp_Bn = femm.mo_lineintegral(0) + # temp_force = femm.mo_lineintegral(3) + # temp_torque = femm.mo_lineintegral(4) femm.mo_clearcontour() femm.mo_close() @@ -167,18 +187,23 @@ def B_r_B_t(D_a, l_s, p1, delta_em, theta_p_r, I_s, theta_b_t, theta_b_s, layer_ B_t_1 = np.loadtxt("B_t_1.csv") B_r_2 = np.loadtxt("B_r_2.csv") B_t_2 = np.loadtxt("B_t_2.csv") - delta_L = np.diff(B_r_1[:, 0])[0] circ = B_r_1[-1, 0] - force = np.array([np.sum(B_r_1[:, 1] * B_t_1[:, 1]), np.sum(B_r_2[:, 1] * B_t_2[:, 1])]) * delta_L + force = np.array( + [np.trapz(B_r_1[:, 1] * B_t_1[:, 1], B_r_1[:, 0]), np.trapz(B_r_2[:, 1] * B_t_2[:, 1], B_r_2[:, 0])] + ) sigma_t = abs(1 / (4 * np.pi * 1e-7) * force) / circ - torque = np.pi / 2 * sigma_t * D_a ** 2 * l_s - - #print("Torque values in MNm:") - #print(torque[0] / 1e6, torque[1] / 1e6) + torque = np.pi / 2 * sigma_t * D_a ** 2 * l_s * 1.08440860215053764 + # print(force[-1], sigma_t[-1], sigma_t[-1]*circ, torque[-1]) + # print(temp_force, np.linalg.norm(temp_force)) + # print(temp_torque) + # print(temp_Bn, np.trapz(B_r_2[:, 1], B_r_2[:, 0])) + # print("Torque values in MNm:") + # print(torque[0] / 1e6, torque[1] / 1e6) # Air gap electro-magnetic torque for the full machine # Average shear stress for the full machine + # print (torque[0],torque[1]) return torque.mean(), sigma_t.mean() @@ -210,6 +235,16 @@ def setup(self): self.add_input("N_c", 0.0, desc="Number of turns per coil") self.add_input("delta_em", 0.0, units="m", desc="airgap length ") self.add_input("I_s", 0.0, units="A", desc="Generator output phase current") + self.add_input("con_angle", 0.0, units="deg", desc="Geometry constraint status") + self.add_input("a_m", 0.0, units="m", desc="Coil separation distance") + self.add_input("W_sc", 0.0, units="m", desc="SC coil width") + self.add_input("Outer_width", 0.0, units="m", desc="Coil outer width") + self.add_input("field_coil_x", np.zeros(8), desc="Field coil points") + self.add_input("field_coil_y", np.zeros(8), desc="Field coil points") + self.add_input("field_coil_xlabel", np.zeros(2), desc="Field coil label points") + self.add_input("field_coil_ylabel", np.zeros(2), desc="Field coil label points") + self.add_input("Slots", 0.0, desc="Stator slots") + self.add_input("y_Q", 0.0, desc="Slots per pole also pole pitch") # Outputs self.add_output("B_g", 0.0, desc="Peak air gap flux density ") @@ -225,9 +260,6 @@ def setup(self): 0.0, desc="Ratio of actual to critical coil flux density, usually constrained to be smaller than 1", ) - self.add_output("a_m", 0.0, units="m", desc="Coil separation distance") - self.add_output("W_sc", 0.0, units="m", desc="SC coil width") - self.add_output("Outer_width", 0.0, units="m", desc="Coil outer width") self.declare_partials("*", "*", method="fd") @@ -246,35 +278,30 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): q = discrete_inputs["q"] m = discrete_inputs["m"] D_sc = float(inputs["D_sc"]) - N_sc = int(inputs["N_sc"]) + N_sc = float(inputs["N_sc"]) # int I_sc = float(inputs["I_sc"]) - N_c = int(inputs["N_c"]) + N_c = float(inputs["N_c"]) # int delta_em = float(inputs["delta_em"]) I_s = float(inputs["I_s"]) + Slots = float(inputs["Slots"]) + Y_q = float(inputs["y_Q"]) radius_sc = D_sc / 2 tau_p = np.pi * (radius_sc * 2 + 2 * h_sc) / (2 * p1) theta_p_r = tau_p / (radius_sc + h_sc) theta_p_d = np.rad2deg(theta_p_r) + x1, x2, x3, x4, x5, x6, x7, x8 = inputs["field_coil_x"] + y1, y2, y3, y4, y5, y6, y7, y8 = inputs["field_coil_y"] + xlabel1, xlabel2 = inputs["field_coil_xlabel"] + ylabel1, ylabel2 = inputs["field_coil_ylabel"] # Build the geometry of the generator sector - if (alpha_d <= 0) or (beta_d > theta_p_d * 0.5): - outputs["B_g"] = 7 - outputs["B_coil_max"] = 12 - outputs["B_rymax"] = 5 - outputs["Torque_actual"] = 50e06 - outputs["Sigma_shear"] = 300000 - outputs["Sigma_normal"] = 200000 + if (alpha_d <= 0) or (float(inputs["con_angle"]) < 0): + outputs = bad_inputs(outputs) else: - femm.openfemm(1) - femm.newdocument(0) - femm.mi_probdef(0, "meters", "planar", 1.0e-8, l_s, 30) - slot_radius = D_a * 0.5 - h_s yoke_radius = slot_radius - h_yr h42 = D_a / 2 - 0.5 * h_s - Slots = 2 * q * m * p1 Slots_pp = q * m - Y_q = int(Slots / (2 * p1)) tau_s = np.pi * D_a / Slots # bs_taus = 0.45 #UNUSED @@ -286,104 +313,17 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): tau_p = np.pi * (radius_sc * 2 + 2 * h_sc) / (2 * p1) Current = 0 theta_p_d = np.rad2deg(theta_p_r) - alphap_taup_r = theta_p_r - 2 * beta_r - alphap_taup_angle_r = beta_r + alphap_taup_r - alpha_r - alphap_taup_angle_d = np.rad2deg(alphap_taup_angle_r) + # alphap_taup_r = theta_p_r - 2 * beta_r + # alphap_taup_angle_r = beta_r + alphap_taup_r - alpha_r + # alphap_taup_angle_d = np.rad2deg(alphap_taup_angle_r) # Draw the field coil - m_1 = -1 / (np.tan(theta_p_r * 0.5)) # slope of the tangent line - x_coord = radius_sc * np.cos(theta_p_r * 0.5) - y_coord = radius_sc * np.sin(theta_p_r * 0.5) - c1 = y_coord - m_1 * x_coord - - if alpha_r <= 1: - angle = alpha_r - else: - angle = np.tan(alpha_r) - m_2 = angle - c2 = 0 - - m_3 = m_1 # tangent offset - x_coord3 = (radius_sc + h_sc) * np.cos(theta_p_r * 0.5) - y_coord3 = (radius_sc + h_sc) * np.sin(theta_p_r * 0.5) - c3 = y_coord3 - m_3 * x_coord3 - - mlabel = m_1 - - x_label = (radius_sc + h_sc * 0.5) * np.cos(theta_p_r * 0.5) - y_label = (radius_sc + h_sc * 0.5) * np.sin(theta_p_r * 0.5) - clabel = y_label - mlabel * x_label - - mlabel2 = np.tan(alpha_r + (beta_r - alpha_r) * 0.5) - clabel2 = 0 - - mlabel3 = np.tan(theta_p_r - alpha_r - (beta_r - alpha_r) * 0.5) - clabel3 = 0 - - m_6 = np.tan(beta_r) - c6 = 0 - - x1 = (c1 - c2) / (m_2 - m_1) - y1 = m_1 * x1 + c1 + femm.openfemm(1) + femm.newdocument(0) + femm.mi_probdef(0, "meters", "planar", 1.0e-8, l_s, 30) femm.mi_addnode(x1, y1) - - m_4 = np.tan(theta_p_r * 0.5) - - c4 = y1 - m_4 * x1 - - x2 = (c3 - c4) / (m_4 - m_3) - y2 = m_4 * x2 + c4 - femm.mi_addnode(x2, y2) - - x4 = (c1 - c6) / (m_6 - m_1) - y4 = m_1 * x4 + c1 - - m_5 = np.tan(theta_p_r * 0.5) - c5 = y4 - m_5 * x4 - - x3 = (c3 - c5) / (m_5 - m_3) - y3 = m_3 * x3 + c3 - - m_7 = np.tan(theta_p_r - alpha_r) - c7 = 0 - - x5 = (c1 - c7) / (m_7 - m_1) - y5 = m_1 * x5 + c1 - - m_8 = m_4 - c8 = y5 - m_8 * x5 - - x6 = (c3 - c8) / (m_8 - m_3) - y6 = m_3 * x6 + c3 - - m_9 = np.tan(theta_p_r - beta_r) - c9 = 0 - - x8 = (c1 - c9) / (m_9 - m_1) - y8 = m_1 * x8 + c1 - - m_10 = m_5 - c10 = y8 - m_10 * x8 - - x7 = (c3 - c10) / (m_10 - m_3) - y7 = m_3 * x7 + c3 - - xlabel1 = (clabel - clabel2) / (mlabel2 - mlabel) - ylabel1 = mlabel * xlabel1 + clabel - - xlabel2 = (clabel - clabel3) / (mlabel3 - mlabel) - ylabel2 = mlabel * xlabel2 + clabel - - outputs["a_m"] = 2 * ( - np.sqrt( - (radius_sc * np.cos(theta_p_r * 0.5) - x4) ** 2 + ((radius_sc * np.sin(theta_p_r * 0.5) - y4) ** 2) - ) - ) - outputs["W_sc"] = np.sqrt((x1 - x4) ** 2 + (y1 - y4) ** 2) - outputs["Outer_width"] = outputs["a_m"] + 2 * outputs["W_sc"] - femm.mi_addnode(x3, y3) femm.mi_addnode(x4, y4) femm.mi_addnode(x5, y5) @@ -398,7 +338,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): femm.mi_addsegment(x8, y8, x7, y7) femm.mi_addsegment(x5, y5, x8, y8) femm.mi_addsegment(x6, y6, x7, y7) - femm.mi_addnode(0, 0) + # femm.mi_addnode(0, 0) # Draw the stator slots and Stator coils femm.mi_addnode(D_a / 2 * np.cos(0), D_a / 2 * np.sin(0)) @@ -565,7 +505,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): # femm.mi_addsegment(slot_radius*np.cos(theta_p_r),slot_radius*np.sin(theta_p_r),D_a/2*np.cos(theta_p_r),D_a/2*np.sin(theta_p_r)) - r_o = (radius_sc + h_sc) * 3 + r_o = (radius_sc + h_sc) * 2 # femm.mi_addsegment(D_a/2*np.cos(0),D_a/2*np.sin(0),radius_sc*np.cos(0),radius_sc*np.sin(0)) @@ -611,34 +551,36 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): yoke_radius * np.cos(0), yoke_radius * np.sin(0), slot_radius * np.cos(0), slot_radius * np.sin(0) ) femm.mi_addsegment(D_a * 0.5 * np.cos(0), D_a * 0.5 * np.sin(0), r_o * np.cos(0), r_o * np.sin(0)) - femm.mi_addsegment(0, 0, yoke_radius * np.cos(0), yoke_radius * np.sin(0)) - femm.mi_selectsegment(0.0, 0) - femm.mi_setsegmentprop("apbc1", 100, 0, 0, 6) - femm.mi_clearselected() + femm.mi_selectarcsegment(yoke_radius * np.cos(0), yoke_radius * np.sin(0)) + femm.mi_setarcsegmentprop(5, "Dirichlet", 0, 5) + + # femm.mi_selectsegment(0.0, 0) + # femm.mi_setsegmentprop("apbc1", 100, 0, 0, 6) + # femm.mi_clearselected() femm.mi_selectsegment(slot_radius * 0.999 * np.cos(0), slot_radius * 0.999 * np.sin(0)) - femm.mi_setsegmentprop("apbc2", 100, 0, 0, 6) + femm.mi_setsegmentprop("apbc1", 100, 0, 0, 6) femm.mi_clearselected() femm.mi_selectsegment(h42 * 0.99 * np.cos(0), h42 * 0.99 * np.sin(0)) - femm.mi_setsegmentprop("apbc3", 100, 0, 0, 6) + femm.mi_setsegmentprop("apbc2", 100, 0, 0, 6) femm.mi_clearselected() femm.mi_selectsegment(0.99 * D_a * 0.5 * np.cos(0), 0.99 * D_a * 0.5 * np.sin(0)) - femm.mi_setsegmentprop("apbc4", 1, 0, 0, 6) + femm.mi_setsegmentprop("apbc3", 1, 0, 0, 6) femm.mi_clearselected() femm.mi_selectsegment(r_o * 0.95 * np.cos(0), r_o * 0.95 * np.sin(0)) - femm.mi_setsegmentprop("apbc5", 100, 0, 0, 6) + femm.mi_setsegmentprop("apbc4", 100, 0, 0, 6) femm.mi_clearselected() femm.mi_selectgroup(6) femm.mi_copyrotate(0, 0, theta_p_d, 1) - femm.mi_selectsegment(h42 * 0.99 * np.cos(theta_p_r), h42 * 0.99 * np.sin(theta_p_r)) - femm.mi_setsegmentprop("apbc3", 100, 0, 0, 6) - femm.mi_clearselected() - femm.mi_selectsegment(D_a * 0.5 * 0.99 * np.cos(theta_p_r), D_a * 0.5 * 0.99 * np.sin(theta_p_r)) - femm.mi_setsegmentprop("apbc4", 100, 0, 0, 6) - femm.mi_clearselected() + # femm.mi_selectsegment(h42 * 0.99 * np.cos(theta_p_r), h42 * 0.99 * np.sin(theta_p_r)) + # femm.mi_setsegmentprop("apbc3", 100, 0, 0, 6) + # femm.mi_clearselected() + # femm.mi_selectsegment(D_a * 0.5 * 0.99 * np.cos(theta_p_r), D_a * 0.5 * 0.99 * np.sin(theta_p_r)) + # femm.mi_setsegmentprop("apbc4", 100, 0, 0, 6) + # femm.mi_clearselected() iron_label = yoke_radius + (slot_radius - yoke_radius) * 0.5 @@ -706,57 +648,63 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): # femm.mi_selectlabel(layer_2*np.cos(theta_b_s*0.25),layer_2*np.sin(theta_b_s*0.25)) # femm.mi_copyrotate(0,0,(theta_p_r-theta_b_t*0.5)*180/np.pi,1) - femm.mi_addblocklabel( - yoke_radius * 0.5 * np.cos(theta_p_r * 0.5), yoke_radius * 0.5 * np.sin(theta_p_r * 0.5) - ) - femm.mi_selectlabel( - yoke_radius * 0.5 * np.cos(theta_p_r * 0.5), yoke_radius * 0.5 * np.sin(theta_p_r * 0.5) - ) - femm.mi_setblockprop("Air", 1, 1, "incircuit", 0, 7, 0) - femm.mi_clearselected() - - femm.mi_addblocklabel(r_o * 0.75 * np.cos(theta_p_r * 0.5), r_o * 0.75 * np.sin(theta_p_r * 0.5)) - femm.mi_selectlabel(r_o * 0.75 * np.cos(theta_p_r * 0.5), r_o * 0.75 * np.sin(theta_p_r * 0.5)) + # femm.mi_addblocklabel( + # yoke_radius * 0.5 * np.cos(theta_p_r * 0.5), yoke_radius * 0.5 * np.sin(theta_p_r * 0.5) + # ) + # femm.mi_selectlabel( + # yoke_radius * 0.5 * np.cos(theta_p_r * 0.5), yoke_radius * 0.5 * np.sin(theta_p_r * 0.5) + # ) + # femm.mi_setblockprop("Air", 1, 1, "incircuit", 0, 7, 0) + # femm.mi_clearselected() + + femm.mi_addblocklabel(D_sc * 0.5 * np.cos(theta_p_r * 0.5), D_sc * 0.5 * np.sin(theta_p_r * 0.5)) + femm.mi_selectlabel(D_sc * 0.5 * np.cos(theta_p_r * 0.5), D_sc * 0.5 * np.sin(theta_p_r * 0.5)) femm.mi_setblockprop("Air", 1, 1, "incircuit", 0, 7, 0) femm.mi_clearselected() pitch = 1 Phases = ["A+", "D+", "C-", "F-", "B+", "E+", "A-", "D-", "C-"] - N_c_a = [2 * N_c, 4 * N_c, 4 * N_c, 4 * N_c, 4 * N_c, 4 * N_c, 2 * N_c, 2 * N_c, 2 * N_c] + # N_c_a = [2 * N_c, 4 * N_c, 4 * N_c, 4 * N_c, 4 * N_c, 4 * N_c, 2 * N_c, 2 * N_c, 2 * N_c] count = 0 angle_r = theta_b_t * 0.5 + theta_b_s * 0.5 delta_theta = theta_b_t + theta_b_s - for pitch in range(1, Y_q, 2): - femm.mi_selectlabel( - layer_2 * np.cos(angle_r + (pitch - 1) * (delta_theta)), - layer_2 * np.sin(angle_r + (pitch - 1) * (delta_theta)), - ) - femm.mi_selectlabel( - layer_2 * np.cos(angle_r + pitch * delta_theta), layer_2 * np.sin(angle_r + pitch * delta_theta) - ) - femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c) - # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c_a[count]) - femm.mi_clearselected() - count = count + 1 + for pitch in range(1, int(np.ceil(Y_q)), 2): + try: + femm.mi_selectlabel( + layer_2 * np.cos(angle_r + (pitch - 1) * (delta_theta)), + layer_2 * np.sin(angle_r + (pitch - 1) * (delta_theta)), + ) + femm.mi_selectlabel( + layer_2 * np.cos(angle_r + pitch * delta_theta), layer_2 * np.sin(angle_r + pitch * delta_theta) + ) + femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c) + # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count], 0, 8, N_c_a[count]) + femm.mi_clearselected() + count = count + 1 + except: + continue count = 0 angle_r = theta_b_t * 0.5 + theta_b_s * 0.5 delta_theta = theta_b_t + theta_b_s - for pitch in range(1, Y_q, 2): - femm.mi_selectlabel( - layer_1 * np.cos(angle_r + (pitch - 1) * (delta_theta)), - layer_1 * np.sin(angle_r + (pitch - 1) * (delta_theta)), - ) - femm.mi_selectlabel( - layer_1 * np.cos(angle_r + pitch * delta_theta), layer_1 * np.sin(angle_r + pitch * delta_theta) - ) - femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c) - # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c_a[count+1]) - femm.mi_clearselected() - count = count + 1 + for pitch in range(1, int(np.ceil(Y_q)), 2): + try: + femm.mi_selectlabel( + layer_1 * np.cos(angle_r + (pitch - 1) * (delta_theta)), + layer_1 * np.sin(angle_r + (pitch - 1) * (delta_theta)), + ) + femm.mi_selectlabel( + layer_1 * np.cos(angle_r + pitch * delta_theta), layer_1 * np.sin(angle_r + pitch * delta_theta) + ) + femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c) + # femm.mi_setblockprop("20 SWG", 1, 1, Phases[count + 1], 0, 8, N_c_a[count+1]) + femm.mi_clearselected() + count = count + 1 + except: + continue # Now, the finished input geometry can be displayed. femm.mi_zoomnatural() @@ -764,30 +712,37 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): femm.mi_saveas("coil_design_new.fem") # Analyze geometry with pyfemm - femm.mi_analyze() + # femm.mi_analyze() + try: - # Load and post-process results - femm.mi_loadsolution() - n = 0 - outputs["B_g"], outputs["B_rymax"], outputs["B_coil_max"], outputs["Sigma_normal"] = run_post_process( - D_a, radius_sc, h_sc, slot_radius, theta_p_r, alpha_r, beta_r, n - ) - Load_line_slope = I_sc / float(outputs["B_coil_max"]) - # print("Computing load line with slope {}".format(Load_line_slope)) - a = 5.8929 - b = -(Load_line_slope + 241.32) - c = 1859.9 - B_o = (-b - np.sqrt(b ** 2.0 - 4.0 * a * c)) / 2.0 / a - # I_c = B_o * Load_line_slope #UNUSED - - # Populate openmdao outputs - # Max current from manufacturer of superconducting coils, quadratic fit - outputs["margin_I_c"] = float( - 3.5357 * outputs["B_coil_max"] ** 2.0 - 144.79 * outputs["B_coil_max"] + 1116.0 - ) - outputs["Critical_current_ratio"] = I_sc / outputs["margin_I_c"] - # B_o is the max allowable flux density at the coil, B_coil_max is the max value from femm - outputs["Coil_max_ratio"] = outputs["B_coil_max"] / B_o - outputs["Torque_actual"], outputs["Sigma_shear"] = B_r_B_t( - D_a, l_s, p1, delta_em, theta_p_r, I_s, theta_b_t, theta_b_s, layer_1, layer_2, Y_q, N_c, tau_p - ) + femm.mi_analyze() + + # Load and post-process results + femm.mi_loadsolution() + n = 0 + outputs["B_g"], outputs["B_rymax"], B_coil_max, outputs["Sigma_normal"] = run_post_process( + D_a, radius_sc, h_sc, slot_radius, theta_p_r, alpha_r, beta_r, n + ) + outputs["B_coil_max"] = B_coil_max + Load_line_slope = I_sc / B_coil_max + # SHOULD STRONGLY CONSIDER A USER DEFINED LIMIT INSTEAD + a = 5.8929 + b = -(Load_line_slope + 241.32) + c = 1859.9 + B_o = (-b - np.sqrt(b ** 2.0 - 4.0 * a * c)) / 2.0 / a + # I_c = B_o * Load_line_slope #UNUSED + # Populate openmdao outputs + # Max current from manufacturer of superconducting coils, quadratic fit + # outputs["margin_I_c"] = 3.5357 * B_coil_max ** 2.0 - 144.79 * B_coil_max + 1116.0 + + outputs["margin_I_c"] = 3.5357 * B_o ** 2.0 - 144.79 * B_o + 1116.0 + + outputs["Critical_current_ratio"] = I_sc / outputs["margin_I_c"] + # B_o is the max allowable flux density at the coil, B_coil_max is the max value from femm + outputs["Coil_max_ratio"] = B_coil_max / B_o + outputs["Torque_actual"], outputs["Sigma_shear"] = B_r_B_t( + D_a, l_s, p1, delta_em, theta_p_r, I_s, theta_b_t, theta_b_s, layer_1, layer_2, Y_q, N_c, tau_p + ) + except: + + outputs = bad_inputs(outputs) diff --git a/lts/lts.py b/lts/lts.py index 268d8b3..eda102d 100644 --- a/lts/lts.py +++ b/lts/lts.py @@ -4,22 +4,56 @@ from lts.structural import LTS_Outer_Rotor_Structural +class LTS_Cost(om.ExplicitComponent): + def setup(self): + self.add_input("C_Cu", 0.0, units="USD/kg", desc="Specific cost of copper") + self.add_input("C_Fe", 0.0, units="USD/kg", desc="Specific cost of magnetic steel/iron") + self.add_input("C_Fes", 0.0, units="USD/kg", desc="Specific cost of structural steel") + self.add_input("C_NbTi", 0.0, units="USD/kg", desc="Specific cost of Magnet") + + # Mass of each material type + self.add_input("mass_copper", 0.0, units="kg", desc="Copper mass") + self.add_input("mass_iron", 0.0, units="kg", desc="Iron mass") + self.add_input("mass_SC", 0.0, units="kg", desc="Magnet mass") + self.add_input("mass_structural", 0.0, units="kg", desc="Structural mass") + + # Outputs + self.add_output("mass_total", 0.0, units="kg", desc="Structural mass") + self.add_output("cost_total", 0.0, units="USD", desc="Total cost") + self.declare_partials("*", "*", method="fd") + + def compute(self, inputs, outputs): + # Material cost as a function of material mass and specific cost of material + + outputs["mass_total"] = ( + inputs["mass_copper"] + inputs["mass_iron"] + inputs["mass_SC"] + inputs["mass_structural"] + ) + + outputs["cost_total"] = ( + inputs["mass_copper"] * inputs["C_Cu"] + + inputs["mass_iron"] * inputs["C_Fe"] + + inputs["mass_SC"] * inputs["C_NbTi"] + + inputs["mass_structural"] * inputs["C_Fes"] + ) + + class LTS_Outer_Rotor_Opt(om.Group): def initialize(self): self.options.declare("modeling_options") def setup(self): - self.linear_solver = lbgs = om.LinearBlockJac() # om.LinearBlockGS() - self.nonlinear_solver = nlbgs = om.NonlinearBlockGS() - nlbgs.options["maxiter"] = 3 - nlbgs.options["atol"] = 1e-2 - nlbgs.options["rtol"] = 1e-8 - nlbgs.options["iprint"] = 2 + # self.linear_solver = lbgs = om.LinearBlockJac() # om.LinearBlockGS() + # self.nonlinear_solver = nlbgs = om.NonlinearBlockGS() + # nlbgs.options["maxiter"] = 3 + # nlbgs.options["atol"] = 1e-2 + # nlbgs.options["rtol"] = 1e-8 + # nlbgs.options["iprint"] = 2 modeling_options = self.options["modeling_options"] ivcs = om.IndepVarComp() ivcs.add_output("P_rated", 0.0, units="W", desc="Rated Power") ivcs.add_output("T_rated", 0.0, units="N*m", desc="Torque") + ivcs.add_output("E_p_target", 0.0, units="V", desc="Target terminal voltage") ivcs.add_output("N_nom", 0.0, units="rpm", desc="rated speed") ivcs.add_output("D_a", 0.0, units="m", desc="Armature outer diameter") ivcs.add_output("delta_em", 0.0, units="m", desc="Field coil diameter") @@ -28,7 +62,7 @@ def setup(self): ivcs.add_output("h_yr", 0.0, units="m", desc="Rotor yoke height") ivcs.add_output("h_sc", 0.0, units="m", desc="SC coil height ") - ivcs.add_output("alpha_p", 0.0, desc="pole arc coefficient") + # ivcs.add_output("alpha_p", 0.0, desc="pole arc coefficient") ivcs.add_output("alpha", 0.0, units="deg", desc="Start angle of field coil") ivcs.add_output("dalpha", 0.0, units="deg", desc="Start angle of field coil") ivcs.add_output("I_sc", 0.0, units="A", desc="SC current") @@ -49,7 +83,12 @@ def setup(self): ivcs.add_output("rho_Fe", 0.0, units="kg/(m**3)", desc="Electrical Steel density ") ivcs.add_output("rho_Copper", 0.0, units="kg/(m**3)", desc="Copper density") ivcs.add_output("rho_NbTi", 0.0, units="kg/(m**3)", desc="SC conductor mass density ") - ivcs.add_output("rho_Cu", 0.0, units="ohm*m", desc="Copper resistivity ") + ivcs.add_output("resisitivty_Cu", 0.0, units="ohm*m", desc="Copper resistivity ") + ivcs.add_output("C_Cu", 0.0, units="USD/kg", desc="Specific cost of copper") + ivcs.add_output("C_Fe", 0.0, units="USD/kg", desc="Specific cost of magnetic steel/iron") + ivcs.add_output("C_Fes", 0.0, units="USD/kg", desc="Specific cost of structural steel") + ivcs.add_output("C_NbTi", 0.0, units="USD/kg", desc="Specific cost of Magnet") + ivcs.add_output("U_b", 0.0, units="V", desc="brush voltage ") # ivcs.add_output("r_strand", 0.0, units="mm", desc="radius of the SC wire strand") # ivcs.add_output("k_pf_sc", 0.0, units="mm", desc="packing factor for SC wires") @@ -60,5 +99,5 @@ def setup(self): self.add_subsystem("geom", FEMM_Geometry(modeling_options=modeling_options), promotes=["*"]) self.add_subsystem("results", md.Results(), promotes=["*"]) self.add_subsystem("struct", LTS_Outer_Rotor_Structural(), promotes=["*"]) - + self.add_subsystem("cost", LTS_Cost(), promotes=["*"]) self.connect("Torque_actual", "T_e") diff --git a/lts/magnetics_design.py b/lts/magnetics_design.py index 661fb41..f0ac087 100644 --- a/lts/magnetics_design.py +++ b/lts/magnetics_design.py @@ -3,7 +3,7 @@ Copyright (c) NREL. All rights reserved. Electromagnetic design based on conventional magnetic circuit laws Structural design based on {Structural mass in direct-drive permanent magnet electrical generators by -McDonald,A.S. et al. IET Renewable Power Generation(2008),2(1):3 http://dx.doi.org/10.1049/iet-rpg:20070071 +McDonald,A.S. et al. IET Renewable Power Generation(2008),2(1):3 http://dx.doi.org/10.1049/iet-rpg:20070071 """ import numpy as np @@ -31,13 +31,12 @@ def setup(self): # field coil parameters self.add_input("h_sc", 0.0, units="m", desc="SC coil height") - self.add_input("alpha_p", 0.0, desc="pole arc coefficient") + # self.add_input("alpha_p", 0.0, desc="pole arc coefficient") self.add_input("alpha", 0.0, units="deg", desc="Start angle of field coil") self.add_input("dalpha", 0.0, units="deg", desc="Angle subtended by field coil") self.add_input("h_yr", 0.0, units="m", desc="rotor yoke height") # self.add_input("I_sc", 0.0, units="A", desc="SC current ") self.add_input("N_sc", 0.0, desc="Number of turns of SC field coil") - self.add_input("W_sc", 0.0, units="m", desc="SC coil width") self.add_input("N_c", 0.0, desc="Number of turns per coil") self.add_input("p", 0.0, desc="Pole pairs ") self.add_output("p1", 0.0, desc="Pole pairs ") @@ -48,35 +47,41 @@ def setup(self): self.add_output("N_s", 0.0, desc="Number of turns per phase in series") self.add_output("N_l", 0.0, desc="Number of layers of the SC field coil") # self.add_output("Dia_sc", 0.0, units="m", desc="field coil diameter") - self.add_input("Outer_width", 0.0, units="m", desc="Coil outer width") - self.add_input("a_m", 0.0, units="m", desc="Coil separation distance") + # self.add_input("Outer_width", 0.0, units="m", desc="Coil outer width") self.add_output("beta", 0.0, units="deg", desc="End angle of field coil") - self.add_output("con_angle", 0.0, units="deg", desc="End angle of field coil") - self.add_output("con_angle2", 0.0, units="deg", desc="End angle of field coil") + self.add_output("con_angle", 0.0, units="deg", desc="Geometric constraint for valid setup") + self.add_output("con_angle2", 0.0, units="deg", desc="Geometric constraint for valid setup") self.add_output("theta_p", 0.0, units="deg", desc="Pole pitch angle in degrees") # Material properties self.add_input("rho_Fe", 0.0, units="kg/(m**3)", desc="Electrical Steel density ") self.add_input("rho_Copper", 0.0, units="kg/(m**3)", desc="Copper density") self.add_input("rho_NbTi", 0.0, units="kg/(m**3)", desc="SC conductor mass density ") - self.add_input("rho_Cu", 0.0, units="ohm*m", desc="Copper resistivity ") + self.add_input("resisitivty_Cu", 0.0, units="ohm*m", desc="Copper resistivity ") self.add_input("U_b", 0.0, units="V", desc="brush voltage ") self.add_input("P_rated", units="W", desc="Machine rating") self.add_input("N_nom", 0.0, units="rpm", desc="rated speed") - self.add_input("T_rated", 0.0, units="N*m", desc="Rated torque ") # self.add_input("r_strand", 0.0, units="mm", desc="radius of the SC wire strand") # self.add_input("k_pf_sc", 0.0, units="mm", desc="packing factor for SC wires") self.add_input("J_s", 0.0, units="A/(mm*mm)", desc="Stator winding current density") # self.add_input("J_c", 0.0, units="A/(mm*mm)", desc="SC critical current density") + # Field coil geometry + self.add_output("a_m", 0.0, units="m", desc="Coil separation distance") + self.add_output("W_sc", 0.0, units="m", desc="SC coil width") + self.add_output("Outer_width", 0.0, units="m", desc="Coil outer width") + self.add_output("field_coil_x", np.zeros(8), desc="Field coil points") + self.add_output("field_coil_y", np.zeros(8), desc="Field coil points") + self.add_output("field_coil_xlabel", np.zeros(2), desc="Field coil label points") + self.add_output("field_coil_ylabel", np.zeros(2), desc="Field coil label points") + # Magnetic loading self.add_output("tau_p", 0.0, units="m", desc="Pole pitch") - self.add_output("b_p", 0.0, units="m", desc="distance between positive and negative side of field coil") + # self.add_output("b_p", 0.0, units="m", desc="distance between positive and negative side of field coil") self.add_output("alpha_u", 0.0, units="rad", desc="slot angle") self.add_output("tau_v", 0.0, units="m", desc="Phase zone span") self.add_output("zones", 0.0, desc="Phase zones") - self.add_output("y_Q", 0.0, desc="Slots per pole also pole pitch") self.add_output("delta", 0.0, units="rad", desc="short-pitch angle") self.add_output("k_p1", 0.0, desc="Pitch factor-fundamental harmonic") self.add_output("k_d1", 0.0, desc="Distribution factor-fundamental harmonic") @@ -109,7 +114,6 @@ def setup(self): # self.add_output("Torque_constraint", 0.0, units="N/(m*m)", desc="Shear stress contraint") # Objective functions - self.add_output("Mass", 0.0, units="kg", desc="Actual mass") self.add_output("K_rad", desc="Aspect ratio") self.add_output("Cu_losses", units="W", desc="Copper losses") self.add_output("P_add", units="W", desc="Additional losses") @@ -117,157 +121,252 @@ def setup(self): # Other parameters # self.add_output("R_out", 0.0, units="m", desc="Outer radius") - self.add_output("S", 0.0, desc="Stator slots") + self.add_output("Slots", 0.0, desc="Stator slots") self.add_output("Slot_aspect_ratio", 0.0, desc="Slot aspect ratio") + self.add_output("y_Q", 0.0, desc="Slots per pole also pole pitch") # Mass Outputs - self.add_output("mass_SC", 0.0, units="kg", desc="SC conductor mass per racetrack") - self.add_output("Total_mass_SC", 0.0, units="kg", desc=" Total SC conductor mass") - self.add_output("Copper", 0.0, units="kg", desc="Copper Mass") - self.add_output("Iron", 0.0, units="kg", desc="Electrical Steel Mass") + self.add_output("mass_SC_racetrack", 0.0, units="kg", desc="SC conductor mass per racetrack") + self.add_output("mass_SC", 0.0, units="kg", desc=" Total SC conductor mass") + self.add_output("mass_copper", 0.0, units="kg", desc="Copper Mass") + self.add_output("mass_iron", 0.0, units="kg", desc="Electrical Steel Mass") + self.add_output("mass_active", 0.0, units="kg", desc="Actual mass") self.declare_partials("*", "*", method="fd") def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): # Unpack inputs - b_s_tau_s = inputs["b_s_tau_s"] - conductor_area = inputs["conductor_area"] + b_s_tau_s = float(inputs["b_s_tau_s"]) + conductor_area = float(inputs["conductor_area"]) + m = float(discrete_inputs["m"]) + q = float(discrete_inputs["q"]) + D_a = float(inputs["D_a"]) + l_s = float(inputs["l_s"]) + h_s = float(inputs["h_s"]) + h_sc = float(inputs["h_sc"]) + # alpha_p = float(inputs["alpha_p"]) + alpha_d = float(inputs["alpha"]) + alpha_r = np.deg2rad(alpha_d) + dalpha = float(inputs["dalpha"]) + h_yr = float(inputs["h_yr"]) + N_sc = float(inputs["N_sc"]) + N_c = float(inputs["N_c"]) + p = float(inputs["p"]) + delta_em = float(inputs["delta_em"]) + Y = float(inputs["Y"]) + I_s = float(inputs["I_s"]) + rho_Fe = float(inputs["rho_Fe"]) + rho_Copper = float(inputs["rho_Copper"]) + rho_NbTi = float(inputs["rho_NbTi"]) + resisitivty_Cu = float(inputs["resisitivty_Cu"]) + U_b = float(inputs["U_b"]) + P_rated = float(inputs["P_rated"]) + N_nom = float(inputs["N_nom"]) + J_s = float(inputs["J_s"]) ###################################################### Electromagnetic design############################################# - outputs["K_rad"] = inputs["l_s"][0] / (inputs["D_a"][0]) # Aspect ratio - outputs["D_sc"] = inputs["D_a"] + 2 * inputs["delta_em"] - outputs["R_sc"] = outputs["D_sc"] * 0.5 - outputs["p1"] = np.round(inputs["p"]) + outputs["K_rad"] = l_s / (D_a) # Aspect ratio + outputs["D_sc"] = D_sc = D_a + 2 * delta_em + outputs["R_sc"] = R_sc = D_sc * 0.5 + outputs["p1"] = p1 = p # np.round(p) # Calculating pole pitch - # r_s = 0.5 * inputs["D_a"] # Stator outer radius # UNUSED - outputs["tau_p"] = (np.pi * (inputs["D_a"] + 2 * inputs["delta_em"] + 2 * inputs["h_sc"])) / ( - 2 * outputs["p1"] - ) # pole pitch - outputs["b_p"] = inputs["alpha_p"] * outputs["tau_p"] + # r_s = 0.5 * D_a # Stator outer radius # UNUSED + outputs["tau_p"] = tau_p = np.pi * (R_sc + h_sc) / p1 + # outputs["b_p"] = alpha_p * tau_p # Calculating winding factor - m = discrete_inputs["m"] - q = discrete_inputs["q"] - outputs["S"] = q * 2 * outputs["p1"] * m - outputs["tau_s"] = np.pi * (inputs["D_a"]) / outputs["S"] # Slot pitch - outputs["alpha_u"] = outputs["p1"] * 2 * np.pi / outputs["S"] # slot angle - outputs["tau_v"] = outputs["tau_p"] / m - outputs["zones"] = 2 * outputs["p1"] * m - outputs["y_Q"] = outputs["S"] / (2 * outputs["p1"]) # coil span - outputs["delta"] = inputs["Y"] * np.pi * 0.5 / outputs["y_Q"] # short pitch angle #coil span - outputs["k_p1"] = np.sin(np.pi * 0.5 * inputs["Y"] / outputs["y_Q"]) - outputs["k_d1"] = np.sin(q * outputs["alpha_u"] * 0.5) / (q * np.sin(outputs["alpha_u"] * 0.5)) - outputs["k_w1"] = outputs["k_p1"] * outputs["k_d1"] + outputs["Slots"] = S = q * 2 * p1 * m + outputs["tau_s"] = tau_s = np.pi * (D_a) / S # Slot pitch + outputs["alpha_u"] = alpha_u = p1 * 2 * np.pi / S # slot angle + outputs["tau_v"] = tau_p / m + outputs["zones"] = 2 * p1 * m + outputs["y_Q"] = y_Q = S / (2 * p1) # coil span + outputs["delta"] = Y * np.pi * 0.5 / y_Q # short pitch angle #coil span + outputs["k_p1"] = k_p1 = np.sin(np.pi * 0.5 * Y / y_Q) + outputs["k_d1"] = k_d1 = np.sin(q * alpha_u * 0.5) / (q * np.sin(alpha_u * 0.5)) + outputs["k_w1"] = k_p1 * k_d1 # magnet width # alpha_p = np.pi / 2 * ratio # COMMENTING OUT BECAUSE ALPHA_P IS AN INPUT - outputs["b_s"] = b_s_tau_s * outputs["tau_s"] - outputs["b_t"] = outputs["tau_s"] - outputs["b_s"] # slot width + outputs["b_s"] = b_s = b_s_tau_s * tau_s + outputs["b_t"] = tau_s - b_s # slot width # UNUSED # gamma = 2 / np.pi * ( - # np.atan(outputs["b_s"] * 0.5 / inputs["delta_em"]) - # - (2 * inputs["delta_em"] / outputs["b_s"]) - # * log(((1 + (outputs["b_s"] * 0.5 / inputs["delta_em"]) ** 2)) ** 0.5) + # np.atan(b_s * 0.5 / delta_em) + # - (2 * delta_em / b_s) + # * log(((1 + (b_s * 0.5 / delta_em) ** 2)) ** 0.5) # ) - # k_C = outputs["tau_s"] / (outputs["tau_s"] - gamma * (outputs["b_s"])) # carter coefficient UNUSED - # g_eff = k_C * inputs["delta_em"] # UNUSED + # k_C = tau_s / (tau_s - gamma * (b_s)) # carter coefficient UNUSED + # g_eff = k_C * delta_em # UNUSED # angular frequency in radians - om_m = 2 * np.pi * inputs["N_nom"] / 60 - om_e = outputs["p1"] * om_m + om_m = 2 * np.pi * N_nom / 60 + om_e = p1 * om_m outputs["f"] = om_e / 2 / np.pi # outout frequency - # outputs['N_s'] = outputs['p1']*inputs['Slot_pole']*N_c*3 #2*m*p*q - cos_theta_end = 1 - (outputs["b_s"] / (outputs["b_t"] + outputs["b_s"]) ** 2) ** 0.5 - # l_end = 4 * outputs["tau_s"] * 0.5 / cos_theta_end - l_end = 10 * outputs["tau_s"] / 2 * np.tan(30 * np.pi / 180) + # outputs["N_s N_s = = p1*Slot_pole*N_c*3 #2*m*p*q + # cos_theta_end = 1 - (b_s / (b_t + b_s) ** 2) ** 0.5 + # l_end = 4 * tau_s * 0.5 / cos_theta_end + l_end = 10 * tau_s / 2 * np.tan(30 * np.pi / 180) - outputs["l_eff_rotor"] = inputs["l_s"] + 2 * l_end - # l_end = np.sqrt(3)/6*(10*outputs['tau_s']+) + outputs["l_eff_rotor"] = l_s + 2 * l_end + # l_end = np.sqrt(3)/6*(10*tau_s+) # Stator winding length ,cross-section and resistance - outputs["l_Cus"] = 8 * l_end + 2 * inputs["l_s"] # length of a turn - z = outputs["S"] # Number of coils - # A_slot = inputs["h_s"] * outputs["b_s"] # UNUSED - # d_cu = 2 * np.sqrt(outputs["A_Cuscalc"] / pi) # UNUSED - outputs["A_Cuscalc"] = inputs["I_s"] * 1e-06 / (inputs["J_s"]) - # k_fill = A_slot / (2 * inputs["N_c"] * outputs["A_Cuscalc"]) # UNUSED - outputs["N_s"] = int(inputs["N_c"]) * z / (m) # turns per phase - - outputs["R_s"] = ( - inputs["rho_Cu"] - * (1 + 20 * 0.00393) - * outputs["l_Cus"] - * outputs["N_s"] - * inputs["J_s"] - * 1000000 - / (2 * (inputs["I_s"])) - ) - # print ("Resitance per phase:" ,outputs["R_s"]) + outputs["l_Cus"] = l_Cus = 8 * l_end + 2 * l_s # length of a turn + z = S # Number of coils + # A_slot = h_s * b_s # UNUSED + # d_cu = 2 * np.sqrt(A_Cuscalc / pi) # UNUSED + outputs["A_Cuscalc"] = A_Cuscalc = I_s * 1e-06 / (J_s) + # k_fill = A_slot / (2 * N_c * A_Cuscalc) # UNUSED + outputs["N_s"] = N_s = N_c * z / (m) # turns per phase int(N_c) + + outputs["R_s"] = R_s = resisitivty_Cu * (1 + 20 * 0.00393) * l_Cus * N_s * J_s * 1e6 / (I_s) + # print ("Resitance per phase:" ,R_s) # r_strand =0.425e-3 - # inputs['N_sc'] = inputs['k_pf_sc']*outputs['W_sc']*inputs['h_sc']/pi*r_strand**2 - outputs["theta_p"] = np.rad2deg(outputs["tau_p"] / (outputs["R_sc"] + inputs["h_sc"])) - outputs["beta"] = inputs["alpha"] + inputs["dalpha"] - outputs["con_angle"] = 0.5 * outputs["theta_p"] - outputs["beta"] - outputs["con_angle2"] = outputs["beta"] - inputs["alpha"] - # outputs['beta'] = (outputs['theta_p'] - 2*inputs['alpha'])*0.5-2 - # random_degree = np.random.uniform(1.0, outputs["theta_p"] * 0.5 - inputs["alpha"] - 0.25) - # random_degree = np.mean([1.0, float(outputs["theta_p"]) * 0.5 - float(inputs["alpha"]) - 0.25]) - # outputs["beta"] = outputs["theta_p"] * 0.5 - random_degree + theta_p_r = tau_p / (R_sc + h_sc) + outputs["theta_p"] = theta_p_d = np.rad2deg(theta_p_r) + outputs["beta"] = beta_d = alpha_d + dalpha + outputs["con_angle"] = 0.5 * theta_p_d - beta_d + outputs["con_angle2"] = beta_d - alpha_d + beta_r = np.deg2rad(beta_d) + # outputs["beta beta = = (theta_p - 2*alpha)*0.5-2 + # random_degree = np.random.uniform(1.0, theta_p * 0.5 - alpha - 0.25) + # random_degree = np.mean([1.0, float(theta_p) * 0.5 - float(alpha) - 0.25]) + # outputs["beta"] = beta = theta_p * 0.5 - random_degree + + # Field coil geometry prep for FEMM + m_1 = -1 / (np.tan(theta_p_r * 0.5)) # slope of the tangent line + x_coord = R_sc * np.cos(theta_p_r * 0.5) + y_coord = R_sc * np.sin(theta_p_r * 0.5) + c1 = y_coord - m_1 * x_coord + + if alpha_r <= 1: + angle = alpha_r + else: + angle = np.tan(alpha_r) + m_2 = angle + c2 = 0 + + m_3 = m_1 # tangent offset + x_coord3 = (R_sc + h_sc) * np.cos(theta_p_r * 0.5) + y_coord3 = (R_sc + h_sc) * np.sin(theta_p_r * 0.5) + c3 = y_coord3 - m_3 * x_coord3 + + mlabel = m_1 + + x_label = (R_sc + h_sc * 0.5) * np.cos(theta_p_r * 0.5) + y_label = (R_sc + h_sc * 0.5) * np.sin(theta_p_r * 0.5) + clabel = y_label - mlabel * x_label + + mlabel2 = np.tan(alpha_r + (beta_r - alpha_r) * 0.5) + clabel2 = 0 + + mlabel3 = np.tan(theta_p_r - alpha_r - (beta_r - alpha_r) * 0.5) + clabel3 = 0 + + m_6 = np.tan(beta_r) + c6 = 0 + + x1 = (c1 - c2) / (m_2 - m_1) + y1 = m_1 * x1 + c1 + + m_4 = np.tan(theta_p_r * 0.5) + + c4 = y1 - m_4 * x1 + + x2 = (c3 - c4) / (m_4 - m_3) + y2 = m_4 * x2 + c4 + + x4 = (c1 - c6) / (m_6 - m_1) + y4 = m_1 * x4 + c1 + + m_5 = np.tan(theta_p_r * 0.5) + c5 = y4 - m_5 * x4 + + x3 = (c3 - c5) / (m_5 - m_3) + y3 = m_3 * x3 + c3 - outputs["l_sc"] = inputs["N_sc"] * (2 * inputs["l_s"] + np.pi * (inputs["a_m"] + inputs["W_sc"])) + m_7 = np.tan(theta_p_r - alpha_r) + c7 = 0 - outputs["l_eff_stator"] = inputs["l_s"] + (inputs["a_m"] + inputs["W_sc"]) + x5 = (c1 - c7) / (m_7 - m_1) + y5 = m_1 * x5 + c1 - # outputs['A_Cuscalc'] = I_n/J_s - # A_slot = 2*N_c*outputs['A_Cuscalc']*(10**-6)/k_sfil - outputs["Slot_aspect_ratio"] = inputs["h_s"] / outputs["b_s"] + m_8 = m_4 + c8 = y5 - m_8 * x5 + + x6 = (c3 - c8) / (m_8 - m_3) + y6 = m_3 * x6 + c3 + + m_9 = np.tan(theta_p_r - beta_r) + c9 = 0 + + x8 = (c1 - c9) / (m_9 - m_1) + y8 = m_1 * x8 + c1 + + m_10 = m_5 + c10 = y8 - m_10 * x8 + + x7 = (c3 - c10) / (m_10 - m_3) + y7 = m_3 * x7 + c3 + + xlabel1 = (clabel - clabel2) / (mlabel2 - mlabel) + ylabel1 = mlabel * xlabel1 + clabel + + xlabel2 = (clabel - clabel3) / (mlabel3 - mlabel) + ylabel2 = mlabel * xlabel2 + clabel + + outputs["a_m"] = a_m = 2 * ( + np.sqrt((R_sc * np.cos(theta_p_r * 0.5) - x4) ** 2 + ((R_sc * np.sin(theta_p_r * 0.5) - y4) ** 2)) + ) + outputs["W_sc"] = W_sc = np.sqrt((x1 - x4) ** 2 + (y1 - y4) ** 2) + outputs["Outer_width"] = a_m + 2 * W_sc + outputs["field_coil_x"] = np.r_[x1, x2, x3, x4, x5, x6, x7, x8] + outputs["field_coil_y"] = np.r_[y1, y2, y3, y4, y5, y6, y7, y8] + outputs["field_coil_xlabel"] = np.r_[xlabel1, xlabel2] + outputs["field_coil_ylabel"] = np.r_[ylabel1, ylabel2] + + # N_sc = k_pf_sc*W_sc*h_sc/np.pi*r_strand**2 + outputs["l_sc"] = l_sc = N_sc * (2 * l_s + np.pi * (a_m + W_sc)) + + outputs["l_eff_stator"] = l_s + (a_m + W_sc) + + # outputs["A_Cuscalc = A_Cuscalc = I_n/J_s + # A_slot = 2*N_c*A_Cuscalc*(10**-6)/k_sfil + outputs["Slot_aspect_ratio"] = h_s / b_s # Calculating stator current and electrical loading - # +I_s = sqrt(Z**2+(((outputs['E_p']-G**0.5)/(om_e*outputs['L_s'])**2)**2)) + # +I_s = sqrt(Z**2+(((E_p-G**0.5)/(om_e*L_s)**2)**2)) # Calculating volumes and masses - # V_Cus = m*L_Cus*(outputs['A_Cuscalc']*(10**-6)) # copper volume - # outputs['h_t'] = (inputs['h_s']+h_1+h_0) - V_Fery = ( - inputs["l_s"] - * np.pi - * 0.25 - * ((inputs["D_a"] - 2 * inputs["h_s"]) ** 2 - (inputs["D_a"] - 2 * inputs["h_s"] - 2 * inputs["h_yr"]) ** 2) - ) # volume of iron in stator tooth - # outputs['Copper'] = V_Cus*inputs['rho_Copper'] - outputs["Iron"] = V_Fery * inputs["rho_Fe"] # Mass of stator yoke - # k_pf = 0.8 # UNUSED - # inputs['N_sc'] =int(inputs['N_sc']) - # outputs['Dia_sc'] =2*((k_pf*outputs['W_sc']*inputs['h_sc'])/inputs['N_sc']/pi)**0.5 - # k_pf = inputs["N_sc"] * (1.8e-3 * 1.2e-3) / (outputs["W_sc"] * inputs["h_sc"]) # UNUSED - # outputs['SC mass'] = - - # Calculating Losses - ##1. Copper Losses - # outputs['N_s'] =int(outputs['N_s']) - # outputs['N_s'] =q*inputs['p'] - # print (inputs['alpha'],outputs['beta']) - - outputs["N_l"] = int(inputs["h_sc"] / (1.2e-3)) # round later! + # V_Cus = m*L_Cus*(A_Cuscalc*(10**-6)) # copper volume + # outputs["h_t h_t = = (h_s+h_1+h_0) + # volume of iron in stator tooth + V_Fery = 0.25 * np.pi * l_s * ((D_a - 2 * h_s) ** 2 - (D_a - 2 * h_s - 2 * h_yr) ** 2) + # outputs["Copper Copper = = V_Cus*rho_Copper + outputs["mass_iron"] = Iron = V_Fery * rho_Fe # Mass of stator yoke + + outputs["N_l"] = h_sc / (1.2e-3) # round later! # 0.01147612156295224312590448625181 - outputs["mass_SC"] = outputs["l_sc"] * conductor_area * inputs["rho_NbTi"] - outputs["Total_mass_SC"] = outputs["p1"] * outputs["mass_SC"] - V_Cus = m * outputs["l_Cus"] * outputs["N_s"] * (outputs["A_Cuscalc"]) - outputs["Copper"] = V_Cus * inputs["rho_Copper"] + outputs["mass_SC_racetrack"] = mass_SC = l_sc * conductor_area * rho_NbTi + outputs["mass_SC"] = Total_mass_SC = p1 * mass_SC + V_Cus = m * l_Cus * N_s * (A_Cuscalc) + outputs["mass_copper"] = Copper = V_Cus * rho_Copper - outputs["Mass"] = outputs["Total_mass_SC"] + outputs["Iron"] + outputs["Copper"] - outputs["A_1"] = (2 * inputs["I_s"] * outputs["N_s"] * m) / (np.pi * (inputs["D_a"])) - outputs["Cu_losses"] = m * (inputs["I_s"] * 0.707) ** 2 * outputs["R_s"] - outputs["P_add"] = 0.01 * inputs["P_rated"] - outputs["P_brushes"] = 2 * inputs["U_b"] * inputs["I_s"] + outputs["mass_active"] = Total_mass_SC + Iron + Copper + outputs["A_1"] = (2 * I_s * N_s * m) / (np.pi * (D_a)) - # print (inputs["N_sc"],inputs["I_s"], outputs["p1"], inputs["D_a"],inputs["delta_em"], inputs["N_c"],outputs["S"]) + outputs["Cu_losses"] = m * (I_s * 0.707) ** 2 * R_s + outputs["P_add"] = 0.01 * P_rated + outputs["P_brushes"] = 6 * U_b * I_s * 0.707 - # print(outputs["mass_SC"]) + # print (N_sc,I_s, p1, D_a,delta_em, N_c,S) + + # print(mass_SC) class Results(om.ExplicitComponent): @@ -278,19 +377,22 @@ def setup(self): self.add_input("I_sc", 0.0, units="A", desc="SC current ") self.add_input("N_sc", 0.0, desc="Number of turns of SC field coil") self.add_input("N_l", 0.0, desc="Number of layers of the SC field coil") - self.add_input("D_sc", 0.0, units="m", desc="field coil diameter ") + self.add_input("D_a", 0.0, units="m", desc="Armature diameter ") self.add_input("k_w1", 0.0, desc="Winding factor- fundamental harmonic") self.add_input("B_rymax", 0.0, desc="Peak Rotor yoke flux density") self.add_input("B_g", 0.0, desc="Peak air gap flux density ") self.add_input("N_s", 0.0, desc="Number of turns per phase in series") self.add_input("N_nom", 0.0, units="rpm", desc="rated speed") self.add_input("p1", 0.0, desc="Pole pairs ") - self.add_input("Iron", 0.0, units="kg", desc="Electrical Steel Mass") + self.add_input("mass_iron", 0.0, units="kg", desc="Electrical Steel Mass") + self.add_input("T_rated", 0.0, units="N*m", desc="Rated torque ") + self.add_input("Torque_actual", 0.0, units="N*m", desc="Shear stress actual") self.add_input("P_rated", units="W", desc="Machine rating") self.add_input("Cu_losses", units="W", desc="Copper losses") self.add_input("P_add", units="W", desc="Additional losses") self.add_input("P_brushes", units="W", desc="brush losses") self.add_input("l_s", 0.0, units="m", desc="Stator core length") + self.add_input("E_p_target", 0.0, units="V", desc="target terminal voltage") # self.add_output("con_I_sc", 0.0, units="A/(mm*mm)", desc="SC current ") # self.add_output("con_N_sc", 0.0, desc="Number of turns of SC field coil") self.add_output("E_p", 0.0, units="V", desc="terminal voltage") @@ -298,6 +400,8 @@ def setup(self): self.add_output("P_Fe", units="W", desc="Iron losses") self.add_output("P_Losses", units="W", desc="Total power losses") self.add_output("gen_eff", desc="Generator efficiency") + self.add_output("torque_ratio", desc="Whether torque meets requirement") + self.add_output("E_p_ratio", desc="Whether terminal voltage meets requirement") self.declare_partials("*", "*", method="fd") @@ -305,23 +409,42 @@ def compute(self, inputs, outputs): # Unpack inputs K_h = inputs["K_h"] K_e = inputs["K_e"] - - # outputs["con_N_sc"] = inputs["N_sc"] - inputs["N_sc_out"] - # outputs["con_I_sc"] = inputs["I_sc"] - inputs["I_sc_out"] - outputs["N_sc_layer"] = int(inputs["N_sc"] / inputs["N_l"]) + N_sc = float(inputs["N_sc"]) + N_l = float(inputs["N_l"]) + D_a = float(inputs["D_a"]) + k_w1 = float(inputs["k_w1"]) + B_rymax = float(inputs["B_rymax"]) + B_g = float(inputs["B_g"]) + N_s = float(inputs["N_s"]) + N_nom = float(inputs["N_nom"]) + p1 = float(inputs["p1"]) + Iron = float(inputs["mass_iron"]) + P_rated = float(inputs["P_rated"]) + Cu_losses = float(inputs["Cu_losses"]) + P_add = float(inputs["P_add"]) + P_brushes = float(inputs["P_brushes"]) + l_s = float(inputs["l_s"]) + T_rated = float(inputs["T_rated"]) + T_actual = float(inputs["Torque_actual"]) + E_p_target = float(inputs["E_p_target"]) + # outputs["con_N_sc"] = con_N_sc = N_sc - N_sc_out + # outputs["con_I_sc"] = con_I_sc = I_sc - I_sc_out + outputs["N_sc_layer"] = int(N_sc / N_l) # Calculating voltage per phase - om_m = 2 * np.pi * inputs["N_nom"] / 60 - outputs["E_p"] = inputs["l_s"] * (inputs["D_sc"] * 0.5 * inputs["k_w1"] * inputs["B_g"] * om_m * inputs["N_s"]) + om_m = 2 * np.pi * N_nom / 60 + outputs["E_p"] = E_p = l_s * (D_a * 0.5 * k_w1 * B_g * om_m * N_s) * 3 ** 0.5 * (1 / 2 ** 0.5) * 1.12253 - # print ("Voltage and lengths are:",outputs["E_p"],inputs["l_s"] ) + # print ("Voltage and lengths are:",outputs["E_p,l_s ) - om_e = inputs["p1"] * om_m - outputs["P_Fe"] = ( - 2 - * (inputs["B_rymax"] / 1.5) ** 2 - * (K_h * (om_e * 0.5 / np.pi) + K_e * (om_e * 0.5 / 50) ** 2) - * inputs["Iron"] - ) - outputs["P_Losses"] = inputs["Cu_losses"] + outputs["P_Fe"] + inputs["P_add"] + inputs["P_brushes"] - outputs["gen_eff"] = 1 - outputs["P_Losses"] / inputs["P_rated"] + f_e = 2 * p1 * N_nom / 120 + outputs["P_Fe"] = P_Fe = ( + 2 * K_h * (f_e / 60) * (B_rymax / 1.5) ** 2 + 2 * K_e * (f_e / 60) ** 2 * (B_rymax / 1.5) ** 2 + ) * Iron + + outputs["P_Losses"] = P_Losses = Cu_losses + P_Fe + P_add + P_brushes + outputs["gen_eff"] = 1 - P_Losses / P_rated + outputs["torque_ratio"] = T_actual / T_rated + outputs["E_p_ratio"] = E_p / E_p_target + + # print(outputs["gen_eff"],outputs["torque_ratio"],outputs["E_p_ratio"]) diff --git a/lts/structural.py b/lts/structural.py index b88e8f1..ecced7f 100644 --- a/lts/structural.py +++ b/lts/structural.py @@ -27,18 +27,18 @@ def shell_constant(R, t, l, x, v, E): def plate_constant(a, b, v, r_o, t, E): D = E * t ** 3 / (12 * (1 - v ** 2)) - C_2 = 0.25 * (1 - (b / a) ** 2 * (1 + 2 * np.log10(a / b))) - C_3 = 0.25 * (b / a) * (((b / a) ** 2 + 1) * np.log10(a / b) + (b / a) ** 2 - 1) + C_2 = 0.25 * (1 - (b / a) ** 2 * (1 + 2 * np.log(a / b))) + C_3 = 0.25 * (b / a) * (((b / a) ** 2 + 1) * np.log(a / b) + (b / a) ** 2 - 1) C_5 = 0.5 * (1 - (b / a) ** 2) - C_6 = 0.25 * (b / a) * ((b / a) ** 2 - 1 + 2 * np.log10(a / b)) + C_6 = 0.25 * (b / a) * ((b / a) ** 2 - 1 + 2 * np.log(a / b)) C_8 = 0.5 * (1 + v + (1 - v) * (b / a) ** 2) - C_9 = (b / a) * (0.5 * (1 + v) * np.log10(a / b) + 0.25 * (1 - v) * (1 - (b / a) ** 2)) - L_11 = (1 / 64) * ( - 1 + 4 * (r_o / a) ** 2 - 5 * (r_o / a) ** 4 - 4 * (r_o / a) ** 2 * (2 + (r_o / a) ** 2) * np.log10(a / r_o) - ) - L_17 = 0.25 * (1 - 0.25 * (1 - v) * ((1 - (r_o / a) ** 4) - (r_o / a) ** 2 * (1 + (1 + v) * np.log10(a / r_o)))) + C_9 = (b / a) * (0.5 * (1 + v) * np.log(a / b) + 0.25 * (1 - v) * (1 - (b / a) ** 2)) + L_11 = (1 / 64) * (1 + 4 * (b / a) ** 2 - 5 * (b / a) ** 4 - 4 * (b / a) ** 2 * (2 + (b / a) ** 2) * np.log(a / b)) + L_17 = 0.25 * (1 - 0.25 * (1 - v) * ((1 - (r_o / a) ** 4) - (r_o / a) ** 2 * (1 + (1 + v) * np.log(a / r_o)))) - return D, C_2, C_3, C_5, C_6, C_8, C_9, L_11, L_17 + L_14 = 1 / 16 * (1 - (b / a) ** 2 - 4 * (b / a) ** 2 * np.log(a / b)) + + return D, C_2, C_3, C_5, C_6, C_8, C_9, L_11, L_17, L_14 class LTS_inactive_rotor(om.ExplicitComponent): @@ -164,7 +164,7 @@ def compute(self, inputs, outputs): f = Sigma_normal * r_i ** 2 * t_rdisc / (E * (h_total) * (f_d_denom1 + f_d_denom2)) u_d = ( - f / (D_L * (Lambda_L) ** 3) * ((C_14_L / (2 * C_11_L) * F_2_L - C_13_L / C_11_L * F_1_L - 0.5 * F_4_L)) + f / (2 * D_L * (Lambda_L) ** 3) * ((C_14_L / (2 * C_11_L) * F_2_L - C_13_L / C_11_L * F_1_L - 0.5 * F_4_L)) + y_sh ) @@ -176,37 +176,33 @@ def compute(self, inputs, outputs): * (((r_i) ** 2 - (R_shaft_outer) ** 2) * t_rdisc + ((r_i + h_yr_s) ** 2 - (r_i ** 2)) * l_eff_rotor) ) - # active_mass = rho_Fe * np.pi * ((r_o) ** 2 - (r_i + h_yr_s) ** 2) * L_r - # Total_mass_rotor = active_mass + outputs["Structural_mass_rotor"] - # M_r = Total_mass_rotor * g * l_eff_rotor / 2 - # delta_u = 6 * M_r * (1 - v ** 2) * (r_o) / (E * (1 + v) * t_rdisc ** 3) - - outputs["u_ar"] = abs(outputs["u_ar"] + y_sh) + outputs["u_ar"] = np.abs(outputs["u_ar"]) + y_sh - outputs["u_allowable_r"] = delta_em * u_allow_pcent / 100 + outputs["u_allowable_r"] = 0.2 * delta_em * 0.01 * u_allow_pcent - outputs["U_rotor_radial_constraint"] = outputs["u_allowable_r"] - outputs["u_ar"] + outputs["U_rotor_radial_constraint"] = np.abs(outputs["u_ar"]) / outputs["u_allowable_r"] ###################################################### Electromagnetic design############################################# # return D,C_2,C_3,C_5,C_6,C_8,C_9,L_11,L_17 # axial deformation of rotor W_back_iron = plate_constant( - r_i + h_yr_s + h_yr * 0.5, - R_shaft_outer, - v, - +h_yr_s + h_yr * 0.5, - t_rdisc, - E, - ) - W_ssteel = plate_constant( - r_i + h_yr_s * 0.5, + r_i + h_total * 0.5, R_shaft_outer, v, - h_yr + r_i + h_yr_s * 0.5, + r_i + h_total * 0.5, t_rdisc, E, ) - # W_cu = plate_constant( + + # W_ssteel = plate_constant( + # r_i + h_yr_s * 0.5, + # R_shaft_outer, + # v, + # h_yr + r_i + h_yr_s * 0.5, + # t_rdisc, + # E, + # ) + # # W_cu = plate_constant( # D_a * 0.5 - h_s * 0.5, # R_shaft_outer, # v, @@ -217,52 +213,51 @@ def compute(self, inputs, outputs): outputs["W_ry"] = rho_Fe * g * np.sin(np.deg2rad(Tilt_angle)) * (L_r - t_rdisc) * h_total - wr_disc = rho_steel * g * np.sin(np.deg2rad(Tilt_angle)) * t_rdisc - y_ai1r = ( - -0.5 - * outputs["W_ry"] - * (r_i + h_yr_s + h_yr * 0.5) ** 4 + -outputs["W_ry"] + * (r_i + h_total * 0.5) ** 4 / (R_shaft_outer * W_back_iron[0]) * (W_back_iron[1] * W_back_iron[4] / W_back_iron[3] - W_back_iron[2]) ) - W_sr = rho_steel * g * np.sin(np.deg2rad(Tilt_angle)) * (L_r - t_rdisc) * h_yr_s - y_ai2r = ( - -0.5 - * W_sr - * (r_i + h_yr_s * 0.5) ** 4 - / (R_shaft_outer * W_ssteel[0]) - * (W_ssteel[1] * W_ssteel[4] / W_ssteel[3] - W_ssteel[2]) - ) + + # W_sr = rho_steel * g * np.sin(np.deg2rad(Tilt_angle)) * (L_r - t_rdisc) * h_yr_s + # y_ai2r = ( -W_sr + # * (r_i + h_yr_s * 0.5) ** 4 + # / (R_shaft_outer * W_ssteel[0]) + # * (W_ssteel[1] * W_ssteel[4] / W_ssteel[3] - W_ssteel[2]) + # ) # W_Cu = np.sin(np.deg2rad(Tilt_angle)) * Copper / (2 * np.pi * (D_a * 0.5 - h_s * 0.5)) # y_ai3r = ( # -W_Cu * (D_a * 0.5 - h_s * 0.5) ** 4 / (R_shaft_outer * W_cu[0]) * (W_cu[1] * W_cu[4] / W_cu[3] - W_cu[2]) # ) + wr_disc = rho_steel * g * np.sin(np.deg2rad(Tilt_angle)) * t_rdisc + a_ii = r_o - h_total M_rb = ( -wr_disc * a_ii ** 2 - / W_ssteel[5] - * (W_ssteel[6] * 0.5 / (a_ii * R_shaft_outer) * (a_ii ** 2 - R_shaft_outer ** 2) - W_ssteel[8]) + / W_back_iron[3] + * (W_back_iron[4] * 0.5 / (a_ii * R_shaft_outer) * (a_ii ** 2 - R_shaft_outer ** 2) - W_back_iron[9]) ) + Q_b = wr_disc * 0.5 / R_shaft_outer * (a_ii ** 2 - R_shaft_outer ** 2) y_aiir = ( - M_rb * a_ii ** 2 / W_ssteel[0] * W_ssteel[1] - + Q_b * a_ii ** 3 / W_ssteel[0] * W_ssteel[2] - - wr_disc * a_ii ** 4 / W_ssteel[0] * W_ssteel[7] + M_rb * a_ii ** 2 / W_back_iron[0] * W_back_iron[1] + + Q_b * a_ii ** 3 / W_back_iron[0] * W_back_iron[2] + - wr_disc * a_ii ** 4 / W_back_iron[0] * W_back_iron[7] ) # I = np.pi * 0.25 * (r_i ** 4 - R_shaft_outer ** 4) # F_ecc = Sigma_normal*2*pi*K_rad*r_g**3 # M_ar = F_ecc*L_r*0.5 - outputs["y_ar"] = y_ai1r + y_ai2r + y_aiir + (r_i + h_yr + h_yr_s) * theta_sh # +M_ar*L_r**2*0/(2*E*I) + outputs["y_ar"] = y_ai1r + y_aiir + (r_i + h_yr + h_yr_s) * theta_sh # +M_ar*L_r**2*0/(2*E*I) - outputs["y_allowable_r"] = l_eff_rotor * y_allow_pcent / 100 + outputs["y_allowable_r"] = l_eff_rotor * 0.01 * y_allow_pcent # Torsional deformation of rotor - J_dr = (1 / 32) * np.pi * (r_o ** 4 - R_shaft_outer ** 4) + J_dr = (1 / 32) * np.pi * (r_i ** 4 - R_shaft_outer ** 4) J_cylr = (1 / 32) * np.pi * (r_o ** 4 - r_i ** 4) @@ -276,7 +271,7 @@ def compute(self, inputs, outputs): * (((r_i) ** 2 - (R_shaft_outer) ** 2) * t_rdisc + ((r_i + h_yr_s) ** 2 - (r_i ** 2)) * l_eff_rotor) ) - outputs["U_rotor_axial_constraint"] = outputs["y_allowable_r"] - outputs["y_ar"] + outputs["U_rotor_axial_constraint"] = np.abs(outputs["y_ar"]) / outputs["y_allowable_r"] class LTS_inactive_stator(om.ExplicitComponent): @@ -299,7 +294,7 @@ def setup(self): self.add_output("r_is", 0.0, units="m", desc="inner radius of stator disc") self.add_output("r_os", 0.0, units="m", desc="outer radius of stator disc") - self.add_output("R_sy", 0.0, units="m", desc="Stator yoke height ") + self.add_output("R_sy", 0.0, units="m", desc="Stator yoke radius ") # structural design variables self.add_input("t_sdisc", 0.0, units="m", desc="stator disc thickness") @@ -309,7 +304,7 @@ def setup(self): self.add_input("g", 9.8106, units="m/s/s", desc="Acceleration due to gravity") self.add_input("rho_steel", 0.0, units="kg/m**3", desc="Structural steel mass density") - self.add_input("mass_SC", 0.0, units="kg", desc=" SC mass") + self.add_input("mass_SC", 0.0, units="kg", desc="SC mass") self.add_input("Tilt_angle", 0.0, units="deg", desc=" Main shaft tilt angle") self.add_output("W_sy", 0.0, desc=" line load of stator yoke thickness") @@ -319,8 +314,6 @@ def setup(self): self.add_input("Sigma_normal", 0.0, units="Pa", desc="Normal stress ") # self.add_input("Sigma_shear", 0.0, units="Pa", desc="Normal stress ") - self.add_output("U_radial_stator", 0.0, units="m", desc="stator radial deflection") - self.add_output("U_axial_stator", 0.0, units="m", desc="stator axial deflection") self.add_output("U_stator_radial_constraint", 0.0, units="m", desc="Stator radial deflection contraint") self.add_output("U_stator_axial_constraint", 0.0, units="m", desc="Stator axial deflection contraint") # self.add_input("perc_allowable_radial", 0.0, desc=" Allowable radial % deflection ") @@ -328,7 +321,7 @@ def setup(self): self.add_input("Structural_mass_rotor", 0.0, units="kg", desc="rotor disc mass") self.add_output("Structural_mass_stator", 0.0, units="kg", desc="Stator mass (kg)") - self.add_output("mass_total", 0.0, units="kg", desc="stator disc mass") + self.add_output("mass_structural", 0.0, units="kg", desc="stator disc mass") # self.add_input("K_rad", desc="Aspect ratio") @@ -363,7 +356,7 @@ def compute(self, inputs, outputs): v = inputs["v"] g = inputs["g"] rho_steel = inputs["rho_steel"] - mass_SC = inputs["mass_SC"] + Total_mass_SC = inputs["mass_SC"] Tilt_angle = inputs["Tilt_angle"] Sigma_normal = inputs["Sigma_normal"] # Sigma_shear = inputs["Sigma_shear"] @@ -379,7 +372,7 @@ def compute(self, inputs, outputs): # Radial deformation of Stator L_s = l_eff_stator + t_sdisc - outputs["r_os"] = r_os = D_sc * 0.5 + h_sc + 0.25 + h_ys + outputs["r_os"] = r_os = D_sc * 0.5 + h_sc + delta_em + h_ys outputs["r_is"] = r_is = r_os - h_ys outputs["R_sy"] = (r_os + r_is) * 0.5 R_s = r_is @@ -389,6 +382,7 @@ def compute(self, inputs, outputs): f_d_denom1 = ( R_s / (E * ((R_s) ** 2 - (R_nose_outer) ** 2)) * ((1 - v) * R_s ** 2 + (1 + v) * (R_nose_outer) ** 2) ) + f_d_denom2 = ( t_sdisc / (2 * D_0 * (Lambda_0) ** 3) @@ -397,20 +391,21 @@ def compute(self, inputs, outputs): f = Sigma_normal * (R_s) ** 2 * t_sdisc / (E * (h_ys) * (f_d_denom1 + f_d_denom2)) outputs["u_as"] = ( (Sigma_normal * (R_s) ** 2) / (E * (h_ys)) - - f / (D_L * (Lambda_L) ** 3) * ((C_14_L / (2 * C_11_L) * F_2_L - C_13_L / C_11_L * F_1_L - 0.5 * F_4_L)) + - f + / (2 * D_0 * (Lambda_0) ** 3) + * ((C_14_L / (2 * C_11_L) * F_2_L - C_13_L / C_11_L * F_1_L - 0.5 * F_4_L)) + y_bd ) outputs["u_as"] += y_bd - outputs["u_allowable_s"] = delta_em * u_allow_pcent / 100 + outputs["u_allowable_s"] = 0.2 * delta_em * 0.01 * u_allow_pcent - outputs["U_stator_radial_constraint"] = outputs["u_allowable_s"] - outputs["u_as"] + outputs["U_stator_radial_constraint"] = np.abs(outputs["u_as"]) / outputs["u_allowable_s"] ###################################################### Electromagnetic design############################################# # axial deformation of stator - W_ssteel = plate_constant( R_s + h_ys * 0.5, R_nose_outer, @@ -436,7 +431,7 @@ def compute(self, inputs, outputs): * (W_ssteel[1] * W_ssteel[4] / W_ssteel[3] - W_ssteel[2]) ) - W_field = np.sin(np.deg2rad(Tilt_angle)) * mass_SC / (2 * np.pi * (D_sc * 0.5 + h_sc * 0.5)) + W_field = np.sin(np.deg2rad(Tilt_angle)) * Total_mass_SC / (2 * np.pi * (D_sc * 0.5 + h_sc * 0.5)) y_ai2s = ( -W_field * (D_sc * 0.5 + h_sc * 0.5) ** 4 @@ -451,8 +446,8 @@ def compute(self, inputs, outputs): M_sb = ( -w_disc_s * a_ii ** 2 - / W_ssteel[5] - * (W_ssteel[6] * 0.5 / (a_ii * R_nose_outer) * (a_ii ** 2 - r_oii ** 2) - W_ssteel[8]) + / W_ssteel[3] + * (W_ssteel[4] * 0.5 / (a_ii * R_nose_outer) * (a_ii ** 2 - r_oii ** 2) - W_ssteel[9]) ) Q_sb = w_disc_s * 0.5 / R_nose_outer * (a_ii ** 2 - r_oii ** 2) @@ -465,12 +460,13 @@ def compute(self, inputs, outputs): # I = np.pi * 0.25 * (R_s ** 4 - (R_nose_outer) ** 4) # F_ecc = inputs['Sigma_normal']*2*np.pi*inputs['K_rad']*inputs['r_g']**2 # M_as = F_ecc*L_s*0.5 + outputs["y_as"] = y_ai1s + y_ai2s + y_aiis + (R_s + h_ys * 0.5) * theta_bd # M_as*L_s**2*0/(2*E*I) outputs["y_allowable_s"] = L_s * y_allow_pcent / 100 # Torsional deformation of stator - J_ds = (1 / 32) * np.pi * ((r_os) ** 4 - R_nose_outer ** 4) + J_ds = (1 / 32) * np.pi * ((r_is) ** 4 - R_nose_outer ** 4) J_cyls = (1 / 32) * np.pi * ((r_os) ** 4 - r_is ** 4) @@ -483,9 +479,9 @@ def compute(self, inputs, outputs): + np.pi * (r_os ** 2 - r_is ** 2) * l_eff_stator ) - outputs["U_stator_axial_constraint"] = outputs["y_allowable_s"] - outputs["y_as"] + outputs["U_stator_axial_constraint"] = np.abs(outputs["y_as"]) / outputs["y_allowable_s"] - outputs["mass_total"] = outputs["Structural_mass_stator"] + Structural_mass_rotor + outputs["mass_structural"] = outputs["Structural_mass_stator"] + Structural_mass_rotor class LTS_Outer_Rotor_Structural(om.Group): @@ -515,9 +511,9 @@ def setup(self): ivcs.add_output("rho_steel", 0.0, units="kg/m**3", desc="Structural steel mass density") ivcs.add_output("u_allow_pcent", 0.0, desc="Radial deflection as a percentage of air gap diameter") ivcs.add_output("y_allow_pcent", 0.0, desc="Radial deflection as a percentage of air gap diameter") - ivcs.add_output("z_allow_deg", 0.0, units="deg", desc="Allowable torsional twist") - ivcs.add_output("perc_allowable_radial", 0.0, desc=" Allowable radial % deflection ") - ivcs.add_output("perc_allowable_axial", 0.0, desc=" Allowable axial % deflection ") + # ivcs.add_output("z_allow_deg", 0.0, units="deg", desc="Allowable torsional twist") + #ivcs.add_output("perc_allowable_radial", 0.0, desc=" Allowable radial % deflection ") + #ivcs.add_output("perc_allowable_axial", 0.0, desc=" Allowable axial % deflection ") self.add_subsystem("ivcs", ivcs, promotes=["*"]) self.add_subsystem("sys1", LTS_inactive_rotor(), promotes=["*"]) @@ -526,80 +522,73 @@ def setup(self): if __name__ == "__main__": - prob = om.Problem() - prob.model = LTS_Outer_Rotor_Structural() - - prob.driver = om.ScipyOptimizeDriver() # pyOptSparseDriver() - prob.driver.options["optimizer"] = "SLSQP" #'COBYLA' - prob.driver.options["maxiter"] = 10 - # prob.driver.opt_settings['IPRINT'] = 4 - # prob.driver.opt_settings['ITRM'] = 3 - # prob.driver.opt_settings['ITMAX'] = 10 - # prob.driver.opt_settings['DELFUN'] = 1e-3 - # prob.driver.opt_settings['DABFUN'] = 1e-3 - # prob.driver.opt_settings['IFILE'] = 'CONMIN_LST.out' - # prob.root.deriv_options['type']='fd' - - recorder = om.SqliteRecorder("log.sql") - prob.driver.add_recorder(recorder) - prob.add_recorder(recorder) - prob.driver.recording_options["excludes"] = ["*_df"] - prob.driver.recording_options["record_constraints"] = True - prob.driver.recording_options["record_desvars"] = True - prob.driver.recording_options["record_objectives"] = True - - prob.model.add_design_var("h_yr_s", lower=0.0250, upper=0.5, ref=0.3) - prob.model.add_design_var("h_ys", lower=0.025, upper=0.6, ref=0.35) - prob.model.add_design_var("t_rdisc", lower=0.025, upper=0.5, ref=0.3) - prob.model.add_design_var("t_sdisc", lower=0.025, upper=0.5, ref=0.3) - prob.model.add_objective("mass_total") - - prob.model.add_constraint("U_rotor_radial_constraint", lower=0.01) - prob.model.add_constraint("U_rotor_axial_constraint", lower=0.01) - prob.model.add_constraint("U_stator_radial_constraint", lower=0.01) - prob.model.add_constraint("U_stator_axial_constraint", lower=0.01) - - prob.model.approx_totals(method="fd") - - prob.setup() + prob_struct = om.Problem() + prob_struct.model = LTS_Outer_Rotor_Structural() + + prob_struct.driver = om.ScipyOptimizeDriver() + prob_struct.driver.options["optimizer"] = "SLSQP" + prob_struct.driver.options["maxiter"] = 100 + + #recorder = om.SqliteRecorder("log.sql") + #prob_struct.driver.add_recorder(recorder) + #prob_struct.add_recorder(recorder) + #prob_struct.driver.recording_options["excludes"] = ["*_df"] + #prob_struct.driver.recording_options["record_constraints"] = True + #prob_struct.driver.recording_options["record_desvars"] = True + #prob_struct.driver.recording_options["record_objectives"] = True + + prob_struct.model.add_design_var("h_yr_s", lower=0.0250, upper=0.5, ref=0.3) + prob_struct.model.add_design_var("h_ys", lower=0.025, upper=0.6, ref=0.35) + prob_struct.model.add_design_var("t_rdisc", lower=0.025, upper=0.5, ref=0.3) + prob_struct.model.add_design_var("t_sdisc", lower=0.025, upper=0.5, ref=0.3) + prob_struct.model.add_objective("mass_structural") + + prob_struct.model.add_constraint("U_rotor_radial_constraint", upper=1.0) + prob_struct.model.add_constraint("U_rotor_axial_constraint", upper=1.0) + prob_struct.model.add_constraint("U_stator_radial_constraint", upper=1.0) + prob_struct.model.add_constraint("U_stator_axial_constraint", upper=1.0) + + prob_struct.model.approx_totals(method="fd") + + prob_struct.setup() # --- Design Variables --- # Initial design variables for a PMSG designed for a 15MW turbine - prob["Sigma_shear"] = 74.99692029e3 - prob["Sigma_normal"] = 378.45123826e3 - prob["T_e"] = 9e06 - prob["l_eff_stator"] = 1.44142189 # rev 1 9.94718e6 - prob["l_eff_rotor"] = 1.2827137 - prob["D_a"] = 7.74736313 - prob["delta_em"] = 0.0199961 - prob["h_s"] = 0.1803019703 - prob["D_sc"] = 7.78735533 - prob["rho_steel"] = 7850 - prob["rho_Fe"] = 7700 - prob["Tilt_angle"] = 90.0 - prob["R_shaft_outer"] = 1.25 - prob["R_nose_outer"] = 0.95 - prob["u_allow_pcent"] = 50 - prob["y_allow_pcent"] = 20 - prob["h_yr"] = 0.1254730934 - prob["h_yr_s"] = 0.025 - prob["h_ys"] = 0.050 - prob["t_rdisc"] = 0.05 - prob["t_sdisc"] = 0.100 - prob["y_bd"] = 0.00 - prob["theta_bd"] = 0.00 - prob["y_sh"] = 0.00 - prob["theta_sh"] = 0.00 - - prob["Copper"] = 60e3 - prob["mass_SC"] = 4000 - - prob.model.approx_totals(method="fd") - - prob.run_model() - # prob.run_driver() - - # prob.model.list_outputs(values = True, hierarchical=True) + #prob_struct["Sigma_shear"] = 74.99692029e3 + prob_struct["Sigma_normal"] = 378.45123826e3 + prob_struct["T_e"] = 9e06 + prob_struct["l_eff_stator"] = 1.44142189 # rev 1 9.94718e6 + prob_struct["l_eff_rotor"] = 1.2827137 + prob_struct["D_a"] = 7.74736313 + prob_struct["delta_em"] = 0.0199961 + prob_struct["h_s"] = 0.1803019703 + prob_struct["D_sc"] = 7.78735533 + prob_struct["rho_steel"] = 7850 + prob_struct["rho_Fe"] = 7700 + prob_struct["Tilt_angle"] = 90.0 + prob_struct["R_shaft_outer"] = 1.25 + prob_struct["R_nose_outer"] = 0.95 + prob_struct["u_allow_pcent"] = 50 + prob_struct["y_allow_pcent"] = 20 + prob_struct["h_yr"] = 0.1254730934 + prob_struct["h_yr_s"] = 0.025 + prob_struct["h_ys"] = 0.050 + prob_struct["t_rdisc"] = 0.05 + prob_struct["t_sdisc"] = 0.1 + prob_struct["y_bd"] = 0.00 + prob_struct["theta_bd"] = 0.00 + prob_struct["y_sh"] = 0.00 + prob_struct["theta_sh"] = 0.00 + + #prob_struct["mass_copper"] = 60e3 + prob_struct["mass_SC"] = 4000 + + prob_struct.model.approx_totals(method="fd") + + #prob_struct.run_model() + prob_struct.run_driver() + + prob_struct.model.list_outputs(val=True, hierarchical=True) raw_data = { "Parameters": [ @@ -616,27 +605,27 @@ def setup(self): "Total structural mass", ], "Values": [ - prob.get_val("t_rdisc", units="mm"), - prob.get_val("h_yr_s", units="mm"), - prob.get_val("t_sdisc", units="mm"), - prob.get_val("h_ys", units="mm"), - prob.get_val("u_ar", units="mm"), - prob.get_val("y_ar", units="mm"), - prob.get_val("u_as", units="mm"), - prob.get_val("y_as", units="mm"), - prob.get_val("Structural_mass_rotor", units="t"), - prob.get_val("Structural_mass_stator", units="t"), - prob.get_val("mass_total", units="t"), + prob_struct.get_val("t_rdisc", units="mm"), + prob_struct.get_val("h_yr_s", units="mm"), + prob_struct.get_val("t_sdisc", units="mm"), + prob_struct.get_val("h_ys", units="mm"), + prob_struct.get_val("u_ar", units="mm"), + prob_struct.get_val("y_ar", units="mm"), + prob_struct.get_val("u_as", units="mm"), + prob_struct.get_val("y_as", units="mm"), + prob_struct.get_val("Structural_mass_rotor", units="t"), + prob_struct.get_val("Structural_mass_stator", units="t"), + prob_struct.get_val("mass_structural", units="t"), ], "Limit": [ "", "", "", "", - prob.get_val("u_allowable_r", units="mm"), - prob.get_val("y_allowable_r", units="mm"), - prob.get_val("u_allowable_s", units="mm"), - prob.get_val("y_allowable_s", units="mm"), + prob_struct.get_val("u_allowable_r", units="mm"), + prob_struct.get_val("y_allowable_r", units="mm"), + prob_struct.get_val("u_allowable_s", units="mm"), + prob_struct.get_val("y_allowable_s", units="mm"), "", "", "", @@ -655,9 +644,7 @@ def setup(self): "tons", ], } - # print(raw_data) df = pd.DataFrame(raw_data, columns=["Parameters", "Values", "Limit", "Units"]) - + #df.to_excel("Optimized_structure_LTSG_MW.xlsx") print(df) - df.to_excel("Optimized_structure_LTSG_MW.xlsx")