diff --git a/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml b/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml index d48d9f0c6..2763e0fac 100644 --- a/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml +++ b/examples/09_floating/IEA-15-240-RWT_VolturnUS-S.yaml @@ -729,10 +729,7 @@ components: line_types: - name: main diameter: 0.185 - mass_density: 685.0 - stiffness: 3270e6 - breaking_load: 1e8 - cost: 100.0 + type: chain transverse_added_mass: 1.0 tangential_added_mass: 0.0 transverse_drag: 1.6 @@ -740,10 +737,7 @@ components: anchor_types: - name: drag_embedment - mass: 1e3 - cost: 1e4 - max_vertical_load: 0.0 - max_lateral_load: 1e5 + type: drag_embedment airfoils: - name: circular diff --git a/examples/09_floating/nrel5mw-semi_oc4.yaml b/examples/09_floating/nrel5mw-semi_oc4.yaml index aa030486d..ca8055e08 100644 --- a/examples/09_floating/nrel5mw-semi_oc4.yaml +++ b/examples/09_floating/nrel5mw-semi_oc4.yaml @@ -634,10 +634,7 @@ components: line_types: - name: main diameter: 0.0766 - mass_density: 113.35 - stiffness: 7.536e8 - breaking_load: 1e8 - cost: 100.0 + type: chain transverse_added_mass: 0.8 tangential_added_mass: 0.25 transverse_drag: 2.0 @@ -645,10 +642,7 @@ components: anchor_types: - name: drag_embedment - mass: 1e3 - cost: 1e4 - max_vertical_load: 0.0 - max_lateral_load: 1e5 + type: drag_embedment airfoils: - aerodynamic_center: 0.275 diff --git a/examples/09_floating/nrel5mw-spar_oc3.yaml b/examples/09_floating/nrel5mw-spar_oc3.yaml index 5b95b4218..0dd230411 100644 --- a/examples/09_floating/nrel5mw-spar_oc3.yaml +++ b/examples/09_floating/nrel5mw-spar_oc3.yaml @@ -405,10 +405,7 @@ components: line_types: - name: main diameter: 0.09 - mass_density: 77.7066 - stiffness: 384.243e6 - breaking_load: 1e8 - cost: 100.0 + type: chain transverse_added_mass: 1.0 tangential_added_mass: 0.0 transverse_drag: 1.6 @@ -416,10 +413,7 @@ components: anchor_types: - name: drag_embedment - mass: 1e3 - cost: 1e4 - max_vertical_load: 0.0 - max_lateral_load: 1e5 + type: drag_embedment airfoils: - aerodynamic_center: 0.275 @@ -709,11 +703,14 @@ control: maxOmega: 1.26711 pitch: PC_zeta: 0.7 # Pitch controller desired damping ratio [-] - PC_omega: 0.5 # Pitch controller desired natural frequency [rad/s] - ps_percent: 0.8 # Percent peak shaving [%, <= 1 ], {default = 80%} + PC_omega: 0.35 # Pitch controller desired natural frequency [rad/s] + ps_percent: 0.9 # Percent peak shaving [%, <= 1 ], {default = 80%} max_pitch: 1.57 # Maximum pitch angle [rad], {default = 90 degrees} max_pitch_rate: 0.1745 min_pitch: 0 # Minimum pitch angle [rad], {default = 0 degrees} + fl_feedback: + twr_freq: 100 # Frequency of notch in floating feedback (F_NotchCornerFreq) + ptfm_freq: .2 # Frequency of LPF in floating feedback (F_FlCornerFreq) torque: control_type: tsr_tracking tsr: 7.01754386 @@ -738,8 +735,8 @@ environment: water_density: 1025.0 water_dyn_viscosity: 1.3351e-3 water_depth: 320.0 - significant_wave_height: 9.14 - significant_wave_period: 13.6 + significant_wave_height: 6.0 + significant_wave_period: 10.0 soil_shear_modulus: 140.e+6 soil_poisson: 0.4 diff --git a/examples/09_floating/tlp_example.py b/examples/09_floating/tlp_example.py index 76c30a56f..8416ba662 100644 --- a/examples/09_floating/tlp_example.py +++ b/examples/09_floating/tlp_example.py @@ -41,6 +41,8 @@ opt["mooring"] = {} opt["mooring"]["n_attach"] = 3 opt["mooring"]["n_anchors"] = 3 +opt["mooring"]["line_anchor"] = ["custom"] * 3 +opt["mooring"]["line_material"] = ["custom"] * 3 opt["materials"] = {} opt["materials"]["n_mat"] = 2 diff --git a/wisdem/commonse/frustum.py b/wisdem/commonse/frustum.py index cf29f1df8..192a66a32 100644 --- a/wisdem/commonse/frustum.py +++ b/wisdem/commonse/frustum.py @@ -34,7 +34,7 @@ def frustum(Db, Dt, H): return vol, cm -def frustumVol(rb, rt, h, diamFlag=False): +def frustumVol(rb_0, rt_0, h, diamFlag=False): """This function returns a frustum's volume with radii or diameter inputs. INPUTS: @@ -51,12 +51,13 @@ def frustumVol(rb, rt, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 return np.pi * (h / 3.0) * (rb * rb + rt * rt + rb * rt) -def frustumCG(rb, rt, h, diamFlag=False): +def frustumCG(rb_0, rt_0, h, diamFlag=False): """This function returns a frustum's center of mass/gravity (centroid) with radii or diameter inputs. NOTE: This is for a SOLID frustum, not a shell @@ -74,12 +75,13 @@ def frustumCG(rb, rt, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 return 0.25 * h * (rb ** 2 + 2.0 * rb * rt + 3.0 * rt ** 2) / (rb ** 2 + rb * rt + rt ** 2) -def frustumIzz(rb, rt, h, diamFlag=False): +def frustumIzz(rb_0, rt_0, h, diamFlag=False): """This function returns a frustum's mass-moment of inertia (divided by density) about the central (axial) z-axis with radii or diameter inputs. NOTE: This is for a SOLID frustum, not a shell @@ -98,15 +100,16 @@ def frustumIzz(rb, rt, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 # Integrate 2*pi*r*r^2 dr dz from r=0 to r(z), z=0 to h # Also equals 0.3*Vol * (rt**5.0 - rb**5.0) / (rt**3.0 - rb**3.0) # Also equals (0.1*np.pi*h * (rt**5.0 - rb**5.0) / (rt - rb) ) return 0.1 * np.pi * h * (rt ** 4.0 + rb * rt ** 3 + rb ** 2 * rt ** 2 + rb ** 3 * rt + rb ** 4.0) -def frustumIxx(rb, rt, h, diamFlag=False): +def frustumIxx(rb_0, rt_0, h, diamFlag=False): """This function returns a frustum's mass-moment of inertia (divided by density) about the transverse x/y-axis passing through the center of mass with radii or diameter inputs. NOTE: This is for a SOLID frustum, not a shell @@ -125,10 +128,11 @@ def frustumIxx(rb, rt, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 # Integrate pi*r(z)^4/4 + pi*r(z)^2*(z-z_cg)^2 dz from z=0 to h - A = 0.5 * frustumIzz(rb, rt, h) + A = 0.5 * frustumIzz(rb_0, rt_0, h) B = ( np.pi * h ** 3 @@ -141,7 +145,7 @@ def frustumIxx(rb, rt, h, diamFlag=False): return A + B -def frustumShellVol(rb, rt, t, h, diamFlag=False): +def frustumShellVol(rb_0, rt_0, t, h, diamFlag=False): """This function returns a frustum shell's volume (for computing mass with density) with radii or diameter inputs. NOTE: This is for a frustum SHELL, not a solid @@ -160,8 +164,9 @@ def frustumShellVol(rb, rt, t, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 # Integrate 2*pi*r*dr*dz from r=ri(z) to ro(z), z=0 to h rb_o = rb rb_i = rb - t @@ -171,7 +176,7 @@ def frustumShellVol(rb, rt, t, h, diamFlag=False): return frustumVol(rb_o, rt_o, h) - frustumVol(rb_i, rt_i, h) -def frustumShellCG(rb, rt, t, h, diamFlag=False): +def frustumShellCG(rb_0, rt_0, t, h, diamFlag=False): """This function returns a frustum's center of mass/gravity (centroid) with radii or diameter inputs. NOTE: This is for a frustum SHELL, not a solid @@ -190,8 +195,9 @@ def frustumShellCG(rb, rt, t, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 # Integrate 2*pi*r*z*dr*dz/V from r=ri(z) to ro(z), z=0 to h rb_o = rb rb_i = rb - t @@ -202,7 +208,7 @@ def frustumShellCG(rb, rt, t, h, diamFlag=False): return h * A / 4.0 / B -def frustumShellIzz(rb, rt, t, h, diamFlag=False): +def frustumShellIzz(rb_0, rt_0, t, h, diamFlag=False): """This function returns a frustum's mass-moment of inertia (divided by density) about the central (axial) z-axis with radii or diameter inputs. NOTE: This is for a frustum SHELL, not a solid @@ -222,8 +228,9 @@ def frustumShellIzz(rb, rt, t, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 # Integrate 2*pi*r*dr*dz from r=ri(z) to ro(z), z=0 to h rb_o = rb rb_i = rb - t @@ -232,7 +239,7 @@ def frustumShellIzz(rb, rt, t, h, diamFlag=False): return frustumIzz(rb_o, rt_o, h) - frustumIzz(rb_i, rt_i, h) -def frustumShellIxx(rb, rt, t, h, diamFlag=False): +def frustumShellIxx(rb_0, rt_0, t, h, diamFlag=False): """This function returns a frustum's mass-moment of inertia (divided by density) about the transverse x/y-axis passing through the center of mass with radii or diameter inputs. NOTE: This is for a frustum SHELL, not a solid @@ -252,8 +259,9 @@ def frustumShellIxx(rb, rt, t, h, diamFlag=False): """ if diamFlag: # Convert diameters to radii - rb *= 0.5 - rt *= 0.5 + rb, rt = 0.5 * rb_0, 0.5 * rt_0 + else: + rb, rt = rb_0, rt_0 # Integrate 2*pi*r*dr*dz from r=ri(z) to ro(z), z=0 to h rb_o = rb rb_i = rb - t diff --git a/wisdem/floatingse/floating.py b/wisdem/floatingse/floating.py index a579dcd07..1bd6db57c 100644 --- a/wisdem/floatingse/floating.py +++ b/wisdem/floatingse/floating.py @@ -110,8 +110,3 @@ def setup(self): for k in range(n_member): for var in mem_vars: self.connect(f"member{k}." + var, f"member{k}:" + var) - - """ - self.connect("max_offset_restoring_force", "mooring_surge_restoring_force") - self.connect("operational_heel_restoring_force", "mooring_pitch_restoring_force") - """ diff --git a/wisdem/floatingse/floating_frame.py b/wisdem/floatingse/floating_frame.py index c971d97b5..b12ad99b6 100644 --- a/wisdem/floatingse/floating_frame.py +++ b/wisdem/floatingse/floating_frame.py @@ -470,6 +470,7 @@ def setup(self): self.add_output("system_center_of_mass", np.zeros(3), units="m") self.add_output("system_mass", 0.0, units="kg") self.add_output("variable_ballast_mass", 0.0, units="kg") + self.add_output("variable_center_of_mass", val=np.zeros(3), units="m") 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)) @@ -619,6 +620,7 @@ def compute(self, inputs, outputs): outputs["member_variable_height"][k] = s_end - spts[0] cg_variable = np.dot(V_variable_member, cg_variable_member) / V_variable + outputs["variable_center_of_mass"] = cg_variable # Now find total system mass outputs["system_mass"] = m_sys + m_variable @@ -643,6 +645,9 @@ def setup(self): 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_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") self.add_input("tower_nodes", NULL * np.ones((MEMMAX, 3)), units="m") self.add_input("tower_Fnode", NULL * np.ones((MEMMAX, 3)), units="N") @@ -702,6 +707,9 @@ def setup(self): self.add_input("rna_I", np.zeros(6), units="kg*m**2") self.add_input("mooring_neutral_load", np.zeros((n_attach, 3)), units="N") self.add_input("mooring_fairlead_joints", np.zeros((n_attach, 3)), units="m") + self.add_input("mooring_stiffness", np.zeros((6, 6)), units="N/m") + self.add_input("variable_ballast_mass", 0.0, units="kg") + self.add_input("variable_center_of_mass", val=np.zeros(3), units="m") NFREQ2 = int(NFREQ / 2) self.add_output("tower_freqs", val=np.zeros(NFREQ), units="Hz") @@ -737,11 +745,16 @@ def compute(self, inputs, outputs): n_attach = opt["mooring"]["n_attach"] m_rna = float(inputs["rna_mass"]) cg_rna = inputs["rna_cg"] + cb = inputs["platform_center_of_buoyancy"] + V_total = inputs["platform_displacement"] I_rna = inputs["rna_I"] I_trans = inputs["transition_piece_I"] + m_variable = float(inputs["variable_ballast_mass"]) + cg_variable = inputs["variable_center_of_mass"] fairlead_joints = inputs["mooring_fairlead_joints"] mooringF = inputs["mooring_neutral_load"] + mooringK = np.abs(np.diag(inputs["mooring_stiffness"])) # Create frame3dd instance: nodes, elements, reactions, and options for frame in ["tower", "system"]: @@ -776,9 +789,20 @@ def compute(self, inputs, outputs): ielem = np.arange(nelem) + 1 elem_obj = pyframe3dd.ElementData(ielem, N1 + 1, N2 + 1, A, Asx, Asy, Izz, Ixx, Iyy, E, G, roll, rho) - # TODO: Hydro_K + Mooring_K for tower (system too?) - rid = np.array([itrans]) # np.array([np.argmin(nodes[:, 2])]) - Rx = Ry = Rz = Rxx = Ryy = Rzz = np.array([RIGID]) + # Use Mooring stiffness (TODO Hydro_K too) + if frame == "tower": + rid = np.array([itrans]) # np.array([np.argmin(nodes[:, 2])]) + else: + ind = [] + for k in range(n_attach): + ind.append(util.closest_node(nodes, fairlead_joints[k, :])) + rid = np.array([ind]) # np.array([np.argmin(nodes[:, 2])]) + + Rx = Ry = Rz = Rxx = Ryy = Rzz = RIGID * np.ones(rid.size) + # Rx, Ry, Rz = [mooringK[0]], [mooringK[1]], [mooringK[2]] + # Only this solution works and there isn't much different with fully rigid + # Rx, Ry, Rz = [RIGID], [RIGID], [mooringK[2]] + # Rxx, Ryy, Rzz = [RIGID], [RIGID], [RIGID] react_obj = pyframe3dd.ReactionData(rid + 1, Rx, Ry, Rz, Rxx, Ryy, Rzz, rigid=RIGID) frame3dd_opt = opt["WISDEM"]["FloatingSE"]["frame3dd"] @@ -789,10 +813,11 @@ def compute(self, inputs, outputs): # Added mass m_trans = float(inputs["transition_piece_mass"]) if frame == "tower": - # TODO: Added mass and stiffness - m_trans += float(inputs["platform_mass"]) + m_trans += float(inputs["platform_mass"]) + inputs["platform_added_mass"][0] + m_variable cg_trans = inputs["transition_node"] - inputs["platform_center_of_mass"] + I_trans[:3] += inputs["platform_added_mass"][3:] else: + m_trans += m_variable cg_trans = np.zeros(3) add_gravity = True mID = np.array([itrans, ihub], dtype=np.int_).flatten() @@ -830,6 +855,11 @@ def compute(self, inputs, outputs): for k in range(n_attach): ind = util.closest_node(nodes, fairlead_joints[k, :]) Fnode[ind, :] += mooringF[k, :] + else: + # Combine all buoyancy forces into one + ind = util.closest_node(nodes, cb) + Fnode[ind, -1] += V_total * 1025 * gravity + Fnode[ihub, :] += inputs["rna_F"] Mnode[ihub, :] += inputs["rna_M"] nF = np.where(np.abs(Fnode).sum(axis=1) > 0.0)[0] diff --git a/wisdem/floatingse/member.py b/wisdem/floatingse/member.py index d0dd1fbf3..feaa4e90e 100644 --- a/wisdem/floatingse/member.py +++ b/wisdem/floatingse/member.py @@ -1651,7 +1651,9 @@ def setup(self): ) self.add_subsystem( - "gc", GeometricConstraints(nPoints=n_height, diamFlag=True), promotes=["constr_taper", "constr_d_to_t"] + "gc", + GeometricConstraints(nPoints=n_height, diamFlag=True), + promotes=["constr_taper", "constr_d_to_t", "slope"], ) self.connect("outer_diameter", "gc.d") self.connect("wall_thickness", "gc.t") diff --git a/wisdem/floatingse/mooring.py b/wisdem/floatingse/mooring.py index e2c7a5330..61e4e1bea 100644 --- a/wisdem/floatingse/mooring.py +++ b/wisdem/floatingse/mooring.py @@ -1,6 +1,7 @@ import numpy as np import openmdao.api as om import wisdem.moorpy as mp +import wisdem.moorpy.MoorProps as props NLINES_MAX = 15 NPTS_PLOT = 101 @@ -25,12 +26,30 @@ class Mooring(om.ExplicitComponent): Depth below water for mooring line attachment line_length : float, [m] Unstretched total mooring line length - anchor_radius : float, [m] - radius from center of spar to mooring anchor point line_diameter : float, [m] diameter of mooring line + line_type : string + Type of mooring line: chain, chain_stud, nylon, polyester, polypropylene, fiber, wire, iwrc, or custom + line_mass_density_coeff: float, [kg/m/m^2] + Line mass density per unit length per diameter^2, if line type is custom + line_stiffness_coeff: float, [N/m^2] + Line stiffness (E*A) per diameter^2, if line type is custom + line_breaking_load_coeff: float, [N/m^2] + Line minimumum breaking load (MBL) per diameter^2, if line type is custom + line_cost_rate_coeff: float, [USD/m/m^2] + Line cost per unit length per diameter^2, if line type is custom + anchor_radius : float, [m] + radius from center of spar to mooring anchor point anchor_type : string - SUCTIONPILE or DRAGEMBEDMENT + Type of anchor for sizing: drag_embedment, suction, plate, micropile, sepla, or custom + anchor_mass : float, [USD] + mass for one anchor, if anchor type is custom + anchor_cost : float, [USD] + cost for one anchor, if anchor type is custom + anchor_max_vertical_load : float, [N] + Maximum tolerable vertical force on anchor, if anchor type is custom + anchor_max_lateral_load : float, [N] + Maximum tolerable lateral force on anchor, if anchor type is custom max_surge_fraction : float Maximum allowable surge offset as a fraction of water depth (0-1) operational_heel : float, [deg] @@ -48,8 +67,6 @@ class Mooring(om.ExplicitComponent): total cost for anchor + legs + miscellaneous costs mooring_stiffness : numpy array[6, 6], [N/m] Linearized stiffness matrix of mooring system at neutral (no offset) conditions. - anchor_cost : float, [USD] - total cost for anchor mooring_neutral_load : numpy array[NLINES_MAX, 3], [N] mooring vertical load in all mooring lines max_offset_restoring_force : float, [N] @@ -60,10 +77,14 @@ class Mooring(om.ExplicitComponent): forces for all mooring lines after max survival heel mooring_plot_matrix : numpy array[NLINES_MAX, NPTS_PLOT, 3], [m] data matrix for plotting - constr_axial_load : float, [m] + constr_axial_load : float range of damaged mooring constr_mooring_length : float mooring line length ratio to nodal distance + constr_anchor_vertical : numpy array[n_lines] + Maximum allowable vertical anchor force minus vertical line tension times safety factor (must be >= 0) + constr_anchor_lateral : numpy array[n_lines] + Maximum allowable lateral anchor force minus lateral line tension times safety factor (must be >= 0) """ @@ -86,8 +107,12 @@ def setup(self): self.add_input("fairlead", 0.0, units="m") self.add_input("line_length", 0.0, units="m") self.add_input("line_diameter", 0.0, units="m") + self.add_input("anchor_radius", 0.0, units="m") + self.add_input("anchor_mass", 0.0, units="kg") self.add_input("anchor_cost", 0.0, units="USD") + self.add_input("anchor_max_vertical_load", 1e30, units="N") + self.add_input("anchor_max_lateral_load", 1e30, units="N") self.add_input("line_mass_density_coeff", 0.0, units="kg/m**3") self.add_input("line_stiffness_coeff", 0.0, units="N/m**2") @@ -95,8 +120,6 @@ def setup(self): self.add_input("line_cost_rate_coeff", 0.0, units="USD/m**3") # User inputs (could be design variables) - # self.add_discrete_input("mooring_type", "CHAIN") - # self.add_discrete_input("anchor_type", "DRAGEMBEDMENT") self.add_input("max_surge_fraction", 0.1) self.add_input("operational_heel", 0.0, units="rad") self.add_input("survival_heel", 0.0, units="rad") @@ -112,12 +135,12 @@ def setup(self): self.add_output("mooring_plot_matrix", np.zeros((n_lines, NPTS_PLOT, 3)), units="m") # Constraints - self.add_output("constr_axial_load", 0.0, units="m") + self.add_output("constr_axial_load", 0.0) self.add_output("constr_mooring_length", 0.0) + self.add_output("constr_anchor_vertical", np.zeros(n_lines)) + self.add_output("constr_anchor_lateral", np.zeros(n_lines)) def compute(self, inputs, outputs): - # Set characteristics based on regressions / empirical data - # self.set_properties(inputs, discrete_inputs) # Write MAP input file and analyze the system at every angle self.evaluate_mooring(inputs, outputs) @@ -133,16 +156,32 @@ def evaluate_mooring(self, inputs, outputs): R_anchor = float(inputs["anchor_radius"]) heel = float(inputs["operational_heel"]) max_heel = float(inputs["survival_heel"]) - d = inputs["line_diameter"] + d = float(inputs["line_diameter"]) L_mooring = inputs["line_length"] - min_break_load = inputs["line_breaking_load_coeff"] * d ** 2 gamma = self.options["gamma"] n_attach = self.options["options"]["n_attach"] n_lines = self.options["options"]["n_anchors"] offset = float(inputs["max_surge_fraction"]) * water_depth n_anchors = self.options["options"]["n_anchors"] ratio = int(n_anchors / n_attach) - gamma = self.options["gamma"] + + line_obj = None + line_mat = self.options["options"]["line_material"][0] + if line_mat == "custom": + min_break_load = float(inputs["line_breaking_load_coeff"]) * d ** 2 + mass_den = float(inputs["line_mass_density_coeff"]) * d ** 2 + ea_stiff = float(inputs["line_stiffness_coeff"]) * d ** 2 + cost_rate = float(inputs["line_cost_rate_coeff"]) * d ** 2 + elif line_mat == "chain_stud": + line_obj = props.getLineProps(1e3 * d, type="chain", stud="stud") + else: + line_obj = props.getLineProps(1e3 * d, type=line_mat) + + if not line_obj is None: + min_break_load = line_obj.MBL + mass_den = line_obj.mlin + ea_stiff = line_obj.EA + cost_rate = line_obj.cost # Geometric constraints on line length if L_mooring > (water_depth - fairlead_depth): @@ -191,10 +230,10 @@ def evaluate_mooring(self, inputs, outputs): config["line_types"] = [{}] config["line_types"][0]["name"] = "myline" config["line_types"][0]["diameter"] = d - config["line_types"][0]["mass_density"] = inputs["line_mass_density_coeff"] * d ** 2 - config["line_types"][0]["stiffness"] = inputs["line_stiffness_coeff"] * d ** 2 - config["line_types"][0]["breaking_load"] = inputs["line_breaking_load_coeff"] * d ** 2 - config["line_types"][0]["cost"] = inputs["line_cost_rate_coeff"] * d ** 2 + config["line_types"][0]["mass_density"] = mass_den + config["line_types"][0]["stiffness"] = ea_stiff + config["line_types"][0]["breaking_load"] = min_break_load + config["line_types"][0]["cost"] = cost_rate config["line_types"][0]["transverse_added_mass"] = 0.0 config["line_types"][0]["tangential_added_mass"] = 0.0 config["line_types"][0]["transverse_drag"] = 0.0 @@ -237,7 +276,15 @@ def evaluate_mooring(self, inputs, outputs): F_maxheel = ms.mooringEq([0, 0, 0, 0, max_heel, 0], DOFtype="coupled") outputs["survival_heel_restoring_force"] = F_maxheel - # TODO: Vertical loads on anchors? + # Anchor load limits + F_anch = np.zeros((n_anchors, 3)) + for k in range(n_anchors): + if np.abs(ms.pointList[k + n_attach].r[0] + water_depth) < 0.1: + F_anch[k, :] = ms.pointList[k].getForces() + outputs["constr_anchor_lateral"] = ( + inputs["anchor_max_lateral_load"] - np.sqrt(np.sum(F_anch[:, :-1] ** 2, axis=1)) * gamma + ) + outputs["constr_anchor_vertical"] = inputs["anchor_max_vertical_load"] - F_anch[:, -1] * gamma # Get angles by which to find the weakest line dangle = 5.0 @@ -274,27 +321,42 @@ def evaluate_mooring(self, inputs, outputs): def compute_cost(self, inputs, outputs): # Unpack variables L_mooring = float(inputs["line_length"]) - # anchorType = discrete_inputs["anchor_type"] d = float(inputs["line_diameter"]) - cost_per_length = float(inputs["line_cost_rate_coeff"]) * d ** 2 - # min_break_load = inputs['line_breaking_load_coeff'] * d**2 - wet_mass_per_length = float(inputs["line_mass_density_coeff"]) * d ** 2 - anchor_rate = float(inputs["anchor_cost"]) + gamma = self.options["gamma"] + + anchor_type = self.options["options"]["line_anchor"][0] + if anchor_type == "custom": + anchor_rate = float(inputs["anchor_cost"]) + anchor_mass = float(inputs["anchor_mass"]) + else: + # Do empirical sizing with MoorPy + fx = (inputs["anchor_max_lateral_load"] - outputs["constr_anchor_lateral"].min()) / gamma + fz = (inputs["anchor_max_vertical_load"] - outputs["constr_anchor_vertical"].min()) / gamma + anchor_rate, _, _ = props.getAnchorProps(fx, fz, type=anchor_type.replace("_", "-")) + anchor_mass = 0.0 # TODO n_anchors = n_lines = self.options["options"]["n_anchors"] + line_obj = None + line_mat = self.options["options"]["line_material"][0] + if line_mat == "custom": + mass_den = float(inputs["line_mass_density_coeff"]) * d ** 2 + cost_rate = float(inputs["line_cost_rate_coeff"]) * d ** 2 + elif line_mat == "chain_stud": + line_obj = props.getLineProps(1e3 * d, type="chain", stud="stud") + else: + line_obj = props.getLineProps(1e3 * d, type=line_mat) + + if not line_obj is None: + mass_den = line_obj.mlin + cost_rate = line_obj.cost + # Cost of anchors - # if anchorType.upper() == "DRAGEMBEDMENT": - # anchor_rate = 1e-3 * min_break_load / gravity / 20 * 2000 - # elif anchorType.upper() == "SUCTIONPILE": - # anchor_rate = 150000.0 * np.sqrt(1e-3 * min_break_load / gravity / 1250.0) - # else: - # raise ValueError("Anchor Type must be DRAGEMBEDMENT or SUCTIONPILE") anchor_total = anchor_rate * n_anchors # Cost of all of the mooring lines - legs_total = n_lines * cost_per_length * L_mooring + legs_total = n_lines * cost_rate * L_mooring # Total summations outputs["mooring_cost"] = legs_total + anchor_total - outputs["line_mass"] = wet_mass_per_length * L_mooring - outputs["mooring_mass"] = wet_mass_per_length * L_mooring * n_lines + outputs["line_mass"] = mass_den * L_mooring + outputs["mooring_mass"] = (outputs["line_mass"] + anchor_mass) * n_lines diff --git a/wisdem/glue_code/gc_LoadInputs.py b/wisdem/glue_code/gc_LoadInputs.py index 34bb96465..d75a36a04 100644 --- a/wisdem/glue_code/gc_LoadInputs.py +++ b/wisdem/glue_code/gc_LoadInputs.py @@ -499,6 +499,8 @@ def set_openmdao_vectors(self): self.modeling_options["mooring"]["node1"] = [""] * n_lines self.modeling_options["mooring"]["node2"] = [""] * n_lines self.modeling_options["mooring"]["line_type"] = [""] * n_lines + self.modeling_options["mooring"]["line_material"] = [""] * n_lines + self.modeling_options["mooring"]["line_anchor"] = [""] * n_lines fairlead_nodes = [] for i in range(n_lines): self.modeling_options["mooring"]["node1"][i] = self.wt_init["components"]["mooring"]["lines"][i][ @@ -521,17 +523,50 @@ def set_openmdao_vectors(self): fairlead_nodes.append(self.wt_init["components"]["mooring"]["nodes"][node1id]["joint"]) if self.modeling_options["mooring"]["node_type"][node2id] == "vessel": fairlead_nodes.append(self.wt_init["components"]["mooring"]["nodes"][node2id]["joint"]) + # Store the anchor type names to start + if self.modeling_options["mooring"]["node_type"][node1id] == "fixed": + self.modeling_options["mooring"]["line_anchor"][i] = self.modeling_options["mooring"][ + "anchor_type" + ][node1id] + if self.modeling_options["mooring"]["node_type"][node2id] == "fixed": + self.modeling_options["mooring"]["line_anchor"][i] = self.modeling_options["mooring"][ + "anchor_type" + ][node2id] self.modeling_options["mooring"]["line_type_name"] = [""] * n_line_types + self.modeling_options["mooring"]["line_type_type"] = [""] * n_line_types for i in range(n_line_types): self.modeling_options["mooring"]["line_type_name"][i] = self.wt_init["components"]["mooring"][ "line_types" ][i]["name"] + self.modeling_options["mooring"]["line_type_type"][i] = self.wt_init["components"]["mooring"][ + "line_types" + ][i]["type"].lower() + for j in range(n_lines): + if ( + self.modeling_options["mooring"]["line_type"][j] + == self.modeling_options["mooring"]["line_type_name"][i] + ): + self.modeling_options["mooring"]["line_material"][j] = self.modeling_options["mooring"][ + "line_type_type" + ][i] self.modeling_options["mooring"]["anchor_type_name"] = [""] * n_anchor_types + self.modeling_options["mooring"]["anchor_type_type"] = [""] * n_anchor_types for i in range(n_anchor_types): self.modeling_options["mooring"]["anchor_type_name"][i] = self.wt_init["components"]["mooring"][ "anchor_types" ][i]["name"] + self.modeling_options["mooring"]["anchor_type_type"][i] = self.wt_init["components"]["mooring"][ + "anchor_types" + ][i]["type"].lower() + for j in range(n_lines): + if ( + self.modeling_options["mooring"]["line_anchor"][j] + == self.modeling_options["mooring"]["anchor_type_name"][i] + ): + self.modeling_options["mooring"]["line_anchor"][j] = self.modeling_options["mooring"][ + "anchor_type_type" + ][i] self.modeling_options["mooring"]["n_attach"] = len(set(fairlead_nodes)) # Assembly diff --git a/wisdem/glue_code/gc_PoseOptimization.py b/wisdem/glue_code/gc_PoseOptimization.py index f21600b3c..2ddcfb902 100644 --- a/wisdem/glue_code/gc_PoseOptimization.py +++ b/wisdem/glue_code/gc_PoseOptimization.py @@ -409,13 +409,19 @@ def set_objective(self, wt_opt): wt_opt.model.add_objective("tcons.tip_deflection_ratio") elif self.opt["merit_figure"] == "tower_mass": - wt_opt.model.add_objective("towerse.tower_mass", ref=1e6) + if not self.modeling["flags"]["floating"]: + wt_opt.model.add_objective("towerse.tower_mass", ref=1e6) + 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"] == "structural_mass": - wt_opt.model.add_objective("towerse.structural_mass", ref=1e6) + if not self.modeling["flags"]["floating"]: + wt_opt.model.add_objective("towerse.structural_mass", ref=1e6) + else: + wt_opt.model.add_objective("floatingse.system_structural_mass", ref=1e6) elif self.opt["merit_figure"] == "tower_cost": wt_opt.model.add_objective("tcc.tower_cost", ref=1e6) @@ -971,17 +977,30 @@ def set_constraints(self, wt_opt): wt_opt.model.add_constraint("floatingse.constr_tower_shell_buckling", upper=1.0) if tower_constr["d_to_t"]["flag"] or monopile_constr["d_to_t"]["flag"]: - wt_opt.model.add_constraint( - "towerse.constr_d_to_t", - lower=tower_constr["d_to_t"]["lower_bound"], - upper=tower_constr["d_to_t"]["upper_bound"], - ) + if not self.modeling["flags"]["floating"]: + wt_opt.model.add_constraint( + "towerse.constr_d_to_t", + lower=tower_constr["d_to_t"]["lower_bound"], + upper=tower_constr["d_to_t"]["upper_bound"], + ) + else: + wt_opt.model.add_constraint( + "floatingse.tower.constr_d_to_t", + lower=tower_constr["d_to_t"]["lower_bound"], + upper=tower_constr["d_to_t"]["upper_bound"], + ) if tower_constr["taper"]["flag"] or monopile_constr["taper"]["flag"]: - wt_opt.model.add_constraint("towerse.constr_taper", lower=tower_constr["taper"]["lower_bound"]) + if not self.modeling["flags"]["floating"]: + wt_opt.model.add_constraint("towerse.constr_taper", lower=tower_constr["taper"]["lower_bound"]) + else: + wt_opt.model.add_constraint("floatingse.tower.constr_taper", lower=tower_constr["taper"]["lower_bound"]) if tower_constr["slope"]["flag"] or monopile_constr["slope"]["flag"]: - wt_opt.model.add_constraint("towerse.slope", upper=1.0) + if not self.modeling["flags"]["floating"]: + wt_opt.model.add_constraint("towerse.slope", upper=1.0) + else: + wt_opt.model.add_constraint("floatingse.tower.slope", upper=1.0) if monopile_constr["pile_depth"]["flag"]: wt_opt.model.add_constraint("towerse.suctionpile_depth", lower=monopile_constr["pile_depth"]["lower_bound"]) @@ -991,10 +1010,18 @@ def set_constraints(self, wt_opt): wt_opt.model.add_constraint("tcons.constr_tower_f_NPmargin", upper=0.0) elif tower_constr["frequency_1"]["flag"] or monopile_constr["frequency_1"]["flag"]: - for k in range(self.modeling["WISDEM"]["TowerSE"]["nLC"]): - kstr = "" if self.modeling["WISDEM"]["TowerSE"]["nLC"] <= 1 else str(k + 1) + if not self.modeling["flags"]["floating"]: + for k in range(self.modeling["WISDEM"]["TowerSE"]["nLC"]): + kstr = "" if self.modeling["WISDEM"]["TowerSE"]["nLC"] <= 1 else str(k + 1) + wt_opt.model.add_constraint( + "towerse.tower" + kstr + ".structural_frequencies", + indices=[0], + lower=tower_constr["frequency_1"]["lower_bound"], + upper=tower_constr["frequency_1"]["upper_bound"], + ) + else: wt_opt.model.add_constraint( - "towerse.tower" + kstr + ".structural_frequencies", + "floatingse.tower_freqs", indices=[0], lower=tower_constr["frequency_1"]["lower_bound"], upper=tower_constr["frequency_1"]["upper_bound"], @@ -1052,6 +1079,12 @@ def set_constraints(self, wt_opt): if float_constr["mooring_length"]["flag"]: wt_opt.model.add_constraint("floatingse.constr_mooring_length", upper=1.0) + if float_constr["anchor_vertical"]["flag"]: + wt_opt.model.add_constraint("floatingse.constr_anchor_vertical", lower=0.0) + + if float_constr["anchor_lateral"]["flag"]: + wt_opt.model.add_constraint("floatingse.constr_anchor_lateral", lower=0.0) + if float_constr["stress"]["flag"]: wt_opt.model.add_constraint("floatingse.constr_system_stress", upper=1.0) diff --git a/wisdem/glue_code/gc_WT_DataStruc.py b/wisdem/glue_code/gc_WT_DataStruc.py index 9becdabac..8500680fc 100644 --- a/wisdem/glue_code/gc_WT_DataStruc.py +++ b/wisdem/glue_code/gc_WT_DataStruc.py @@ -1,6 +1,8 @@ import copy + import numpy as np 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 @@ -77,19 +79,17 @@ def setup(self): if modeling_options["WISDEM"]["RotorSE"]["inn_af"]: inn_af = om.IndepVarComp() inn_af.add_output( - "s_opt_r_thick", - val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["t/c"]["n_opt"]) + "s_opt_r_thick", val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["t/c"]["n_opt"]) ) inn_af.add_output( - "r_thick_opt", + "r_thick_opt", val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["t/c"]["n_opt"]), ) inn_af.add_output( - "s_opt_L_D", - val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["L/D"]["n_opt"]) + "s_opt_L_D", val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["L/D"]["n_opt"]) ) inn_af.add_output( - "L_D_opt", + "L_D_opt", val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["L/D"]["n_opt"]), ) self.add_subsystem("inn_af", inn_af) @@ -133,7 +133,7 @@ def setup(self): self.connect("airfoils.cd", "blade.interp_airfoils.cd") self.connect("airfoils.cm", "blade.interp_airfoils.cm") - if modeling_options["WISDEM"]["RotorSE"]["inn_af"]: + if modeling_options["WISDEM"]["RotorSE"]["inn_af"]: self.connect("airfoils.aoa", "blade.run_inn_af.aoa") self.connect("inn_af.s_opt_r_thick", "blade.run_inn_af.s_opt_r_thick") @@ -161,7 +161,10 @@ def setup(self): nacelle_ivc.add_output("gear_ratio", val=1.0, desc="Total gear ratio of drivetrain (use 1.0 for direct)") if modeling_options["flags"]["nacelle"]: nacelle_ivc.add_output( - "distance_hub2mb", val=0.0, units="m", desc="Distance from hub flange to first main bearing along shaft" + "distance_hub2mb", + val=0.0, + units="m", + desc="Distance from hub flange to first main bearing along shaft", ) nacelle_ivc.add_output( "distance_mb2mb", val=0.0, units="m", desc="Distance from first to second main bearing along shaft" @@ -196,8 +199,12 @@ def setup(self): units="kg", desc="Override regular regression-based calculation of transformer mass with this value", ) - nacelle_ivc.add_discrete_output("mb1Type", val="CARB", desc="Type of main bearing: CARB / CRB / SRB / TRB") - nacelle_ivc.add_discrete_output("mb2Type", val="SRB", desc="Type of main bearing: CARB / CRB / SRB / TRB") + nacelle_ivc.add_discrete_output( + "mb1Type", val="CARB", desc="Type of main bearing: CARB / CRB / SRB / TRB" + ) + nacelle_ivc.add_discrete_output( + "mb2Type", val="SRB", desc="Type of main bearing: CARB / CRB / SRB / TRB" + ) nacelle_ivc.add_discrete_output( "uptower", val=True, desc="If power electronics are located uptower (True) or at tower base (False)" ) @@ -214,7 +221,10 @@ def setup(self): if modeling_options["WISDEM"]["DriveSE"]["direct"]: # Direct only nacelle_ivc.add_output( - "nose_diameter", val=np.zeros(2), units="m", desc="Diameter of nose (also called turret or spindle)" + "nose_diameter", + val=np.zeros(2), + units="m", + desc="Diameter of nose (also called turret or spindle)", ) nacelle_ivc.add_output( "nose_wall_thickness", @@ -231,11 +241,15 @@ def setup(self): else: # Geared only nacelle_ivc.add_output("hss_length", val=0.0, units="m", desc="Length of high speed shaft") - nacelle_ivc.add_output("hss_diameter", val=np.zeros(2), units="m", desc="Diameter of high speed shaft") + nacelle_ivc.add_output( + "hss_diameter", val=np.zeros(2), units="m", desc="Diameter of high speed shaft" + ) nacelle_ivc.add_output( "hss_wall_thickness", val=np.zeros(2), units="m", desc="Wall thickness of high speed shaft" ) - nacelle_ivc.add_output("bedplate_flange_width", val=0.0, units="m", desc="Bedplate I-beam flange width") + nacelle_ivc.add_output( + "bedplate_flange_width", val=0.0, units="m", desc="Bedplate I-beam flange width" + ) nacelle_ivc.add_output( "bedplate_flange_thickness", val=0.0, units="m", desc="Bedplate I-beam flange thickness" ) @@ -248,7 +262,9 @@ def setup(self): desc="3-letter string of Es or Ps to denote epicyclic or parallel gear configuration", ) nacelle_ivc.add_discrete_output( - "planet_numbers", val=[3, 3, 0], desc="Number of planets for epicyclic stages (use 0 for parallel)" + "planet_numbers", + val=[3, 3, 0], + desc="Number of planets for epicyclic stages (use 0 for parallel)", ) # Mulit-body properties @@ -567,10 +583,12 @@ def setup(self): ) opt_var.add_output("af_position", val=np.ones(rotorse_options["n_af_span"])) opt_var.add_output( - "s_opt_spar_cap_ss", val=np.ones(opt_options["design_variables"]["blade"]["structure"]["spar_cap_ss"]["n_opt"]) + "s_opt_spar_cap_ss", + val=np.ones(opt_options["design_variables"]["blade"]["structure"]["spar_cap_ss"]["n_opt"]), ) opt_var.add_output( - "s_opt_spar_cap_ps", val=np.ones(opt_options["design_variables"]["blade"]["structure"]["spar_cap_ps"]["n_opt"]) + "s_opt_spar_cap_ps", + val=np.ones(opt_options["design_variables"]["blade"]["structure"]["spar_cap_ps"]["n_opt"]), ) opt_var.add_output( "spar_cap_ss_opt", @@ -612,8 +630,10 @@ def setup(self): if rotorse_options["inn_af"]: self.add_subsystem( "run_inn_af", - INN_Airfoils(rotorse_options=rotorse_options, - aero_shape_opt_options=opt_options["design_variables"]["blade"]["aero_shape"]), + INN_Airfoils( + rotorse_options=rotorse_options, + aero_shape_opt_options=opt_options["design_variables"]["blade"]["aero_shape"], + ), ) self.connect("outer_shape_bem.s", "run_inn_af.s") self.connect("pa.chord_param", "run_inn_af.chord") @@ -1128,35 +1148,30 @@ def setup(self): val=np.zeros((n_span, n_xy, 2)), desc="3D array of the non-dimensional x and y airfoil coordinates of the airfoils interpolated along span for n_span stations. The leading edge is place at x=0 and y=0.", ) + self.add_input("s_opt_r_thick", val=np.ones(aero_shape_opt_options["t/c"]["n_opt"])) self.add_input( - "s_opt_r_thick", - val=np.ones(aero_shape_opt_options["t/c"]["n_opt"]) - ) - self.add_input( - "r_thick_opt", + "r_thick_opt", val=np.ones(aero_shape_opt_options["t/c"]["n_opt"]), ) + self.add_input("s_opt_L_D", val=np.ones(aero_shape_opt_options["L/D"]["n_opt"])) self.add_input( - "s_opt_L_D", - val=np.ones(aero_shape_opt_options["L/D"]["n_opt"]) - ) - self.add_input( - "L_D_opt", + "L_D_opt", val=np.ones(aero_shape_opt_options["L/D"]["n_opt"]), ) self.add_input( "chord", val=np.zeros(n_span), units="m", desc="1D array of the chord values defined along blade span." ) self.add_input("rated_TSR", val=0.0, desc="Constant tip speed ratio in region II.") - + def compute(self, inputs, outputs): - + try: from INN_interface.INN import INN except: raise Exception("The INN framework for airfoil design is activated, but not installed correctly") import matplotlib - matplotlib.use('TkAgg') + + matplotlib.use("TkAgg") import matplotlib.pyplot as plt from wisdem.ccblade.Polar import Polar @@ -1175,12 +1190,11 @@ def compute(self, inputs, outputs): for i in range(r_thick_index_start, r_thick_index_end): cd = 0.015 - stall_margin = 5. - Re = 9000000. + stall_margin = 5.0 + Re = 9000000.0 inn = INN() try: - cst, alpha = inn.inverse_design(cd, L_D[i], stall_margin, r_thick[i], Re, - N=1, process_samples=True) + cst, alpha = inn.inverse_design(cd, L_D[i], stall_margin, r_thick[i], Re, N=1, process_samples=True) except: raise Exception("The INN for airfoil design failed in the inverse mode") alpha = np.arange(-4, 20, 0.25) @@ -1188,9 +1202,9 @@ def compute(self, inputs, outputs): cd, cl = inn.generate_polars(cst, Re, alpha=alpha) except: raise Exception("The INN for airfoil design failed in the forward mode") - - inn_polar = Polar(Re, alpha, cl[0,:], cd[0,:], np.zeros_like(cl[0,:])) - polar3d = inn_polar.correction3D(inputs["s"][i], inputs["chord"][i]/103., inputs["rated_TSR"]) + + inn_polar = Polar(Re, alpha, cl[0, :], cd[0, :], np.zeros_like(cl[0, :])) + polar3d = inn_polar.correction3D(inputs["s"][i], inputs["chord"][i] / 103.0, inputs["rated_TSR"]) cdmax = 1.5 polar = polar3d.extrapolate(cdmax) # Extrapolate polars for alpha between -180 deg and 180 deg @@ -1198,22 +1212,25 @@ def compute(self, inputs, outputs): cd_interp = np.interp(np.degrees(inputs["aoa"]), polar.alpha, polar.cd) cm_interp = np.interp(np.degrees(inputs["aoa"]), polar.alpha, polar.cm) - f, ax = plt.subplots(3, 1, figsize=(5.3, 8)) - ax[0].plot(inputs["aoa"] * 180. / np.pi, cl_interp, label="INN") - ax[0].plot(inputs["aoa"] * 180. / np.pi, inputs["cl_interp_yaml"][i,:,0,0], label="yaml") + ax[0].plot(inputs["aoa"] * 180.0 / np.pi, cl_interp, label="INN") + ax[0].plot(inputs["aoa"] * 180.0 / np.pi, inputs["cl_interp_yaml"][i, :, 0, 0], label="yaml") ax[0].grid(color=[0.8, 0.8, 0.8], linestyle="--") ax[0].legend() ax[0].set_ylabel("CL (-)", fontweight="bold") ax[0].set_title("Span Location {:2.2%}".format(inputs["s"][i]), fontweight="bold") ax[0].set_xlim(left=-4, right=20) - ax[1].semilogy(inputs["aoa"] * 180. / np.pi, cd_interp, label="INN") - ax[1].semilogy(inputs["aoa"] * 180. / np.pi, inputs["cd_interp_yaml"][i,:,0,0], label="yaml") + ax[1].semilogy(inputs["aoa"] * 180.0 / np.pi, cd_interp, label="INN") + ax[1].semilogy(inputs["aoa"] * 180.0 / np.pi, inputs["cd_interp_yaml"][i, :, 0, 0], label="yaml") ax[1].grid(color=[0.8, 0.8, 0.8], linestyle="--") ax[1].set_ylabel("CD (-)", fontweight="bold") ax[1].set_xlim(left=-4, right=20) - ax[2].plot(inputs["aoa"] * 180. / np.pi, cl_interp/cd_interp, label="INN") - ax[2].plot(inputs["aoa"] * 180. / np.pi, inputs["cl_interp_yaml"][i,:,0,0]/inputs["cd_interp_yaml"][i,:,0,0], label="yaml") + ax[2].plot(inputs["aoa"] * 180.0 / np.pi, cl_interp / cd_interp, label="INN") + ax[2].plot( + inputs["aoa"] * 180.0 / np.pi, + inputs["cl_interp_yaml"][i, :, 0, 0] / inputs["cd_interp_yaml"][i, :, 0, 0], + label="yaml", + ) ax[2].grid(color=[0.8, 0.8, 0.8], linestyle="--") ax[2].set_ylabel("CL/CD (-)", fontweight="bold") ax[2].set_xlabel("Angles of Attack (deg)", fontweight="bold") @@ -1222,6 +1239,7 @@ def compute(self, inputs, outputs): # plt.show() plt.close() + class Blade_Lofted_Shape(om.ExplicitComponent): # Openmdao component to generate the x, y, z coordinates of the points describing the blade outer shape. def initialize(self): @@ -1919,7 +1937,7 @@ class Hub(om.Group): # Openmdao group with the hub data coming from the input yaml file. def initialize(self): self.options.declare("flags") - + def setup(self): ivc = self.add_subsystem("hub_indep_vars", om.IndepVarComp(), promotes=["*"]) @@ -2299,8 +2317,6 @@ def setup(self): n_nodes = mooring_init_options["n_nodes"] n_lines = mooring_init_options["n_lines"] - n_line_types = mooring_init_options["n_line_types"] - # n_anchor_types = mooring_init_options["n_anchor_types"] n_design = 1 if mooring_init_options["symmetric"] else n_lines @@ -2316,8 +2332,6 @@ def setup(self): ivc.add_discrete_output("nodes_joint_name", val=[""] * n_nodes) ivc.add_output("unstretched_length_in", val=np.zeros(n_design), units="m") ivc.add_discrete_output("line_id", val=[""] * n_lines) - # ivc.add_discrete_output("n_lines", val=n_lines) - ivc.add_discrete_output("line_type_names", val=[""] * n_line_types) ## For MoorDyn ivc.add_output("line_diameter_in", val=np.zeros(n_design), units="m") ivc.add_output("line_mass_density_coeff", val=np.zeros(n_lines), units="kg/m**3") ivc.add_output("line_stiffness_coeff", val=np.zeros(n_lines), units="N/m**2") @@ -2327,11 +2341,10 @@ def setup(self): ivc.add_output("line_tangential_added_mass_coeff", val=np.zeros(n_lines), units="kg/m**3") ivc.add_output("line_transverse_drag_coeff", val=np.zeros(n_lines), units="N/m**2") ivc.add_output("line_tangential_drag_coeff", val=np.zeros(n_lines), units="N/m**2") - # ivc.add_discrete_output("anchor_names", val=[""] * n_anchor_types) ivc.add_output("anchor_mass", val=np.zeros(n_lines), units="kg") ivc.add_output("anchor_cost", val=np.zeros(n_lines), units="USD") - ivc.add_output("anchor_max_vertical_load", val=np.zeros(n_lines), units="N") - ivc.add_output("anchor_max_lateral_load", val=np.zeros(n_lines), units="N") + ivc.add_output("anchor_max_vertical_load", val=1e30 * np.ones(n_lines), units="N") + ivc.add_output("anchor_max_lateral_load", val=1e30 * np.ones(n_lines), units="N") self.add_subsystem("moorprop", MooringProperties(mooring_init_options=mooring_init_options), promotes=["*"]) self.add_subsystem("moorjoint", MooringJoints(options=self.options["options"]), promotes=["*"]) @@ -2370,22 +2383,44 @@ def setup(self): def compute(self, inputs, outputs): n_lines = self.options["mooring_init_options"]["n_lines"] + line_mat = self.options["mooring_init_options"]["line_material"] outputs["line_diameter"] = d = inputs["line_diameter_in"] * np.ones(n_lines) outputs["unstretched_length"] = inputs["unstretched_length_in"] * np.ones(n_lines) - d2 = d * d - varlist = [ - "line_mass_density", - "line_stiffness", - "line_breaking_load", - "line_cost_rate", - "line_transverse_added_mass", - "line_tangential_added_mass", - "line_transverse_drag", - "line_tangential_drag", - ] - for var in varlist: - outputs[var] = d2 * inputs[var + "_coeff"] + + line_obj = None + if line_mat[0] == "custom": + varlist = [ + "line_mass_density", + "line_stiffness", + "line_breaking_load", + "line_cost_rate", + "line_transverse_added_mass", + "line_tangential_added_mass", + "line_transverse_drag", + "line_tangential_drag", + ] + for var in varlist: + outputs[var] = d2 * inputs[var + "_coeff"] + + elif line_mat[0] == "chain_stud": + line_obj = mp.getLineProps(1e3 * d[0], type="chain", stud="stud") + else: + line_obj = mp.getLineProps(1e3 * d[0], type=line_mat[0]) + + if not line_obj is None: + outputs["line_mass_density"] = line_obj.mlin + outputs["line_stiffness"] = line_obj.EA + outputs["line_breaking_load"] = line_obj.MBL + outputs["line_cost_rate"] = line_obj.cost + varlist = [ + "line_transverse_added_mass", + "line_tangential_added_mass", + "line_transverse_drag", + "line_tangential_drag", + ] + for var in varlist: + outputs[var] = d2 * inputs[var + "_coeff"] class MooringJoints(om.ExplicitComponent): diff --git a/wisdem/glue_code/gc_WT_InitModel.py b/wisdem/glue_code/gc_WT_InitModel.py index ab3988e0b..4cc03975a 100644 --- a/wisdem/glue_code/gc_WT_InitModel.py +++ b/wisdem/glue_code/gc_WT_InitModel.py @@ -1042,10 +1042,11 @@ def assign_mooring_values(wt_opt, modeling_options, mooring): d2 = mooring["line_types"][ii]["diameter"] ** 2 if jj < n_design: wt_opt["mooring.line_diameter_in"][jj] = mooring["line_types"][ii]["diameter"] - wt_opt["mooring.line_mass_density_coeff"][jj] = mooring["line_types"][ii]["mass_density"] / d2 - wt_opt["mooring.line_stiffness_coeff"][jj] = mooring["line_types"][ii]["stiffness"] / d2 - wt_opt["mooring.line_breaking_load_coeff"][jj] = mooring["line_types"][ii]["breaking_load"] / d2 - wt_opt["mooring.line_cost_rate_coeff"][jj] = mooring["line_types"][ii]["cost"] / d2 + if mooring_init_options["line_material"][jj] == "custom": + wt_opt["mooring.line_mass_density_coeff"][jj] = mooring["line_types"][ii]["mass_density"] / d2 + wt_opt["mooring.line_stiffness_coeff"][jj] = mooring["line_types"][ii]["stiffness"] / d2 + wt_opt["mooring.line_breaking_load_coeff"][jj] = mooring["line_types"][ii]["breaking_load"] / d2 + wt_opt["mooring.line_cost_rate_coeff"][jj] = mooring["line_types"][ii]["cost"] / d2 wt_opt["mooring.line_transverse_added_mass_coeff"][jj] = ( mooring["line_types"][ii]["transverse_added_mass"] / d2 ) @@ -1058,12 +1059,15 @@ def assign_mooring_values(wt_opt, modeling_options, mooring): if node1 == iname or node2 == iname and mooring["nodes"][ii]["node_type"] == "fixed": for kk, kname in enumerate(anchor_names): if kname == mooring["nodes"][ii]["anchor_type"]: - wt_opt["mooring.anchor_mass"][jj] = mooring["anchor_types"][kk]["mass"] - wt_opt["mooring.anchor_cost"][jj] = mooring["anchor_types"][kk]["cost"] - wt_opt["mooring.anchor_max_vertical_load"][jj] = mooring["anchor_types"][kk][ - "max_vertical_load" - ] - wt_opt["mooring.anchor_max_lateral_load"][jj] = mooring["anchor_types"][kk]["max_lateral_load"] + if mooring_init_options["line_anchor"][jj] == "custom": + wt_opt["mooring.anchor_mass"][jj] = mooring["anchor_types"][kk]["mass"] + wt_opt["mooring.anchor_cost"][jj] = mooring["anchor_types"][kk]["cost"] + wt_opt["mooring.anchor_max_vertical_load"][jj] = mooring["anchor_types"][kk][ + "max_vertical_load" + ] + wt_opt["mooring.anchor_max_lateral_load"][jj] = mooring["anchor_types"][kk][ + "max_lateral_load" + ] # Give warnings if we have different types or asymmetrical lines if ( diff --git a/wisdem/glue_code/glue_code.py b/wisdem/glue_code/glue_code.py index a5665be66..25269e265 100644 --- a/wisdem/glue_code/glue_code.py +++ b/wisdem/glue_code/glue_code.py @@ -468,7 +468,10 @@ def setup(self): "fairlead", "fairlead_radius", "anchor_radius", + "anchor_mass", "anchor_cost", + "anchor_max_vertical_load", + "anchor_max_lateral_load", "line_diameter", "line_mass_density_coeff", "line_stiffness_coeff", diff --git a/wisdem/inputs/analysis_schema.yaml b/wisdem/inputs/analysis_schema.yaml index 91a3e8c5a..e38377355 100644 --- a/wisdem/inputs/analysis_schema.yaml +++ b/wisdem/inputs/analysis_schema.yaml @@ -1067,6 +1067,18 @@ properties: description: Keep the mooring line length within the bounds for catenary hang or TLP tension properties: flag: *flag + anchor_vertical: + type: object + default: {} + description: Ensure that the maximum vertical force on the anchor does not exceed limit + properties: + flag: *flag + anchor_lateral: + type: object + default: {} + description: Ensure that the maximum lateral force on the anchor does not exceed limit + properties: + flag: *flag stress: *vonmises global_buckling: *globalb shell_buckling: *shellb diff --git a/wisdem/inputs/geometry_schema.yaml b/wisdem/inputs/geometry_schema.yaml index 025d850da..4c12f3339 100644 --- a/wisdem/inputs/geometry_schema.yaml +++ b/wisdem/inputs/geometry_schema.yaml @@ -1750,16 +1750,24 @@ properties: required: - name - diameter - - mass_density - - stiffness - - breaking_load - - cost - - transverse_added_mass # Optional? Ask MattH <<< I'm not sure if these should be optional or required... will think/discuss more. - - tangential_added_mass # Optional? Ask MattH - - transverse_drag # Optional? Ask MattH - - tangential_drag # Optional? Ask MattH + - type optional: - damping + - transverse_added_mass + - tangential_added_mass + - transverse_drag + - tangential_drag + if: + # If type == custom, then we need all of the details + properties: + type: + const: custom + then: + required: + - mass_density + - stiffness + - breaking_load + - cost properties: name: type: string @@ -1769,6 +1777,10 @@ properties: units: meter description: the volume-equivalent diameter of the line – the diameter of a cylinder having the same displacement per unit length minimum: 0.0 + type: + type: string + enum: [chain, chain_stud, nylon, polyester, polypropylene, wire_fiber, fiber, wire, wire_wire, iwrc, Chain, Chain_Stud, Nylon, Polyester, Polypropylene, Wire, Wire_Fiber, Fiber, Wire, Wire_Wire, IWRC, CHAIN, CHAIN_STUD, NYLON, POLYESTER, POLYPROPYLENE, WIRE, WIRE_FIBER, FIBER, WIRE, WIRE_WIRE, custom, Custom, CUSTOM] + description: Type of material for property lookup mass_density: type: number unit: kilogram/meter @@ -1793,22 +1805,27 @@ properties: type: number unit: Newton * second description: internal damping (BA) + default: 0.0 transverse_added_mass: type: number description: transverse added mass coefficient (with respect to line displacement) minimum: 0.0 + default: 0.0 tangential_added_mass: type: number description: tangential added mass coefficient (with respect to line displacement) minimum: 0.0 + default: 0.0 transverse_drag: type: number description: transverse drag coefficient (with respect to frontal area, d*l) minimum: 0.0 + default: 0.0 tangential_drag: type: number description: tangential drag coefficient (with respect to surface area, π*d*l) minimum: 0.0 + default: 0.0 anchor_types: type: array description: List of anchor properties used in the system @@ -1816,14 +1833,26 @@ properties: type: object required: - name - - mass - - cost - - max_lateral_load - - max_vertical_load + - type + if: + # If type == custom, then we need all of the details + properties: + type: + const: custom + then: + required: + - mass + - cost + - max_lateral_load + - max_vertical_load properties: name: type: string description: Name of anchor to be referenced by anchor_id in Nodes section + type: + type: string + enum: [drag_embedment, suction, plate, micropile, sepla, Drag_Embedment, Suction, Plate, Micropile, Sepla, DRAG_EMBEDMENT, SUCTION, PLATE, MICROPILE, SEPLA, custom, Custom, CUSTOM] + description: Type of anchor for property lookup mass: type: number unit: kilogram @@ -2403,15 +2432,22 @@ properties: fl_feedback: type: object description: Platform velocity feedback using nacelle IMU - required: - - gain - optional: - - filter + default: {} properties: gain: type: number description: Platform velocity feedback gain used by the ROSCO controller (https://github.com/NREL/ROSCO) unit: seconds + twr_freq: + type: number + description: Tower natural frequency / frequency of notch filter in fl_feedback ROSCO control module + unit: rad/s + default: 100. + ptfm_freq: + type: number + description: Platform natural frequency / frequency of LPF in fl_feedback ROSCO control module + unit: rad/s + default: 0.2 filter: $ref: "#/definitions/filter" description: Nacelle IMU filter, used to remove tower natural frequency. Usually a LPF or Notch Filter diff --git a/wisdem/moorpy/Catenary.py b/wisdem/moorpy/Catenary.py index 18084f744..b32c3b676 100644 --- a/wisdem/moorpy/Catenary.py +++ b/wisdem/moorpy/Catenary.py @@ -332,7 +332,7 @@ def catenary(XF, ZF, L, EA, W, CB=0, HF0=0, VF0=0, Tol=0.000001, nNodes=20, MaxI # check for errors ( WOULD SOME NOT ALREADY HAVE BEEN CAUGHT AND RAISED ALREADY?) if info["error"] == True: - breakpoint() + # breakpoint() # >>>> what about errors for which we can first plot the line profile?? <<<< raise CatenaryError("Error in catenary computations: " + info["message"]) diff --git a/wisdem/moorpy/MoorProps.py b/wisdem/moorpy/MoorProps.py index 9e0f39a94..dcf0d3ffb 100644 --- a/wisdem/moorpy/MoorProps.py +++ b/wisdem/moorpy/MoorProps.py @@ -9,7 +9,7 @@ import wisdem.moorpy as mp -def getLineProps(d, type="chain", stud="studless", source="Orcaflex-altered", name=""): +def getLineProps(dmm, type="chain", stud="studless", source="Orcaflex-altered", name=""): """getLineProps version 3.2: Restructuring v3.1 to 'Orcaflex-original' and 'Orcaflex-altered' Motivation: The existing public, and NREL-internal references for mooring line component property @@ -33,7 +33,7 @@ def getLineProps(d, type="chain", stud="studless", source="Orcaflex-altered", na """ if source == "Orcaflex-original": - d = d / 1000 # orcaflex uses meters https://www.orcina.com/webhelp/OrcaFlex/ + d = dmm / 1000 # orcaflex uses meters https://www.orcina.com/webhelp/OrcaFlex/ if type == "chain": c = 1.96e4 # grade 2=1.37e4; grade 3=1.96e4; ORQ=2.11e4; R4=2.74e4 @@ -89,7 +89,7 @@ def getLineProps(d, type="chain", stud="studless", source="Orcaflex-altered", na raise ValueError("getLineProps error: Linetype not valid. Choose from given rope types or chain ") elif source == "Orcaflex-altered": - d = d / 1000 # orcaflex uses meters https://www.orcina.com/webhelp/OrcaFlex/ + d = dmm / 1000 # orcaflex uses meters https://www.orcina.com/webhelp/OrcaFlex/ if type == "chain": c = 2.74e4 # grade 2=1.37e4; grade 3=1.96e4; ORQ=2.11e4; R4=2.74e4 @@ -169,13 +169,13 @@ def getLineProps(d, type="chain", stud="studless", source="Orcaflex-altered", na # Set up a main identifier for the linetype. Useful for things like "chain_bot" or "chain_top" if name == "": - typestring = type + str(d) + typestring = f"{type}{dmm:.0f}" else: typestring = name notes = f"made with getLineProps - source: {source}" - return mp.LineType(typestring, d_vol, massden, EA, MBL=MBL, cost=cost, notes=notes, input_type=type, input_d=d) + return mp.LineType(typestring, d_vol, massden, EA, MBL=MBL, cost=cost, notes=notes, input_type=type, input_d=dmm) def getAnchorProps(fx, fz, type="drag-embedment", display=0): diff --git a/wisdem/moorpy/body.py b/wisdem/moorpy/body.py index 5942e1e01..e6b18940a 100644 --- a/wisdem/moorpy/body.py +++ b/wisdem/moorpy/body.py @@ -76,6 +76,8 @@ def __init__(self, mooringSys, num, type, r6, m=0, v=0, rCG=np.zeros(3), AWP=0, self.sharedLineTheta = [] self.fairR = 0.0 + self.R = np.eye(3) # body orientation rotation matrix + # print("Created Body "+str(self.number)) def attachPoint(self, pointID, rAttach): @@ -125,9 +127,11 @@ def setPosition(self, r6): f"Body setPosition method requires an argument of size 6, but size {len(r6):d} was provided" ) + self.R = rotationMatrix(self.r6[3], self.r6[4], self.r6[5]) # update body rotation matrix + # update the position of any attached Points for PointID, rPointRel in zip(self.attachedP, self.rPointRel): - rPoint = transformPosition(rPointRel, r6) + rPoint = np.matmul(self.R, rPointRel) + self.r6[:3] # rPoint = transformPosition(rPointRel, r6) self.sys.pointList[PointID - 1].setPosition(rPoint) if self.sys.display > 3: @@ -365,3 +369,6 @@ def redraw(self): linebit[2][0].set_3d_properties([self.r6[2], rz[2]]) return linebit + + +# diff --git a/wisdem/moorpy/helpers.py b/wisdem/moorpy/helpers.py index e18af954b..76ea8c214 100644 --- a/wisdem/moorpy/helpers.py +++ b/wisdem/moorpy/helpers.py @@ -35,6 +35,14 @@ def __init__(self, message): self.message = message +# Generic MoorPy error +class MoorPyError(Error): + """Derived error class for MoorPy. Contains an error message""" + + def __init__(self, message): + self.message = str(message) + + def printMat(mat): """Prints a matrix to a format that is specified @@ -393,9 +401,6 @@ def step_func(X, args, Y, oths, Ytarget, err, tols, iter, maxIter): else: dX = step_func(X, args, Y, oths, Ytarget, err, tols, iter, maxIter) - # if display>2: - # breakpoint() - # Make sure we're not diverging by keeping things from reversing too much. # Track the previous step (dX_last) and if the current step reverses too much, stop it part way. # Stop it at a plane part way between the current X value and the previous X value (using golden ratio, why not). @@ -427,7 +432,6 @@ def step_func(X, args, Y, oths, Ytarget, err, tols, iter, maxIter): ) # set the maximum permissible dx in each direction based an an acceleration limit if dX_max == 0.0: # avoid a divide-by-zero case (if dX[i] was zero to start with) - breakpoint() dX[i] = 0.0 else: a_i = dX[i] / dX_max # calculate ratio of desired dx to max dx @@ -445,8 +449,6 @@ def step_func(X, args, Y, oths, Ytarget, err, tols, iter, maxIter): print(f" now dX will be {dX}") dXlist[iter, :] = dX - if iter == 196: - breakpoint() # enforce bounds for i in range(N): @@ -472,9 +474,6 @@ def step_func(X, args, Y, oths, Ytarget, err, tols, iter, maxIter): ) print("Solution X is " + str(X)) - # if abs(err) > 10: - # breakpoint() - if display > 0: print( "Total run time: {:8.2f} seconds = {:8.2f} minutes".format( diff --git a/wisdem/moorpy/system.py b/wisdem/moorpy/system.py index 5bd4054fe..6e2ea2bfc 100644 --- a/wisdem/moorpy/system.py +++ b/wisdem/moorpy/system.py @@ -7,7 +7,16 @@ from wisdem.moorpy.point import Point # import wisdem.moorpy.MoorSolve as msolve -from wisdem.moorpy.helpers import SolveError, getH, dsolve2, printVec, rotatePosition, rotationMatrix, set_axes_equal +from wisdem.moorpy.helpers import ( + SolveError, + MoorPyError, + getH, + dsolve2, + printVec, + rotatePosition, + rotationMatrix, + set_axes_equal, +) from wisdem.moorpy.lineType import LineType @@ -345,7 +354,7 @@ def load(self, filename): "Generic Fairlead/Vessel-type points aren't supported when bodies are defined." ) if len(self.bodyList) == 0: - print("Adding a body to attach fairlead points to.") + # print("Adding a body to attach fairlead points to.") self.bodyList.append(Body(self, 1, 0, np.zeros(6))) # , m=m, v=v, rCG=rCG) ) rRel = np.array(entries[2:5], dtype=float) @@ -493,7 +502,7 @@ def parseYAML(self, data): if len(self.bodyList) > 1: raise ValueError("Generic Fairlead/Vessel-type points aren't supported when bodies are defined.") if len(self.bodyList) == 0: - print("Adding a body to attach fairlead points to.") + # print("Adding a body to attach fairlead points to.") self.bodyList.append(Body(self, 1, 0, np.zeros(6))) # , m=m, v=v, rCG=rCG) ) rRel = np.array(d["location"], dtype=float) @@ -540,7 +549,7 @@ def parseYAML(self, data): elif "water_density" in data: self.rho = data["water_density"] - def unload(self, fileName, **kwargs): + def unload(self, fileName, MDversion=2, **kwargs): """Unloads a MoorPy system into a MoorDyn-style input file Parameters @@ -595,8 +604,8 @@ def unload(self, fileName, **kwargs): # Point Properties (for each point in pointList) #! Add Comments - CdA = 0 - Ca = 0 + CdA = 0.0 + Ca = 0.0 # Line Properties flag = "p" # "-" @@ -614,36 +623,37 @@ def unload(self, fileName, **kwargs): Outputs = ["FairTen1", "FairTen2", "FairTen3"] #! Standard Option (Fairing Tenstion for num of lines) - print("attempting to write " + fileName + " using MoorDyn v?.??") + print("attempting to write " + fileName + " for MoorDyn v" + str(MDversion)) # Array to add strings to for each line of moordyn input file L = [] # Input File Header - L.append("------------------- MoorDyn v?.??. Input File ------------------------------------") + L.append(f" MoorDyn v{MDversion} Input File ") if "description" in locals(): L.append("MoorDyn input for " + description) else: - L.append("MoorDyn input") - #! - L.append("{:5} Echo - echo the input file data (flag)".format(str(Echo).upper())) + L.append("Generated by MoorPy") + + # L.append("{:5} Echo - echo the input file data (flag)" + # .format(str(Echo).upper())) # Line Dictionary Header - L.append("---------------------- LINE DICTIONARY -----------------------------------------------------") - L.append("LineType Diam MassDenInAir EA cIntDamp EI Can Cat Cdn Cdt") - L.append("(-) (m) (kg/m) (N) (Pa-s) (N-m^2) (-) (-) (-) (-)") + L.append("---------------------- LINE TYPES -----------------------------------------------------") + L.append("LineType Diam MassDen EA cIntDamp EI Can Cat Cdn Cdt") + L.append(" (-) (m) (kg/m) (N) (Pa-s) (N-m^2) (-) (-) (-) (-)") # Line Dicationary Table for key in self.lineTypes: # for key,value in self.lineTypes.items(): (Another way to iterate through dictionary) L.append( - "{:<10}{:<10.4f}{:<10.4f} {:<10.5e} ".format( + "{:<15} {:7.4f} {:8.2f} {:7.3e} ".format( key, self.lineTypes[key].d, self.lineTypes[key].mlin, self.lineTypes[key].EA ) - + "{:<11.4f}{:<7.3f}{:<8.3f}{:<7.3f}{:<8.3f}{:<8.3f}".format(cIntDamp, EI, Can, Cat, Cdn, Cdt) + + "{:7.3e} {:7.3e} {:<7.3f} {:<7.3f} {:<7.3f} {:<7.3f}".format(cIntDamp, EI, Can, Cat, Cdn, Cdt) ) # Rod Dictionary Header - L.append("--------------------- ROD DICTIONARY -----------------------------------------------------") + L.append("--------------------- ROD TYPES -----------------------------------------------------") L.append("RodType Diam MassDenInAir Can Cat Cdn Cdt ") L.append("(-) (m) (kg/m) (-) (-) (-) (-) ") @@ -664,7 +674,7 @@ def unload(self, fileName, **kwargs): # Body List Table for body in self.bodyList: L.append( - " {:<7d} {:<5.2f} {:<5.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<8.2f} {:<7.2f} {:<11.2f}".format( + " {:<4d} {:<5.2f} {:<5.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<6.2f} {:<8.2f} {:<7.2f} {:<11.2f}".format( body.number, body.r6[0], body.r6[1], @@ -678,7 +688,7 @@ def unload(self, fileName, **kwargs): body.m, body.v, ) - + "{:<9d}{:<9d}{:<7d}{:<2d}{:<2d}{:<8d}{:<1d}".format( + + "{:<9d} {:<9d} {:<7d} {:<2d} {:<2d} {:<8d} {:<1d}".format( IX, IY, IZ, CdA_xyz[0], CdA_xyz[1], CdA_xyz[2], Ca_xyz[0], Ca_xyz[1], Ca_xyz[2] ) ) @@ -693,9 +703,9 @@ def unload(self, fileName, **kwargs): """ # Point Properties Header - L.append("---------------------- POINT LIST -----------------------------------------------------") - L.append("Node Type X Y Z M V FX FY FZ CdA Ca ") - L.append("(-) (-) (m) (m) (m) (kg) (m^3) (kN) (kN) (kN) (m2) ()") + L.append("---------------------- POINT LIST ---------------------------------------------------------") + L.append("Node Type X Y Z M V FX FY FZ CdA Ca ") + L.append("(-) (-) (m) (m) (m) (kg) (m^3) (kN) (kN) (kN) (m2) ()") # Point Properties Table for point in self.pointList: @@ -721,7 +731,7 @@ def unload(self, fileName, **kwargs): point_type = "Vessel" L.append( - "{:<9d} {:8} {:<20.15f} {:<20.15f} {:<11.4f} {:<7.2f} {:<7.2f} {:<5.2f} {:<5.2f} {:<7.2f}".format( + "{:<4d} {:12} {:8.2f} {:8.2f} {:8.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f} {:6.2f}".format( point.number, point_type, point_pos[0], @@ -732,14 +742,15 @@ def unload(self, fileName, **kwargs): point.fExt[0], point.fExt[1], point.fExt[2], + CdA, + Ca, ) - + " {:<6d} {:<1d}".format(CdA, Ca) ) # Line Properties Header L.append("---------------------- LINE LIST -----------------------------------------------------") - L.append("Line LineType UnstrLen NumSegs NodeAnch NodeFair Flags/Outputs") - L.append("(-) (-) (m) (-) (-) (-) (-)") + L.append("Line LineType UnstrLen NumSegs AttachA AttachB Outputs") + L.append("(-) (-) (m) (-) (-) (-) (-)") # Line Properties Table # (Create a ix2 array of connection points from a list of m points) @@ -757,7 +768,7 @@ def unload(self, fileName, **kwargs): # Populate text for i in range(len(self.lineList)): L.append( - "{:<9d}{:<11}{:<11.2f}{:<11d}{:<10d}{:<10d}{}".format( + "{:<4d} {:<15} {:8.3f} {:5d} {:7d} {:8d} {}".format( self.lineList[i].number, self.lineList[i].type, self.lineList[i].L, @@ -1234,9 +1245,6 @@ def solveEquilibrium(self, DOFtype="free", plots=0, rmsTol=10, maxIter=200): elif dX[i] < -db[i]: dX[i] = -db[i] - # if iter == 6: - # breakpoint() - X0 = X0 + 1.0 * dX # print(X0) # print(self.mooringEq(X0)) @@ -1292,7 +1300,7 @@ def solveEquilibrium3( if rmsTol != 0.0: tols = np.zeros(len(X0)) + rmsTol print("WHAT IS PASSING rmsTol in to solveEquilibrium3?") - breakpoint() + elif np.isscalar(tol): if tol < 0: tols = -tol * db # tolerances set relative to max step size @@ -1347,9 +1355,6 @@ def step_func_equil(X, args, Y, oths, Ytarget, err, tol_, iter, maxIter): # dX = np.matmul(np.linalg.inv(K), Y) # calculate position adjustment according to Newton's method dX = np.linalg.solve(K, Y) # calculate position adjustment according to Newton's method - if iter == 196: - breakpoint() - # but limit adjustment to keep things under control for i in range(n): if dX[i] > db[i]: @@ -1366,7 +1371,12 @@ def step_func_equil(X, args, Y, oths, Ytarget, err, tol_, iter, maxIter): # Call dsolve function # X, Y, info = msolve.dsolve(eval_func_equil, X0, step_func=step_func_equil, tol=tol, maxIter=maxIter) - X, Y, info = dsolve2(eval_func_equil, X0, step_func=step_func_equil, tol=tols, maxIter=maxIter, display=display) + try: + X, Y, info = dsolve2( + eval_func_equil, X0, step_func=step_func_equil, tol=tols, maxIter=maxIter, display=display + ) + except Exception as e: + raise MoorPyError(e) # Don't need to call Ytarget in dsolve because it's already set to be zeros if self.display > 1: @@ -1394,22 +1404,22 @@ def step_func_equil(X, args, Y, oths, Ytarget, err, tol_, iter, maxIter): print(f"current system stiffness: {K}") print(f"\n Current force {F}") - # plot the convergence failure - if n < 8: - fig, ax = plt.subplots(2 * n, 1, sharex=True) - for i in range(n): - ax[i].plot(info["Xs"][: info["iter"] + 1, i]) - ax[n + i].plot(info["Es"][: info["iter"] + 1, i]) - ax[-1].set_xlabel("iteration") - else: - fig, ax = plt.subplots(n, 2, sharex=True) - for i in range(n): - ax[i, 0].plot(info["Xs"][: info["iter"] + 1, i]) - ax[i, 1].plot(info["Es"][: info["iter"] + 1, i]) - ax[-1, 0].set_xlabel("iteration, X") - ax[-1, 1].set_xlabel("iteration, Error") - plt.show() - breakpoint() + # plot the convergence failure + if n < 8: + fig, ax = plt.subplots(2 * n, 1, sharex=True) + for i in range(n): + ax[i].plot(info["Xs"][: info["iter"] + 1, i]) + ax[n + i].plot(info["Es"][: info["iter"] + 1, i]) + ax[-1].set_xlabel("iteration") + else: + fig, ax = plt.subplots(n, 2, sharex=True) + for i in range(n): + ax[i, 0].plot(info["Xs"][: info["iter"] + 1, i]) + ax[i, 1].plot(info["Es"][: info["iter"] + 1, i]) + ax[-1, 0].set_xlabel("iteration, X") + ax[-1, 1].set_xlabel("iteration, Error") + plt.show() + raise SolveError( f"solveEquilibrium3 failed to find equilibrium after {iter} iterations, with residual forces of {F}" ) @@ -1939,9 +1949,9 @@ def plot(self, bounds="default", ax=None, color=None, hidebox=False, rbound=0, t if ax == None: fig = plt.figure() # fig = plt.figure(figsize=(20/2.54,12/2.54), dpi=300) - ax = fig.gca(projection="3d") + ax = plt.axes(projection="3d") else: - fig = plt.gcf() # will this work like this? <<< + fig = ax.get_figure() # set bounds if rbound == 0: diff --git a/wisdem/pyframe3dd/pyframe3dd.py b/wisdem/pyframe3dd/pyframe3dd.py index e1dd3cfd6..bacd8f6e3 100644 --- a/wisdem/pyframe3dd/pyframe3dd.py +++ b/wisdem/pyframe3dd/pyframe3dd.py @@ -8,14 +8,15 @@ """ from __future__ import print_function -import numpy as np + +import os import math -from ctypes import POINTER, c_int, c_double, Structure, pointer +from sys import platform +from ctypes import POINTER, Structure, c_int, pointer, c_double from collections import namedtuple -import os from distutils.sysconfig import get_config_var -from sys import platform +import numpy as np libext = get_config_var("EXT_SUFFIX") if libext is None or libext == "": @@ -388,16 +389,14 @@ def __init__(self, nodes, reactions, elements, options): rigid = 1 else: self.rnode = np.array(reactions.node).astype(np.int32) - self.rKx = np.array(reactions.Kx).astype( - np.float64 - ) # convert rather than copy to allow old syntax of integers + # convert rather than copy to allow old syntax of integers + self.rKx = np.array(reactions.Kx).astype(np.float64) self.rKy = np.array(reactions.Ky).astype(np.float64) self.rKz = np.array(reactions.Kz).astype(np.float64) self.rKtx = np.array(reactions.Ktx).astype(np.float64) self.rKty = np.array(reactions.Kty).astype(np.float64) self.rKtz = np.array(reactions.Ktz).astype(np.float64) rigid = reactions.rigid - # elements self.eelement = np.array(elements.element).astype(np.int32) self.eN1 = np.array(elements.N1).astype(np.int32) diff --git a/wisdem/test/test_floatingse/test_floating.py b/wisdem/test/test_floatingse/test_floating.py index dfc7217e1..33f726ca8 100644 --- a/wisdem/test/test_floatingse/test_floating.py +++ b/wisdem/test/test_floatingse/test_floating.py @@ -42,6 +42,8 @@ def testMassPropertiesSpar(self): opt["mooring"] = {} opt["mooring"]["n_attach"] = 3 opt["mooring"]["n_anchors"] = 3 + opt["mooring"]["line_anchor"] = ["custom"] * 3 + opt["mooring"]["line_material"] = ["custom"] * 3 opt["materials"] = {} opt["materials"]["n_mat"] = 2 @@ -96,6 +98,7 @@ def testMassPropertiesSpar(self): prob["fairlead"] = 70.0 prob["anchor_radius"] = 853.87 # Distance from centerline to sea floor landing [m] prob["anchor_cost"] = 1e5 + prob["anchor_mass"] = 0.0 # Mooring constraints prob["max_surge_fraction"] = 0.1 # Max surge/sway offset [m] diff --git a/wisdem/test/test_floatingse/test_frame.py b/wisdem/test/test_floatingse/test_frame.py index d79fda50c..bf31b3156 100644 --- a/wisdem/test/test_floatingse/test_frame.py +++ b/wisdem/test/test_floatingse/test_frame.py @@ -137,6 +137,7 @@ def setUp(self): self.inputs["mooring_neutral_load"][:, 1] = [0.0, 50, -50] self.inputs["mooring_neutral_load"][:, 2] = -1e3 self.inputs["mooring_fairlead_joints"] = np.array([[0.0, 0.0, 0.0], [0.5, 1.0, 0.0], [1.0, 0.0, 0.0]]) + self.inputs["mooring_stiffness"] = 5 * np.eye(6) self.inputs["transition_node"] = self.inputs["tower_nodes"][0, :] self.inputs["tower_top_node"] = self.inputs["tower_nodes"][2, :] self.inputs["rna_mass"] = 1e4 diff --git a/wisdem/test/test_floatingse/test_mooring.py b/wisdem/test/test_floatingse/test_mooring.py index 46d7d51c0..82be653ba 100644 --- a/wisdem/test/test_floatingse/test_mooring.py +++ b/wisdem/test/test_floatingse/test_mooring.py @@ -9,13 +9,15 @@ class TestMooring(unittest.TestCase): def setUp(self): self.inputs = {} self.outputs = {} - self.discrete_inputs = {} myones = np.ones(1) self.inputs["fairlead"] = 10.0 * myones self.inputs["fairlead_radius"] = 11.0 * myones self.inputs["anchor_radius"] = 250.0 * myones - self.inputs["anchor_cost"] = 10.0 + self.inputs["anchor_mass"] = 100.0 * myones + self.inputs["anchor_cost"] = 10.0 * myones + self.inputs["anchor_max_vertical_load"] = 1e10 + self.inputs["anchor_max_lateral_load"] = 1e10 self.inputs["rho_water"] = 1e3 self.inputs["water_depth"] = 200.0 # 100.0 @@ -33,6 +35,8 @@ def setUp(self): opt = {} opt["n_attach"] = 3 opt["n_anchors"] = 6 + opt["line_anchor"] = ["custom"] # ["drag_embedment"] + opt["line_material"] = ["custom"] # ["chain"] self.mymap = mm.Mooring(options=opt, gamma=1.35) @@ -66,9 +70,23 @@ def testRunMap(self): self.outputs["survival_heel_restoring_force"][4], self.outputs["operational_heel_restoring_force"][4] ) - self.assertAlmostEqual(self.outputs["mooring_mass"], 6 * 270 * 0.1) + self.assertAlmostEqual(self.outputs["mooring_mass"], 6 * 270 * 0.1 + 6 * 100) self.assertAlmostEqual(self.outputs["mooring_cost"], 6 * 270 * 0.4 + 6 * 10) + def testRunMap_MoorProps(self): + + opt = {} + opt["n_attach"] = 3 + opt["n_anchors"] = 6 + opt["line_anchor"] = ["drag_embedment"] + opt["line_material"] = ["chain"] + + mymap = mm.Mooring(options=opt, gamma=1.35) + + mymap.compute(self.inputs, self.outputs) + self.assertAlmostEqual(self.outputs["mooring_mass"], 6 * 270 * 199.0 + 6 * 0) + self.assertAlmostEqual(self.outputs["mooring_cost"], 6 * 270 * 514.415) + def suite(): suite = unittest.TestSuite()