diff --git a/WISDEM/wisdem/glue_code/gc_PoseOptimization.py b/WISDEM/wisdem/glue_code/gc_PoseOptimization.py index 9b3cc11c4..9efb7b1e4 100644 --- a/WISDEM/wisdem/glue_code/gc_PoseOptimization.py +++ b/WISDEM/wisdem/glue_code/gc_PoseOptimization.py @@ -134,12 +134,12 @@ def get_number_design_variables(self): n_DV += 1 if "axial_joints" in kgrp: n_DV += len(kgrp["axial_joints"]) - - n_design = 1 if self.modeling["mooring"]["symmetric"] else self.modeling["mooring"]["n_lines"] - if mooring_opt["line_length"]["flag"]: - n_DV += n_design - if mooring_opt["line_diameter"]["flag"]: - n_DV += n_design + if self.modeling["flags"]["mooring"]: + n_design = 1 if self.modeling["mooring"]["symmetric"] else self.modeling["mooring"]["n_lines"] + if mooring_opt["line_length"]["flag"]: + n_DV += n_design + if mooring_opt["line_diameter"]["flag"]: + n_DV += n_design # Wrap-up at end with multiplier for finite differencing if self.opt["driver"]["optimization"]["form"] == "central": diff --git a/WISDEM/wisdem/glue_code/gc_WT_InitModel.py b/WISDEM/wisdem/glue_code/gc_WT_InitModel.py index e2c93e873..6e2be56de 100644 --- a/WISDEM/wisdem/glue_code/gc_WT_InitModel.py +++ b/WISDEM/wisdem/glue_code/gc_WT_InitModel.py @@ -416,10 +416,11 @@ def assign_internal_structure_2d_fem_values(wt_opt, modeling_options, internal_s else: definition_layer[i] = 6 flag = False - if layer_name[k] == internal_structure_2d_fem["layers"][i]["end_nd_arc"]["fixed"]: - index_layer_end[i] = k - flag = True - break + for k in range(n_layers): + if layer_name[k] == internal_structure_2d_fem["layers"][i]["end_nd_arc"]["fixed"]: + index_layer_end[i] = k + flag = True + break if flag == False: raise ValueError("Error with layer " + internal_structure_2d_fem["layers"][i]["name"]) else: diff --git a/WISDEM/wisdem/inputs/analysis_schema.yaml b/WISDEM/wisdem/inputs/analysis_schema.yaml index 54954f36a..d601a569e 100644 --- a/WISDEM/wisdem/inputs/analysis_schema.yaml +++ b/WISDEM/wisdem/inputs/analysis_schema.yaml @@ -85,7 +85,7 @@ properties: description: First index of the array of design variables/constraints that is optimized/constrained index_end: &index_end type: integer - default: 0 + default: 8 minimum: 0 unit: none description: Last index of the array of design variables/constraints that is optimized/constrained diff --git a/WISDEM/wisdem/rotorse/rotor_structure.py b/WISDEM/wisdem/rotorse/rotor_structure.py index b3adb1916..ee9a50eef 100644 --- a/WISDEM/wisdem/rotorse/rotor_structure.py +++ b/WISDEM/wisdem/rotorse/rotor_structure.py @@ -94,9 +94,15 @@ def setup(self): ) # Outputs - self.add_output("Px_af", val=np.zeros(n_span), units="N/m", desc="total distributed loads in airfoil x-direction") - self.add_output("Py_af", val=np.zeros(n_span), units="N/m", desc="total distributed loads in airfoil y-direction") - self.add_output("Pz_af", val=np.zeros(n_span), units="N/m", desc="total distributed loads in airfoil z-direction") + self.add_output( + "Px_af", val=np.zeros(n_span), units="N/m", desc="total distributed loads in airfoil x-direction" + ) + self.add_output( + "Py_af", val=np.zeros(n_span), units="N/m", desc="total distributed loads in airfoil y-direction" + ) + self.add_output( + "Pz_af", val=np.zeros(n_span), units="N/m", desc="total distributed loads in airfoil z-direction" + ) def compute(self, inputs, outputs): @@ -163,13 +169,22 @@ def setup(self): # all inputs/outputs in airfoil coordinate system self.add_input( - "Px_af", val=np.zeros(n_span), units="N/m", desc="distributed load (force per unit length) in airfoil x-direction" + "Px_af", + val=np.zeros(n_span), + units="N/m", + desc="distributed load (force per unit length) in airfoil x-direction", ) self.add_input( - "Py_af", val=np.zeros(n_span), units="N/m", desc="distributed load (force per unit length) in airfoil y-direction" + "Py_af", + val=np.zeros(n_span), + units="N/m", + desc="distributed load (force per unit length) in airfoil y-direction", ) self.add_input( - "Pz_af", val=np.zeros(n_span), units="N/m", desc="distributed load (force per unit length) in airfoil z-direction" + "Pz_af", + val=np.zeros(n_span), + units="N/m", + desc="distributed load (force per unit length) in airfoil z-direction", ) self.add_input( @@ -365,7 +380,7 @@ def compute(self, inputs, outputs): EIxy_cs -= x_ec_cs * y_ec_cs * EA # get rotation angle - alpha = 0.5 * np.arctan(2 * EIxy_cs / (EIyy_cs - EIxx_cs)) + alpha = 0.5 * np.arctan2(2 * EIxy_cs, (EIyy_cs - EIxx_cs)) # get moments and positions in principal axes EI11 = EIxx_cs - EIxy_cs * np.tan(alpha) @@ -387,10 +402,10 @@ def rotate(x, y): rad = np.zeros(n) # 'radius' of rigidity at node- set to zero inode = 1 + np.arange(n) # Node numbers (1-based indexing) if self.options["pbeam"]: - nodes = pyframe3dd.NodeData(inode, np.zeros(n), np.zeros(n), r, rad) + nodes = pyframe3dd.NodeData(inode, np.zeros(n), np.zeros(n), -r, rad) L = np.diff(r) else: - nodes = pyframe3dd.NodeData(inode, x_az, y_az, z_az, rad) + nodes = pyframe3dd.NodeData(inode, y_az, x_az, -z_az, rad) L = np.sqrt(np.diff(x_az) ** 2 + np.diff(y_az) ** 2 + np.diff(z_az) ** 2) # ----------------------------------- @@ -471,7 +486,7 @@ def rotate(x, y): Py_af = P.y Pz_af = P.z - Px, Py, Pz = Pz_af, Py_af, -Px_af # switch to local c.s. + Px, Py, Pz = Pz_af, Py_af, Px_af # switch to local c.s. xx1 = xy1 = xz1 = np.zeros(n - 1) xx2 = xy2 = xz2 = L - 1e-6 # subtract small number b.c. of precision wx1 = Px[:-1] @@ -492,11 +507,6 @@ def rotate(x, y): # For now, just 1 load case and blade iCase = 0 - # Displacements in global (blade) c.s. - dx = displacements.dx[iCase, :] - dy = displacements.dy[iCase, :] - dz = displacements.dz[iCase, :] - # Mode shapes and frequencies n_freq2 = int(self.n_freq / 2) freq_x, freq_y, _, mshapes_x, mshapes_y, _ = util.get_xyz_mode_shapes( @@ -509,12 +519,12 @@ def rotate(x, y): # shear and bending, one per element (convert from local to global c.s.) Fz = np.r_[-forces.Nx[iCase, 0], forces.Nx[iCase, 1::2]] - Vy = np.r_[-forces.Vy[iCase, 0], forces.Vy[iCase, 1::2]] - Vx = np.r_[forces.Vz[iCase, 0], -forces.Vz[iCase, 1::2]] + # Vy = np.r_[-forces.Vy[iCase, 0], forces.Vy[iCase, 1::2]] + # Vx = np.r_[-forces.Vz[iCase, 0], forces.Vz[iCase, 1::2]] - Tz = np.r_[-forces.Txx[iCase, 0], forces.Txx[iCase, 1::2]] + # Tz = np.r_[-forces.Txx[iCase, 0], forces.Txx[iCase, 1::2]] My = np.r_[-forces.Myy[iCase, 0], forces.Myy[iCase, 1::2]] - Mx = np.r_[forces.Mzz[iCase, 0], -forces.Mzz[iCase, 1::2]] + Mx = np.r_[-forces.Mzz[iCase, 0], forces.Mzz[iCase, 1::2]] def strain(xu, yu, xl, yl): # use profile c.s. to use Hansen's notation @@ -529,12 +539,10 @@ def strain(xu, yu, xl, yl): # compute strain x, y = rotate(xuu, yuu) - strainU = -( - M1 / EI11 * y - M2 / EI22 * x + Fz / EA - ) # negative sign because Hansen c3 is opposite of Precomp z + strainU = M1 / EI11 * y - M2 / EI22 * x - Fz / EA x, y = rotate(xll, yll) - strainL = -(M1 / EI11 * y - M2 / EI22 * x + Fz / EA) + strainL = M1 / EI11 * y - M2 / EI22 * x - Fz / EA return strainU, strainL @@ -553,9 +561,10 @@ def strain(xu, yu, xl, yl): outputs["edge_mode_freqs"] = freq_y outputs["flap_mode_freqs"] = freq_x outputs["freq_distance"] = freq_y[0] / freq_x[0] - outputs["dx"] = dx - outputs["dy"] = dy - outputs["dz"] = dz + # Displacements in global (blade) c.s. + outputs["dx"] = displacements.dx[iCase, :] + outputs["dy"] = displacements.dy[iCase, :] + outputs["dz"] = displacements.dz[iCase, :] outputs["strainU_spar"] = strainU_spar outputs["strainL_spar"] = strainL_spar outputs["strainU_te"] = strainU_te @@ -756,19 +765,19 @@ class BladeRootSizing(ExplicitComponent): Recommended diameter of the blade root fastener circle ratio : float Ratio of recommended diameter over actual diameter. It can be constrained to be smaller than 1 - + """ def initialize(self): self.options.declare("rotorse_options") def setup(self): - + rotorse_options = self.options["rotorse_options"] self.n_span = n_span = rotorse_options["n_span"] self.n_layers = n_layers = rotorse_options["n_layers"] - self.add_input("rootD", val=0., units="m", desc="Blade root outer diameter / Chord at blade span station 0") + self.add_input("rootD", val=0.0, units="m", desc="Blade root outer diameter / Chord at blade span station 0") self.add_input( "layer_thickness", val=np.zeros((n_layers, n_span)), @@ -787,33 +796,37 @@ def setup(self): ) self.add_input("root_M", val=np.zeros(3), units="N*m", desc="Blade root moment in blade c.s.") self.add_input("s_f", val=rotorse_options["root_fastener_s_f"], desc="Safety factor") - self.add_input("d_f", val=0., units="m", desc="Diameter of the fastener") - self.add_input("sigma_max", val=0., units="Pa", desc="Max stress on bolt") + self.add_input("d_f", val=0.0, units="m", desc="Diameter of the fastener") + self.add_input("sigma_max", val=0.0, units="Pa", desc="Max stress on bolt") - self.add_output("d_r", val=0., units="m", desc="Root fastener circle diameter") - self.add_output("ratio", val=0., desc="Ratio of recommended diameter over actual diameter. It can be constrained to be smaller than 1") + self.add_output("d_r", val=0.0, units="m", desc="Root fastener circle diameter") + self.add_output( + "ratio", + val=0.0, + desc="Ratio of recommended diameter over actual diameter. It can be constrained to be smaller than 1", + ) def compute(self, inputs, outputs): - - Mxy = np.sqrt(inputs["root_M"][0]**2. + inputs["root_M"][1]**2.) - d_r = np.sqrt((48. * Mxy * inputs["s_f"])/(np.pi**2. * inputs["sigma_max"] * inputs["d_f"] )) - + Mxy = np.sqrt(inputs["root_M"][0] ** 2.0 + inputs["root_M"][1] ** 2.0) + + d_r = np.sqrt((48.0 * Mxy * inputs["s_f"]) / (np.pi ** 2.0 * inputs["sigma_max"] * inputs["d_f"])) + sectors = np.array([]) for i in range(self.n_layers): - sectors = np.unique(np.hstack([sectors, inputs["layer_start_nd"][i,0], inputs["layer_end_nd"][i,0]])) + sectors = np.unique(np.hstack([sectors, inputs["layer_start_nd"][i, 0], inputs["layer_end_nd"][i, 0]])) thick = np.zeros(len(sectors)) for j in range(len(sectors)): for i in range(self.n_layers): - if inputs["layer_start_nd"][i,0] <= sectors[j] and inputs["layer_end_nd"][i,0] >= sectors[j]: - thick[j] += inputs["layer_thickness"][i,0] + if inputs["layer_start_nd"][i, 0] <= sectors[j] and inputs["layer_end_nd"][i, 0] >= sectors[j]: + thick[j] += inputs["layer_thickness"][i, 0] # check = np.all(thick == thick[0]) # if not check: # raise Exception('All Values in Array are not same') d_r_actual = inputs["rootD"] - 0.5 * thick[0] - + ratio = d_r / d_r_actual outputs["d_r"] = d_r @@ -1176,4 +1189,4 @@ def setup(self): self.connect("frame.edge_mode_freqs", "constr.edge_mode_freqs") # Blade root moment to blade root sizing - self.connect("frame.root_M", "brs.root_M") \ No newline at end of file + self.connect("frame.root_M", "brs.root_M") diff --git a/WISDEM/wisdem/test/test_gluecode/test_gc_modified_yaml.py b/WISDEM/wisdem/test/test_gluecode/test_gc_modified_yaml.py index 3039cf185..9b06a16dd 100644 --- a/WISDEM/wisdem/test/test_gluecode/test_gc_modified_yaml.py +++ b/WISDEM/wisdem/test/test_gluecode/test_gc_modified_yaml.py @@ -32,7 +32,7 @@ def test15MW(self): self.assertAlmostEqual(wt_opt["re.precomp.blade_mass"][0], 69911.6542917299, 1) self.assertAlmostEqual(wt_opt["rp.AEP"][0] * 1.0e-6, 77.84277082383369, 1) self.assertAlmostEqual(wt_opt["financese.lcoe"][0] * 1.0e3, 64.1372990450, 1) - self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 22.8104130129, 1) + self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 25.4598632918, 1) self.assertAlmostEqual(wt_opt["towerse.z_param"][-1], 144.386, 3) diff --git a/WISDEM/wisdem/test/test_gluecode/test_gluecode.py b/WISDEM/wisdem/test/test_gluecode/test_gluecode.py index 59b404871..44aa18573 100644 --- a/WISDEM/wisdem/test/test_gluecode/test_gluecode.py +++ b/WISDEM/wisdem/test/test_gluecode/test_gluecode.py @@ -27,8 +27,8 @@ def test5MW(self): self.assertAlmostEqual(wt_opt["re.precomp.blade_mass"][0], 16403.682326940743, 2) self.assertAlmostEqual(wt_opt["rp.AEP"][0] * 1.0e-6, 23.8821935913, 2) - self.assertAlmostEqual(wt_opt["financese.lcoe"][0] * 1.0e3, 51.6455656178, 2) - self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 4.2027339083, 1) + self.assertAlmostEqual(wt_opt["financese.lcoe"][0] * 1.0e3, 51.6575941808, 2) + self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 4.4814362481, 1) self.assertAlmostEqual(wt_opt["towerse.z_param"][-1], 87.7, 2) def test15MW(self): @@ -41,7 +41,7 @@ def test15MW(self): self.assertAlmostEqual(wt_opt["re.precomp.blade_mass"][0], 69911.6542917299, 1) self.assertAlmostEqual(wt_opt["rp.AEP"][0] * 1.0e-6, 77.76990259134561, 1) self.assertAlmostEqual(wt_opt["financese.lcoe"][0] * 1.0e3, 64.1977089928, 1) - self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 22.8104130129, 1) + self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 25.4598658124, 1) self.assertAlmostEqual(wt_opt["towerse.z_param"][-1], 144.386, 3) def test3p4MW(self): @@ -54,7 +54,7 @@ def test3p4MW(self): self.assertAlmostEqual(wt_opt["re.precomp.blade_mass"][0], 14555.7435212969, 1) self.assertAlmostEqual(wt_opt["rp.AEP"][0] * 1.0e-6, 13.6037235499, 1) self.assertAlmostEqual(wt_opt["financese.lcoe"][0] * 1.0e3, 37.4970510481, 1) - self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 6.606075926116244, 1) + self.assertAlmostEqual(wt_opt["rs.tip_pos.tip_deflection"][0], 6.9771055306, 1) self.assertAlmostEqual(wt_opt["towerse.z_param"][-1], 108.0, 3)