diff --git a/WISDEM/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml b/WISDEM/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml index 603f2eccd..507675b9f 100644 --- a/WISDEM/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml +++ b/WISDEM/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml @@ -690,6 +690,8 @@ components: joint2: col3_lower_pontoon outer_shape: *pontlow_out internal_structure: *pontlow_int + transition_piece_mass: 100.e3 + transition_piece_cost: 300.e3 mooring: nodes: - name: line1_anchor diff --git a/WISDEM/examples/16_inverse_design/analysis_options.yaml b/WISDEM/examples/16_inverse_design/analysis_options.yaml new file mode 100644 index 000000000..97c348b1a --- /dev/null +++ b/WISDEM/examples/16_inverse_design/analysis_options.yaml @@ -0,0 +1,50 @@ +general: + folder_output: outputs + fname_output: refturb_output + +design_variables: + floating: + joints: + flag: True + z_coordinate: + - names: [spar_keel] + lower_bound: -40.0 + upper_bound: -15.0 + - names: [spar_freeboard] + lower_bound: -40.0 + upper_bound: -15.0 + members: + flag: True + groups: + - names: [spar] + ballast: + lower_bound: 1.0 + upper_bound: 1e4 + + mooring: + line_length: + flag: True + lower_bound: 100.0 + upper_bound: 1000.0 + +merit_figure: inverse_design + +inverse_design: + floatingse.platform_mass: + ref_value: 6.e6 + units: kg + floatingse.mooring_mass: + ref_value: 4.5e5 + units: kg + floatingse.system_structural_center_of_mass: + ref_value: [0., -60.] + indices: [1, 2] + units: m + +driver: + optimization: + flag: True # Flag to enable optimization + tol: 1.e-4 # Optimality tolerance + solver: SLSQP # Optimization solver. Other options are 'SLSQP' - 'CONMIN' + step_size: 1.e-6 # Step size for finite differencing + form: forward # Finite differencing mode, either forward or central diff --git a/WISDEM/examples/16_inverse_design/inverse_spar_design.py b/WISDEM/examples/16_inverse_design/inverse_spar_design.py new file mode 100644 index 000000000..133b42d22 --- /dev/null +++ b/WISDEM/examples/16_inverse_design/inverse_spar_design.py @@ -0,0 +1,16 @@ +import os +import sys +import time + +from wisdem.commonse.mpi_tools import MPI +from wisdem.glue_code.runWISDEM import run_wisdem + +## File management +run_dir = os.path.dirname(os.path.realpath(__file__)) + os.sep +wisdem_examples = os.path.dirname(os.path.dirname(run_dir)) +fname_wt_input_oc3 = os.path.join(wisdem_examples, "09_floating", "nrel5mw-spar_oc3.yaml") +fname_modeling_options = os.path.join(wisdem_examples, "09_floating", "modeling_options_noRNA.yaml") +fname_analysis_options = run_dir + os.sep + "analysis_options.yaml" + + +wt_opt, modeling_options, opt_options = run_wisdem(fname_wt_input_oc3, fname_modeling_options, fname_analysis_options) diff --git a/WISDEM/wisdem/floatingse/floating_frame.py b/WISDEM/wisdem/floatingse/floating_frame.py index 10443c48c..a5c4829c5 100644 --- a/WISDEM/wisdem/floatingse/floating_frame.py +++ b/WISDEM/wisdem/floatingse/floating_frame.py @@ -12,9 +12,6 @@ RIGID = 1e30 EPS = 1e-6 -# TODO: -# - Added mass, hydro stiffness for tower sim - class PlatformFrame(om.ExplicitComponent): def initialize(self): @@ -58,6 +55,12 @@ def setup(self): self.add_input(f"member{k}:Pz", np.zeros(MEMMAX), units="N/m") self.add_input(f"member{k}:qdyn", np.zeros(MEMMAX), units="Pa") + self.add_input("transition_node", np.zeros(3), units="m") + self.add_input("transition_piece_mass", 0.0, units="kg") + self.add_input("transition_piece_cost", 0.0, units="USD") + + self.add_output("transition_piece_I", np.zeros(6), units="kg*m**2") + self.add_output("platform_nodes", NULL * np.ones((NNODES_MAX, 3)), units="m") self.add_output("platform_Fnode", NULL * np.ones((NNODES_MAX, 3)), units="N") self.add_output("platform_Rnode", NULL * np.ones(NNODES_MAX), units="m") @@ -85,12 +88,12 @@ def setup(self): self.add_discrete_output("platform_elem_memid", [-1] * NELEM_MAX) self.add_output("platform_displacement", 0.0, units="m**3") self.add_output("platform_center_of_buoyancy", np.zeros(3), units="m") - self.add_output("platform_center_of_mass", np.zeros(3), units="m") + self.add_output("platform_hull_center_of_mass", np.zeros(3), units="m") self.add_output("platform_centroid", np.zeros(3), units="m") self.add_output("platform_ballast_mass", 0.0, units="kg") self.add_output("platform_hull_mass", 0.0, units="kg") self.add_output("platform_mass", 0.0, units="kg") - self.add_output("platform_I_total", np.zeros(6), units="kg*m**2") + self.add_output("platform_I_hull", np.zeros(6), units="kg*m**2") self.add_output("platform_cost", 0.0, units="USD") self.add_output("platform_Awater", 0.0, units="m**2") self.add_output("platform_Iwater", 0.0, units="m**4") @@ -98,12 +101,11 @@ def setup(self): self.add_output("platform_variable_capacity", np.zeros(n_member), units="m**3") self.node_mem2glob = {} - # self.node_glob2mem = {} def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): - # This shouldn't change during an optimization, so save some time? - if len(self.node_mem2glob) == 0: - self.set_connectivity(inputs, outputs) + # Seems like we have to run this each time as numbering can change during optimization + self.node_mem2glob = {} + self.set_connectivity(inputs, outputs) self.set_node_props(inputs, outputs) self.set_element_props(inputs, outputs, discrete_inputs, discrete_outputs) @@ -157,7 +159,9 @@ def set_node_props(self, inputs, outputs): n_member = opt["floating"]["members"]["n_members"] # Number of valid nodes - nnode = np.where(outputs["platform_nodes"][:, 0] == NULL)[0][0] + node_platform = outputs["platform_nodes"] + nnode = np.where(node_platform[:, 0] == NULL)[0][0] + node_platform = node_platform[:nnode, :] # Find greatest radius of all members at node intersections Rnode = np.zeros(nnode) @@ -175,6 +179,13 @@ def set_node_props(self, inputs, outputs): iglob = self.node_mem2glob[(k, icb)] Fnode[iglob, 2] += inputs[f"member{k}:buoyancy_force"] + # Get transition piece inertial properties + itrans_platform = util.closest_node(node_platform, inputs["transition_node"]) + m_trans = float(inputs["transition_piece_mass"]) + r_trans = Rnode[itrans_platform] + I_trans = m_trans * r_trans ** 2.0 * np.r_[0.5, 0.5, 1.0, np.zeros(3)] + outputs["transition_piece_I"] = I_trans + # Store outputs outputs["platform_Rnode"] = NULL * np.ones(NNODES_MAX) outputs["platform_Rnode"][:nnode] = Rnode @@ -265,13 +276,21 @@ def set_element_props(self, inputs, outputs, discrete_inputs, discrete_outputs): cg_plat += imass * inputs[f"member{k}:center_of_mass"] cb_plat += ivol * inputs[f"member{k}:center_of_buoyancy"] + # Add transition piece + m_trans = inputs["transition_piece_mass"] + cg_trans = inputs["transition_node"] + I_trans = util.assembleI(outputs["transition_piece_I"]) + mass += m_trans + cost += inputs["transition_piece_cost"] + cg_plat += m_trans * cg_trans + # Finalize outputs cg_plat /= mass cb_plat /= volume # With CG known, loop back through to compute platform I unit_z = np.array([0.0, 0.0, 1.0]) - I_total = np.zeros((3, 3)) + I_hull = np.zeros((3, 3)) for k in range(n_member): xyz_k = inputs[f"member{k}:nodes_xyz"] inodes = np.where(xyz_k[:, 0] == NULL)[0][0] @@ -287,10 +306,14 @@ def set_element_props(self, inputs, outputs, discrete_inputs, discrete_outputs): # Rotate member inertia tensor I_k = util.assembleI(inputs[f"member{k}:I_total"]) - I_k2 = T @ I_k @ T.T + I_k_rot = T @ I_k @ T.T # Now do parallel axis theorem - I_total += np.array(I_k2) + imass * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) + I_hull += np.array(I_k_rot) + imass * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) + + # Add in transition piece + R = cg_plat - cg_trans + I_hull += I_trans + m_trans * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) # Store outputs nelem = elem_A.size @@ -340,9 +363,9 @@ def set_element_props(self, inputs, outputs, discrete_inputs, discrete_outputs): outputs["platform_hull_mass"] = mass - m_ball outputs["platform_cost"] = cost outputs["platform_displacement"] = volume - outputs["platform_center_of_mass"] = cg_plat + outputs["platform_hull_center_of_mass"] = cg_plat outputs["platform_center_of_buoyancy"] = cb_plat - outputs["platform_I_total"] = util.unassembleI(I_total) + outputs["platform_I_hull"] = util.unassembleI(I_hull) outputs["platform_Awater"] = Awater outputs["platform_Iwater"] = Iwater outputs["platform_added_mass"] = m_added @@ -357,7 +380,8 @@ def setup(self): def compute(self, inputs, outputs): transition_node = inputs["transition_node"] - tower_top_node = transition_node + tower_top_node = 0 #previous code altered the original definition of transition_node + tower_top_node += transition_node tower_top_node[2] += float(inputs["tower_height"]) outputs["tower_top_node"] = tower_top_node @@ -395,8 +419,9 @@ def setup(self): self.add_input("platform_elem_Pz1", NULL * np.ones(NELEM_MAX), units="N/m") self.add_input("platform_elem_Pz2", NULL * np.ones(NELEM_MAX), units="N/m") self.add_input("platform_elem_qdyn", NULL * np.ones(NELEM_MAX), units="Pa") - self.add_input("platform_center_of_mass", np.zeros(3), units="m") + self.add_input("platform_hull_center_of_mass", np.zeros(3), units="m") self.add_input("platform_mass", 0.0, units="kg") + self.add_input("platform_I_hull", np.zeros(6), units="kg*m**2") self.add_input("platform_displacement", 0.0, units="m**3") self.add_input("tower_nodes", NULL * np.ones((MEMMAX, 3)), units="m") @@ -433,7 +458,6 @@ def setup(self): self.add_input("rho_water", 0.0, units="kg/m**3") self.add_input("tower_top_node", np.zeros(3), units="m") self.add_input("transition_node", np.zeros(3), units="m") - self.add_input("transition_piece_mass", 0.0, units="kg") self.add_input("rna_mass", 0.0, units="kg") self.add_input("rna_cg", np.zeros(3), units="m") self.add_input("mooring_neutral_load", np.zeros((n_attach, 3)), units="N") @@ -478,7 +502,8 @@ def setup(self): self.add_output("constr_variable_margin", val=0.0) self.add_output("member_variable_volume", val=np.zeros(n_member), units="m**3") self.add_output("member_variable_height", val=np.zeros(n_member)) - self.add_output("transition_piece_I", np.zeros(6), units="kg*m**2") + self.add_output("platform_total_center_of_mass", np.zeros(3), units="m") + self.add_output("platform_I_total", np.zeros(6), units="kg*m**2") def compute(self, inputs, outputs): # Combine nodes @@ -576,17 +601,17 @@ def compute(self, inputs, outputs): # Mass summaries m_platform = inputs["platform_mass"] + cg_platform = inputs["platform_hull_center_of_mass"] + I_platform = util.assembleI(inputs["platform_I_hull"]) m_tower = inputs["tower_mass"] m_rna = inputs["rna_mass"] - m_trans = inputs["transition_piece_mass"] - m_sys = m_platform + m_tower + m_rna + m_trans + m_sys = m_platform + m_tower + m_rna outputs["system_structural_mass"] = m_sys outputs["system_structural_center_of_mass"] = ( - m_platform * inputs["platform_center_of_mass"] + m_platform * cg_platform + m_tower * inputs["tower_center_of_mass"] + m_rna * (inputs["rna_cg"] + inputs["tower_top_node"]) - + m_trans * inputs["transition_node"] ) / m_sys # Balance out variable ballast @@ -600,6 +625,7 @@ def compute(self, inputs, outputs): outputs["constr_variable_margin"] = V_variable / capacity_sum V_variable_member = V_variable * capacity / capacity_sum outputs["member_variable_volume"] = V_variable_member + m_variable_member = V_variable_member * rho_water # Now find the CG of the variable mass assigned to each member n_member = capacity.size @@ -632,11 +658,51 @@ def compute(self, inputs, outputs): m_sys * outputs["system_structural_center_of_mass"] + m_variable * cg_variable ) / (m_sys + m_variable) - # Transition piece properties - m_trans = float(inputs["transition_piece_mass"]) - r_trans = inputs["platform_Rnode"][itrans_platform] - I_trans = m_trans * r_trans ** 2.0 * np.r_[0.5, 0.5, 1.0, np.zeros(3)] - outputs["transition_piece_I"] = I_trans + # Compute the total cg for the platform and the variable ballast together using a weighted sum approach + cg_plat_total = (m_variable * cg_variable + m_platform * cg_platform) / (m_variable + m_platform) + outputs["platform_total_center_of_mass"] = cg_plat_total + + # Now loop again to compute variable I + unit_z = np.array([0.0, 0.0, 1.0]) + I_variable = np.zeros((3, 3)) + for k in range(n_member): + if V_variable_member[k] == 0.0: + continue + + xyz = inputs[f"member{k}:nodes_xyz"] + inodes = np.where(xyz[:, 0] == NULL)[0][0] + xyz = xyz[:inodes, :] + vec_k = xyz[-1, :] - xyz[0, :] + + ds = outputs["member_variable_height"][k] + + # Compute I aligned with member + h_k = ds * np.sqrt(np.sum(vec_k ** 2)) + if h_k == 0.0: + continue + r_k = np.sqrt(V_variable_member[k] / h_k / np.pi) + I_k = ( + m_variable_member[k] * np.r_[(3 * r_k ** 2 + h_k ** 2) / 12.0 * np.ones(2), 0.5 * r_k ** 2, np.ones(3)] + ) + + # Rotate I to global c.s. + T = util.rotate_align_vectors(vec_k, unit_z) + I_k_rot = T @ util.assembleI(I_k) @ T.T + + # Now do parallel axis theorem + R = cg_variable - cg_variable_member[k, :] + I_variable += np.array(I_k_rot) + m_variable_member[k] * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) + + # Find platform I with variable contribution + I_total = np.zeros((3, 3)) + + # Compute the full moment of inertia for the platform and variable ballast + R = cg_plat_total - cg_platform + I_total += I_platform + m_platform * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) + + R = cg_plat_total - cg_variable + I_total += I_variable + m_variable * (np.dot(R, R) * np.eye(3) - np.outer(R, R)) + outputs["platform_I_total"] = util.unassembleI(I_total) class FrameAnalysis(om.ExplicitComponent): @@ -648,7 +714,7 @@ def setup(self): n_attach = opt["mooring"]["n_attach"] self.add_input("platform_mass", 0.0, units="kg") - self.add_input("platform_center_of_mass", np.zeros(3), units="m") + self.add_input("platform_hull_center_of_mass", np.zeros(3), units="m") self.add_input("platform_added_mass", np.zeros(6), units="kg") self.add_input("platform_center_of_buoyancy", np.zeros(3), units="m") self.add_input("platform_displacement", 0.0, units="m**3") @@ -818,8 +884,9 @@ def compute(self, inputs, outputs): m_trans = float(inputs["transition_piece_mass"]) if frame == "tower": m_trans += float(inputs["platform_mass"]) + inputs["platform_added_mass"][0] + m_variable - cg_trans = inputs["transition_node"] - inputs["platform_center_of_mass"] + cg_trans = inputs["transition_node"] - inputs["platform_hull_center_of_mass"] I_trans[:3] += inputs["platform_added_mass"][3:] + else: m_trans += m_variable cg_trans = np.zeros(3) @@ -885,6 +952,7 @@ def compute(self, inputs, outputs): # Add the load case and run myframe.addLoadCase(load_obj) # myframe.write(f"{frame}.3dd") + # myframe.draw() displacements, forces, reactions, internalForces, mass, modal = myframe.run() # natural frequncies diff --git a/WISDEM/wisdem/glue_code/gc_PoseOptimization.py b/WISDEM/wisdem/glue_code/gc_PoseOptimization.py index 14b786efc..ef43544b8 100644 --- a/WISDEM/wisdem/glue_code/gc_PoseOptimization.py +++ b/WISDEM/wisdem/glue_code/gc_PoseOptimization.py @@ -432,8 +432,8 @@ def set_objective(self, wt_opt): else: wt_opt.model.add_objective("floatingse.tower_mass", ref=1e6) - elif self.opt["merit_figure"] == "mononpile_mass": - wt_opt.model.add_objective("towerse.mononpile_mass", ref=1e6) + elif self.opt["merit_figure"] == "monopile_mass": + wt_opt.model.add_objective("towerse.monopile_mass", ref=1e6) elif self.opt["merit_figure"] == "structural_mass": if not self.modeling["flags"]["floating"]: @@ -470,6 +470,10 @@ def set_objective(self, wt_opt): wt_opt.model.add_objective("rotorse.rp.powercurve.Cp_regII", ref=-1.0) else: wt_opt.model.add_objective("rotorse.ccblade.CP", ref=-1.0) + + elif self.opt["merit_figure"] == "inverse_design": + wt_opt.model.add_objective("inverse_design.objective") + else: raise ValueError("The merit figure " + self.opt["merit_figure"] + " is unknown or not supported.") diff --git a/WISDEM/wisdem/glue_code/gc_WT_DataStruc.py b/WISDEM/wisdem/glue_code/gc_WT_DataStruc.py index e2e1e0983..1ed5bb7fd 100644 --- a/WISDEM/wisdem/glue_code/gc_WT_DataStruc.py +++ b/WISDEM/wisdem/glue_code/gc_WT_DataStruc.py @@ -1,9 +1,10 @@ import copy import numpy as np +from scipy.interpolate import PchipInterpolator, interp1d + import openmdao.api as om import wisdem.moorpy.MoorProps as mp -from scipy.interpolate import PchipInterpolator, interp1d from wisdem.commonse.utilities import arc_length, arc_length_deriv from wisdem.rotorse.parametrize_rotor import ParametrizeBladeAero, ParametrizeBladeStruct from wisdem.rotorse.geometry_tools.geometry import remap2grid, trailing_edge_smoothing @@ -2128,6 +2129,8 @@ def setup(self): jivc = self.add_subsystem("joints", om.IndepVarComp(), promotes=["*"]) jivc.add_output("location_in", val=np.zeros((n_joints, 3)), units="m") jivc.add_output("transition_node", val=np.zeros(3), units="m") + jivc.add_output("transition_piece_mass", val=0.0, units="kg", desc="point mass of transition piece") + jivc.add_output("transition_piece_cost", val=0.0, units="USD", desc="cost of transition piece") # Additions for optimizing individual nodes or multiple nodes concurrently self.add_subsystem("nodedv", NodeDVs(options=floating_init_options["joints"]), promotes=["*"]) diff --git a/WISDEM/wisdem/glue_code/gc_WT_InitModel.py b/WISDEM/wisdem/glue_code/gc_WT_InitModel.py index 4cc03975a..dd2820e0b 100644 --- a/WISDEM/wisdem/glue_code/gc_WT_InitModel.py +++ b/WISDEM/wisdem/glue_code/gc_WT_InitModel.py @@ -922,7 +922,8 @@ def assign_floating_values(wt_opt, modeling_options, floating): else: itrans = modeling_options["floating"]["transition_joint"] wt_opt["floating.transition_node"] = wt_opt["floating.location_in"][itrans, :] - + wt_opt["floating.transition_piece_mass"] = floating["transition_piece_mass"] + wt_opt["floating.transition_piece_cost"] = floating["transition_piece_cost"] # Make sure IVCs are initialized too for k, linked_node_dict in enumerate(modeling_options["floating"]["joints"]["design_variable_data"]): idx = linked_node_dict["indices"] @@ -1139,7 +1140,7 @@ def assign_environment_values(wt_opt, environment, offshore, blade_flag): wt_opt["env.G_soil"] = environment["soil_shear_modulus"] wt_opt["env.nu_soil"] = environment["soil_poisson"] if blade_flag: - wt_opt['rotorse.wt_class.V_mean_overwrite'] = environment['V_mean'] + wt_opt["rotorse.wt_class.V_mean_overwrite"] = environment["V_mean"] return wt_opt diff --git a/WISDEM/wisdem/glue_code/glue_code.py b/WISDEM/wisdem/glue_code/glue_code.py index aac03bc68..ce2200f43 100644 --- a/WISDEM/wisdem/glue_code/glue_code.py +++ b/WISDEM/wisdem/glue_code/glue_code.py @@ -1,4 +1,5 @@ import numpy as np + import openmdao.api as om from wisdem.rotorse.rotor import RotorSE from wisdem.towerse.tower import TowerSE @@ -423,6 +424,8 @@ def setup(self): self.connect("tower.outfitting_factor", "floatingse.tower.outfitting_factor_in") self.connect("tower.layer_mat", "floatingse.tower.layer_materials") self.connect("floating.transition_node", "floatingse.transition_node") + self.connect("floating.transition_piece_mass", "floatingse.transition_piece_mass") + self.connect("floating.transition_piece_cost", "floatingse.transition_piece_cost") if modeling_options["flags"]["tower"]: self.connect("tower_grid.height", "floatingse.tower_height") if modeling_options["flags"]["nacelle"]: @@ -556,6 +559,87 @@ def setup(self): self.connect("costs.controls_machine_rating_cost_coeff", "tcc.controls_machine_rating_cost_coeff") self.connect("costs.crane_cost", "tcc.crane_cost") + # Final component for inverse design objective + if opt_options["inverse_design"]: + self.add_subsystem("inverse_design", InverseDesign(opt_options=opt_options)) + + for name in opt_options["inverse_design"]: + indices = opt_options["inverse_design"][name]["indices"] + short_name = name.replace(".", "_") + self.connect(name, f"inverse_design.{short_name}", src_indices=indices) + + +class InverseDesign(om.ExplicitComponent): + """ + Component that takes in an arbitrary set of user-defined inputs and computes + the root-mean-square (RMS) difference between the values in the model and + a set of reference values. + + This is useful for inverse design problems where we are trying to design a + wind turbine system that has a certain set of properties. Specifically, we + might be trying to match performance values from a report by allowing the + optimizer to select the design variable values that most closely produce a + system that has those properties. + + """ + + def initialize(self): + self.options.declare("opt_options") + + def setup(self): + opt_options = self.options["opt_options"] + + # Loop through all of the keys in the inverse_design definition + for name in opt_options["inverse_design"]: + item = opt_options["inverse_design"][name] + + indices = item["indices"] + + # Grab the short name for each parameter to match + short_name = name.replace(".", "_") + + # Only apply units if they're provided by the user + if "units" in item: + units = item["units"] + else: + units = None + + self.add_input( + short_name, + val=np.zeros(len(indices)), + units=units, + ) + + # Create a singular output called objective + self.add_output( + "objective", + val=0.0, + ) + + def compute(self, inputs, outputs): + opt_options = self.options["opt_options"] + + total = 0.0 + # Loop through all of the keys in the inverse_design definition + for name in opt_options["inverse_design"]: + + item = opt_options["inverse_design"][name] + + # Grab the short name for each parameter to match + short_name = name.replace(".", "_") + + # Grab the reference value provided by the user + ref_value = item["ref_value"] + + # Compute the mean square difference between the parameter + # value outputted from the model and the reference value. Sum this + # to `total` to get the total across all parameters + total += np.sum(((inputs[short_name] - ref_value) / (np.abs(ref_value) + 1.0)) ** 2) + + # Take the square root of the total + rms_total = np.sqrt(total) + outputs["objective"] = rms_total + class WindPark(om.Group): # Openmdao group to run the cost analysis of a wind park @@ -610,6 +694,9 @@ def setup(self): self.connect("mooring.line_diameter", "orbit.mooring_line_diameter", src_indices=[0]) self.connect("mooring.unstretched_length", "orbit.mooring_line_length", src_indices=[0]) self.connect("mooring.anchor_mass", "orbit.anchor_mass", src_indices=[0]) + self.connect("floating.transition_piece_mass", "orbit.transition_piece_mass") + self.connect("floating.transition_piece_cost", "orbit.transition_piece_cost") + self.connect("floatingse.platform_cost", "orbit.floating_substructure_cost") self.connect("rotorse.re.precomp.blade_mass", "orbit.blade_mass") self.connect("tcc.turbine_cost_kW", "orbit.turbine_capex") if modeling_options["flags"]["nacelle"]: diff --git a/WISDEM/wisdem/inputs/analysis_schema.yaml b/WISDEM/wisdem/inputs/analysis_schema.yaml index 948c80097..98249ecf7 100644 --- a/WISDEM/wisdem/inputs/analysis_schema.yaml +++ b/WISDEM/wisdem/inputs/analysis_schema.yaml @@ -731,20 +731,6 @@ properties: flag: *flag lower_bound: *lboundm upper_bound: *uboundm - line_mass_density_coeff: - type: object - default: {} - properties: - flag: *flag - lower_bound: *lboundm - upper_bound: *uboundm - line_stiffness_coeff: - type: object - default: {} - properties: - flag: *flag - lower_bound: *lboundm - upper_bound: *uboundm constraints: type: object @@ -1181,9 +1167,29 @@ properties: merit_figure: type: string - description: Objective function / merit figure for optimization (not checked via schema, but in python). Choices are LCOE- levelized cost of energy, AEP- turbine annual energy production, Cp- rotor power coefficient, blade_mass, tower_mass, tower_cost, monopile_mass, monopile_cost, structural_mass- tower+monopile mass, structural_cost- tower+monopile cost, blade_tip_deflection- blade tip deflection distance towards tower, My_std- blade flap moment standard deviation, flp1_std- trailing flap standard deviation + description: Objective function / merit figure for optimization (not checked via schema, but in python). Choices are LCOE- levelized cost of energy, AEP- turbine annual energy production, Cp- rotor power coefficient, blade_mass, tower_mass, tower_cost, monopile_mass, monopile_cost, structural_mass- tower+monopile mass, structural_cost- tower+monopile cost, blade_tip_deflection- blade tip deflection distance towards tower, My_std- blade flap moment standard deviation, flp1_std- trailing flap standard deviation, platform_mass- floating platform mass without variable ballast, inverse_design- custom inverse design used to match reference values of any arbitrary output default: LCOE + inverse_design: + type: object + description: For use with the inverse_design merit_figure. Specifies the reference output variable's 'prom_name' name and the desired value, accepts multiple variables. A normalized difference between the actual value and reference value is calculated for each variable. A Root Mean Square (RMS) is calculated with all variables and the optimizer minimizes the RMS. If the refernce output variable is an array, specify the element index number via "idx". + default: {} + additionalProperties: + type: object + required: + - ref_value + optional: + - indices + - units + properties: + ref_value: + type: [number, array] + indices: + type: array + default: [0] + units: + type: string + driver: type: object default: {} diff --git a/WISDEM/wisdem/inputs/geometry_schema.yaml b/WISDEM/wisdem/inputs/geometry_schema.yaml index c793fc885..a1f272df5 100644 --- a/WISDEM/wisdem/inputs/geometry_schema.yaml +++ b/WISDEM/wisdem/inputs/geometry_schema.yaml @@ -1239,6 +1239,8 @@ properties: - members optional: - rigid_bodies + - transition_piece_mass + - transition_piece_cost properties: joints: type: array @@ -1613,6 +1615,18 @@ properties: Ca: *ca Cp: *cp Cd: *cd + transition_piece_mass: + type: number + description: Total mass of transition piece + unit: kg + minimum: 0.0 + default: 0.0 + transition_piece_cost: + type: number + description: Total cost of transition piece + unit: USD + minimum: 0.0 + default: 0.0 mooring: # Will insert mooring schema that we agreed to here or will keep as separate file description: Ontology definition for mooring systems suitable for use with the WEIS co-design analysis tool diff --git a/WISDEM/wisdem/pyframe3dd/pyframe3dd.py b/WISDEM/wisdem/pyframe3dd/pyframe3dd.py index bacd8f6e3..594fd728c 100644 --- a/WISDEM/wisdem/pyframe3dd/pyframe3dd.py +++ b/WISDEM/wisdem/pyframe3dd/pyframe3dd.py @@ -1194,6 +1194,31 @@ def write(self, fname): f.write("# End of input data file\n") f.close() + def draw(self): + # Visualization for debugging + import matplotlib.pyplot as plt + + nnode = len(self.nx) + node_array = np.zeros((nnode, 3)) + mynodes = {} + for k in range(nnode): + temp = np.r_[self.nx[k], self.ny[k], self.nz[k]] + mynodes[self.nnode[k]] = temp + node_array[k, :] = temp + myelem = [] + for k in range(len(self.eN1)): + myelem.append((self.eN1[k], self.eN2[k])) + + fig = plt.figure() + ax = fig.add_subplot(111, projection="3d") + for e in myelem: + xs = np.array([mynodes[e[0]][0], mynodes[e[1]][0]]) + ys = np.array([mynodes[e[0]][1], mynodes[e[1]][1]]) + zs = np.array([mynodes[e[0]][2], mynodes[e[1]][2]]) + ax.plot(xs, ys, zs, "b-") + ax.plot(node_array[:, 0], node_array[:, 1], node_array[:, 2], ".k", markersize=10) + plt.show() + class StaticLoadCase(object): """docstring""" diff --git a/WISDEM/wisdem/test/test_examples/test_examples.py b/WISDEM/wisdem/test/test_examples/test_examples.py index 7a8eddc93..f8fc40952 100644 --- a/WISDEM/wisdem/test/test_examples/test_examples.py +++ b/WISDEM/wisdem/test/test_examples/test_examples.py @@ -49,6 +49,7 @@ "13_design_of_experiments/doe_driver", "14_overridden_values/driver", "15_step_size_study/driver", + "16_inverse_design/inverse_spar_design", ] diff --git a/WISDEM/wisdem/test/test_floatingse/test_frame.py b/WISDEM/wisdem/test/test_floatingse/test_frame.py index e2d97c025..36d8be970 100644 --- a/WISDEM/wisdem/test/test_floatingse/test_frame.py +++ b/WISDEM/wisdem/test/test_floatingse/test_frame.py @@ -146,6 +146,7 @@ def setUp(self): self.inputs["rna_F"] = np.array([1e2, 1e1, 0.0]) self.inputs["rna_M"] = np.array([2e1, 2e2, 0.0]) self.inputs["transition_piece_mass"] = 1e3 + self.inputs["transition_piece_cost"] = 3e3 self.inputs["rho_water"] = 1e3 def testTetrahedron(self): @@ -222,16 +223,18 @@ def testTetrahedron(self): npt.assert_equal(self.outputs["platform_center_of_buoyancy"], centroid) npt.assert_equal(self.outputs["platform_centroid"], centroid) - npt.assert_equal(self.outputs["platform_center_of_mass"], centroid) - self.assertEqual(self.outputs["platform_mass"], 6e3) + cg = (6e3 * centroid + 1e3 * np.array([0.0, 0.0, 1.0])) / 7e3 + npt.assert_equal(self.outputs["platform_hull_center_of_mass"], cg) + self.assertEqual(self.outputs["platform_mass"], 6e3 + 1e3) self.assertEqual(self.outputs["platform_ballast_mass"], 6e2) - self.assertEqual(self.outputs["platform_hull_mass"], 6e3 - 6e2) - self.assertEqual(self.outputs["platform_cost"], 6 * 2e3) + self.assertEqual(self.outputs["platform_hull_mass"], 6e3 + 1e3 - 6e2) + self.assertEqual(self.outputs["platform_cost"], 6 * 2e3 + 3e3) self.assertEqual(self.outputs["platform_Awater"], 30) self.assertEqual(self.outputs["platform_Iwater"], 6 * 15 + 5 * R.sum()) npt.assert_equal(self.outputs["platform_added_mass"], 6 * np.arange(6)) npt.assert_equal(self.outputs["platform_variable_capacity"], 10 + np.arange(6)) - npt.assert_array_less(1e2, self.outputs["platform_I_total"]) + npt.assert_equal(self.outputs["transition_piece_I"], 1e3 * 0.25 * np.r_[0.5, 0.5, 1.0, np.zeros(3)]) + npt.assert_array_less(1e2, self.outputs["platform_I_hull"] - self.outputs["transition_piece_I"]) # Should find a transition mode even though one wasn't set # npt.assert_equal(self.outputs["transition_node"], [0.0, 0.0, 1.0]) @@ -343,7 +346,7 @@ def testPlatformTower(self): ) / self.outputs["system_mass"], ) - npt.assert_equal(self.outputs["transition_piece_I"], 1e3 * 0.25 * np.r_[0.5, 0.5, 1.0, np.zeros(3)]) + npt.assert_array_less(self.outputs["platform_I_hull"], self.outputs["platform_I_total"]) def testRunFrame(self): myobj = frame.PlatformFrame(options=self.opt) diff --git a/WISDEM/wisdem/test/test_gluecode/test_gc_yaml_floating.py b/WISDEM/wisdem/test/test_gluecode/test_gc_yaml_floating.py new file mode 100644 index 000000000..3a6df969d --- /dev/null +++ b/WISDEM/wisdem/test/test_gluecode/test_gc_yaml_floating.py @@ -0,0 +1,49 @@ +""" +A test where we can modify different parts of the modeling_options to ensure +they work with a top-to-bottom WISDEM run. +""" + +import os +import unittest + +from wisdem.glue_code.runWISDEM import run_wisdem + +test_dir = ( + os.path.dirname(os.path.dirname(os.path.dirname(os.path.dirname(os.path.realpath(__file__))))) + + os.sep + + "examples" + + os.sep + + "09_floating" + + os.sep +) +fname_analysis_options = test_dir + "analysis_options.yaml" +fname_modeling_options = test_dir + "modeling_options.yaml" + + +class TestRegression(unittest.TestCase): + def test15MW(self): + ## IEA 15MW + fname_wt_input = test_dir + "IEA-15-240-RWT_VolturnUS-S.yaml" + wt_opt, modeling_options, opt_options = run_wisdem( + fname_wt_input, fname_modeling_options, fname_analysis_options + ) + + self.assertAlmostEqual(wt_opt["rotorse.re.precomp.blade_mass"][0], 69911.65429172986, 1) + self.assertAlmostEqual(wt_opt["rotorse.rp.AEP"][0] * 1.0e-6, 77.96023669686348, 1) + self.assertAlmostEqual(wt_opt["financese.lcoe"][0] * 1.0e3, 88.65572457443253, 1) + self.assertAlmostEqual(wt_opt["rotorse.rs.tip_pos.tip_deflection"][0], 25.96039967525026, 1) + + +def suite(): + suite = unittest.TestSuite() + suite.addTest(unittest.makeSuite(TestRegression)) + return suite + + +if __name__ == "__main__": + result = unittest.TextTestRunner().run(suite()) + + if result.wasSuccessful(): + exit(0) + else: + exit(1) diff --git a/weis/aeroelasticse/openmdao_openfast.py b/weis/aeroelasticse/openmdao_openfast.py index fe20ca40e..c80d05e0d 100644 --- a/weis/aeroelasticse/openmdao_openfast.py +++ b/weis/aeroelasticse/openmdao_openfast.py @@ -247,7 +247,7 @@ def setup(self): self.add_input("platform_elem_E", NULL * np.ones(NELEM_MAX), units="Pa") self.add_input("platform_elem_G", NULL * np.ones(NELEM_MAX), units="Pa") self.add_discrete_input("platform_elem_memid", [0]*NELEM_MAX) - self.add_input("platform_center_of_mass", np.zeros(3), units="m") + self.add_input("platform_total_center_of_mass", np.zeros(3), units="m") self.add_input("platform_mass", 0.0, units="kg") self.add_input("platform_I_total", np.zeros(6), units="kg*m**2") @@ -709,22 +709,22 @@ def update_FAST_model(self, fst_vt, inputs, discrete_inputs): fst_vt['ElastoDyn']['TipMass(2)'] = 0. fst_vt['ElastoDyn']['TipMass(3)'] = 0. - tower_base_height = max(float(inputs['tower_base_height']), float(inputs["platform_center_of_mass"][2])) + tower_base_height = max(float(inputs['tower_base_height']), float(inputs["platform_total_center_of_mass"][2])) fst_vt['ElastoDyn']['TowerBsHt'] = tower_base_height # Height of tower base above ground level [onshore] or MSL [offshore] (meters) fst_vt['ElastoDyn']['TowerHt'] = float(inputs['hub_height']) - float(inputs['distance_tt_hub']) # Height of tower above ground level [onshore] or MSL [offshore] (meters) # TODO: There is some confusion on PtfmRefzt # DZ: based on the openfast r-tests: - # if this is floating, the z ref. point is 0. Is this the reference that platform_center_of_mass is relative to? + # if this is floating, the z ref. point is 0. Is this the reference that platform_total_center_of_mass is relative to? # if fixed bottom, it's the tower base height. if modopt['flags']['floating']: fst_vt['ElastoDyn']['PtfmMass'] = float(inputs["platform_mass"]) fst_vt['ElastoDyn']['PtfmRIner'] = float(inputs["platform_I_total"][0]) fst_vt['ElastoDyn']['PtfmPIner'] = float(inputs["platform_I_total"][1]) fst_vt['ElastoDyn']['PtfmYIner'] = float(inputs["platform_I_total"][2]) - fst_vt['ElastoDyn']['PtfmCMxt'] = float(inputs["platform_center_of_mass"][0]) - fst_vt['ElastoDyn']['PtfmCMyt'] = float(inputs["platform_center_of_mass"][1]) - fst_vt['ElastoDyn']['PtfmCMzt'] = float(inputs["platform_center_of_mass"][2]) + fst_vt['ElastoDyn']['PtfmCMxt'] = float(inputs["platform_total_center_of_mass"][0]) + fst_vt['ElastoDyn']['PtfmCMyt'] = float(inputs["platform_total_center_of_mass"][1]) + fst_vt['ElastoDyn']['PtfmCMzt'] = float(inputs["platform_total_center_of_mass"][2]) fst_vt['ElastoDyn']['PtfmRefzt'] = 0. # Vertical distance from the ground level [onshore] or MSL [offshore] to the platform reference point (meters) else: diff --git a/weis/frequency/raft_wrapper.py b/weis/frequency/raft_wrapper.py index 243d269b1..632020e2f 100644 --- a/weis/frequency/raft_wrapper.py +++ b/weis/frequency/raft_wrapper.py @@ -124,7 +124,7 @@ def setup(self): self.add_input(f"member{k}:ring_stiffener_web_thickness", 0.0, units="m") self.add_input(f"member{k}:ring_stiffener_flange_width", 1e-6, units="m") self.add_input(f"member{k}:ring_stiffener_flange_thickness", 0.0, units="m") - self.add_input(f"member{k}:ring_stiffener_spacing", 1000.0, units="m") + self.add_input(f"member{k}:ring_stiffener_spacing", 0.0) self.add_input(f"platform_member{k+1}_stations", val=np.zeros(n_height)) self.add_output(f"platform_member{k+1}_heading", val=np.zeros(0), units='deg') @@ -228,8 +228,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs): # Ring stiffener discretization conversion if float(inputs[f"member{k}:ring_stiffener_spacing"]) > 0.0: - outputs[f"platform_member{k+1}_ring_spacing"] = (inputs[f"member{k}:ring_stiffener_spacing"] / - inputs[f"member{k}:height"]) + outputs[f"platform_member{k+1}_ring_spacing"] = inputs[f"member{k}:ring_stiffener_spacing"] h_web = inputs[f"member{k}:ring_stiffener_web_height"] t_web = inputs[f"member{k}:ring_stiffener_web_thickness"] t_flange = inputs[f"member{k}:ring_stiffener_flange_thickness"] diff --git a/weis/glue_code/glue_code.py b/weis/glue_code/glue_code.py index 2e5f8e44a..def79f816 100644 --- a/weis/glue_code/glue_code.py +++ b/weis/glue_code/glue_code.py @@ -384,7 +384,7 @@ def setup(self): self.connect("floatingse.platform_elem_G", "aeroelastic.platform_elem_G") self.connect("floatingse.platform_elem_memid", "aeroelastic.platform_elem_memid") self.connect("floatingse.platform_mass", "aeroelastic.platform_mass") - self.connect("floatingse.platform_center_of_mass", "aeroelastic.platform_center_of_mass") + self.connect("floatingse.platform_total_center_of_mass", "aeroelastic.platform_total_center_of_mass") self.connect("floatingse.platform_I_total", "aeroelastic.platform_I_total") self.connect("floatingse.platform_displacement", "aeroelastic.platform_displacement") self.connect("floating.transition_node", "aeroelastic.transition_node") @@ -416,7 +416,7 @@ def setup(self): self.connect('tower.cd', 'aeroelastic.tower_cd') self.connect('tower_grid.height', 'aeroelastic.tower_height') self.connect('tower_grid.foundation_height', 'aeroelastic.tower_base_height') - #self.connect('monopile.transition_piece_mass', 'aeroelastic.transition_piece_mass') ## TODO + self.connect('floating.transition_piece_mass', 'aeroelastic.transition_piece_mass') self.connect('floatingse.transition_piece_I', 'aeroelastic.transition_piece_I', src_indices=[0,1,2]) self.connect('airfoils.aoa', 'aeroelastic.airfoils_aoa')