Merge 40ea16a into 42ad2bb

wisdem/glue_code/gc_LoadInputs.py

@@ -674,18 +674,18 @@ def recursive_flag(d):
                 blade_opt_options["aero_shape"]["chord"]["n_opt"],
             )
 
-        if not blade_opt_options["aero_shape"]["t/c"]["flag"]:
-            blade_opt_options["aero_shape"]["t/c"]["n_opt"] = self.modeling_options["WISDEM"]["RotorSE"]["n_span"]
-        elif blade_opt_options["aero_shape"]["t/c"]["n_opt"] > self.modeling_options["WISDEM"]["RotorSE"]["n_span"]:
+        if not blade_opt_options["aero_shape"]["rthick"]["flag"]:
+            blade_opt_options["aero_shape"]["rthick"]["n_opt"] = self.modeling_options["WISDEM"]["RotorSE"]["n_span"]
+        elif blade_opt_options["aero_shape"]["rthick"]["n_opt"] > self.modeling_options["WISDEM"]["RotorSE"]["n_span"]:
             raise ValueError("you are attempting to do an analysis using fewer analysis points than control points.")
-        elif blade_opt_options["aero_shape"]["t/c"]["n_opt"] < 4:
-            raise ValueError("Cannot optimize t/c with less than 4 control points along blade span")
-        elif blade_opt_options["aero_shape"]["t/c"]["n_opt"] > self.modeling_options["WISDEM"]["RotorSE"]["n_span"]:
+        elif blade_opt_options["aero_shape"]["rthick"]["n_opt"] < 4:
+            raise ValueError("Cannot optimize rthick with less than 4 control points along blade span")
+        elif blade_opt_options["aero_shape"]["rthick"]["n_opt"] > self.modeling_options["WISDEM"]["RotorSE"]["n_span"]:
             raise ValueError(
                 """Please set WISDEM->RotorSE->n_span in the modeling options yaml larger
-                than aero_shape->t/c->n_opt in the analysis options yaml. n_span and t/c n_opt are """,
+                than aero_shape->rthick->n_opt in the analysis options yaml. n_span and rthick n_opt are """,
                 self.modeling_options["WISDEM"]["RotorSE"]["n_span"],
-                blade_opt_options["aero_shape"]["t/c"]["n_opt"],
+                blade_opt_options["aero_shape"]["rthick"]["n_opt"],
             )
 
         if not blade_opt_options["aero_shape"]["L/D"]["flag"]:
@@ -853,13 +853,14 @@ def write_ontology(self, wt_opt, fname_output):
             self.wt_init["components"]["blade"]["outer_shape_bem"]["pitch_axis"]["values"] = wt_opt[
                 "blade.outer_shape_bem.pitch_axis"
             ].tolist()
+            self.wt_init["components"]["blade"]["outer_shape_bem"]["rthick"] = {}
+            self.wt_init["components"]["blade"]["outer_shape_bem"]["rthick"]["grid"] = wt_opt[
+                "blade.outer_shape_bem.s"
+            ].tolist()
+            self.wt_init["components"]["blade"]["outer_shape_bem"]["rthick"]["values"] = wt_opt[
+                "blade.interp_airfoils.r_thick_interp"
+            ].tolist()
             if self.modeling_options["WISDEM"]["RotorSE"]["inn_af"]:
-                self.wt_init["components"]["blade"]["outer_shape_bem"]["t/c"]["grid"] = wt_opt[
-                    "blade.outer_shape_bem.s"
-                ].tolist()
-                self.wt_init["components"]["blade"]["outer_shape_bem"]["t/c"]["values"] = wt_opt[
-                    "blade.interp_airfoils.r_thick_interp"
-                ].tolist()
                 self.wt_init["components"]["blade"]["outer_shape_bem"]["L/D"]["grid"] = wt_opt[
                     "blade.outer_shape_bem.s"
                 ].tolist()

wisdem/glue_code/gc_PoseOptimization.py

@@ -634,14 +634,14 @@ def set_design_variables(self, wt_opt, wt_init):
                 upper=init_stall_opt[indices] + stall_options["max_increase"],
             )
 
-        t_c_options = blade_opt["aero_shape"]["t/c"]
+        t_c_options = blade_opt["aero_shape"]["rthick"]
         if t_c_options["flag"]:
             n_opt = t_c_options["n_opt"]
             indices = range(t_c_options["index_start"], t_c_options["index_end"])
             s_opt_t_c = np.linspace(0.0, 1.0, n_opt)
             t_c_interpolator = PchipInterpolator(
-                wt_init["components"]["blade"]["outer_shape_bem"]["t/c"]["grid"],
-                wt_init["components"]["blade"]["outer_shape_bem"]["t/c"]["values"],
+                wt_init["components"]["blade"]["outer_shape_bem"]["rthick"]["grid"],
+                wt_init["components"]["blade"]["outer_shape_bem"]["rthick"]["values"],
             )
             init_t_c_opt = t_c_interpolator(s_opt_t_c)
 
@@ -1553,10 +1553,10 @@ def set_initial(self, wt_opt, wt_init):
             init_chord_opt = chord_interpolator(wt_opt["blade.opt_var.s_opt_chord"])
             wt_opt["blade.opt_var.chord_opt"] = init_chord_opt
             if self.modeling["WISDEM"]["RotorSE"]["inn_af"]:
-                wt_opt["inn_af.s_opt_r_thick"] = np.linspace(0.0, 1.0, blade_opt["aero_shape"]["t/c"]["n_opt"])
+                wt_opt["inn_af.s_opt_r_thick"] = np.linspace(0.0, 1.0, blade_opt["aero_shape"]["rthick"]["n_opt"])
                 r_thick_interpolator = PchipInterpolator(
-                    wt_init["components"]["blade"]["outer_shape_bem"]["t/c"]["grid"],
-                    wt_init["components"]["blade"]["outer_shape_bem"]["t/c"]["values"],
+                    wt_init["components"]["blade"]["outer_shape_bem"]["rthick"]["grid"],
+                    wt_init["components"]["blade"]["outer_shape_bem"]["rthick"]["values"],
                 )
                 init_r_thick_opt = r_thick_interpolator(wt_opt["inn_af.s_opt_r_thick"])
                 wt_opt["inn_af.r_thick_opt"] = init_r_thick_opt

wisdem/glue_code/gc_WT_DataStruc.py

@@ -93,11 +93,11 @@ def setup(self):
 
                 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"]["rthick"]["n_opt"])
                 )
                 inn_af.add_output(
                     "r_thick_opt",
-                    val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["t/c"]["n_opt"]),
+                    val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["rthick"]["n_opt"]),
                 )
                 inn_af.add_output(
                     "s_opt_L_D", val=np.ones(opt_options["design_variables"]["blade"]["aero_shape"]["L/D"]["n_opt"])
@@ -197,7 +197,7 @@ def setup(self):
                 ),
             )
             self.connect("airfoils.name", "blade.interp_airfoils.name")
-            self.connect("airfoils.r_thick", "blade.interp_airfoils.r_thick")
+            self.connect("airfoils.r_thick", "blade.interp_airfoils.r_thick_discrete")
             self.connect("airfoils.ac", "blade.interp_airfoils.ac")
             self.connect("airfoils.coord_xy", "blade.interp_airfoils.coord_xy")
             self.connect("airfoils.aoa", "blade.interp_airfoils.aoa")
@@ -621,7 +621,7 @@ def setup(self):
             self.connect("blade.interp_airfoils.cm_interp", "af_3d.cm")
             self.connect("blade.high_level_blade_props.rotor_radius", "af_3d.rotor_radius")
             self.connect("blade.high_level_blade_props.r_blade", "af_3d.r_blade")
-            self.connect("blade.interp_airfoils.r_thick_interp", "af_3d.r_thick_interp")
+            self.connect("blade.interp_airfoils.r_thick_interp", "af_3d.r_thick")
             self.connect("blade.pa.chord_param", "af_3d.chord")
             self.connect("control.rated_TSR", "af_3d.rated_TSR")
         if modeling_options["flags"]["tower"]:
@@ -722,6 +722,7 @@ def setup(self):
 
         # Connections from oute_shape_bem to interp_airfoils
         self.connect("outer_shape_bem.s", "interp_airfoils.s")
+        self.connect("outer_shape_bem.r_thick_yaml_interp", "interp_airfoils.r_thick_yaml")
         self.connect("pa.chord_param", ["interp_airfoils.chord", "compute_coord_xy_dim.chord"])
         self.connect("outer_shape_bem.pitch_axis", ["interp_airfoils.pitch_axis", "compute_coord_xy_dim.pitch_axis"])
         self.connect("opt_var.af_position", "interp_airfoils.af_position")
@@ -742,7 +743,7 @@ def setup(self):
             )
             self.connect("outer_shape_bem.s", "run_inn_af.s")
             self.connect("pa.chord_param", "run_inn_af.chord")
-            self.connect("interp_airfoils.r_thick_interp", "run_inn_af.r_thick_interp_yaml")
+            self.connect("interp_airfoils.r_thick_interp", "run_inn_af.r_thick")
             self.connect("interp_airfoils.cl_interp", "run_inn_af.cl_interp_yaml")
             self.connect("interp_airfoils.cd_interp", "run_inn_af.cd_interp_yaml")
             self.connect("interp_airfoils.cm_interp", "run_inn_af.cm_interp_yaml")
@@ -860,6 +861,9 @@ def setup(self):
             units="m",
             desc="2D array of the coordinates (x,y,z) of the blade reference axis, defined along blade span. The coordinate system is the one of BeamDyn: it is placed at blade root with x pointing the suction side of the blade, y pointing the trailing edge and z along the blade span. A standard configuration will have negative x values (prebend), if swept positive y values, and positive z values.",
         )
+        ivc.add_output(
+            "r_thick_yaml", val=np.zeros(n_span), desc="1D array of the relative thickness values defined along blade span."
+        )
 
         self.add_subsystem(
             "compute_blade_outer_shape_bem",
@@ -896,6 +900,11 @@ def setup(self):
             units="rad",
             desc="1D array of the twist values defined along blade span. The twist is defined positive for negative rotations around the z axis (the same as in BeamDyn).",
         )
+        self.add_input(
+            "r_thick_yaml",
+            val=np.zeros(n_span),
+            desc="1D array of the relative thickness values defined along blade span.",
+        )
         self.add_input(
             "pitch_axis_yaml",
             val=np.zeros(n_span),
@@ -938,6 +947,11 @@ def setup(self):
             val=np.zeros(n_span),
             desc="1D array of the chordwise position of the pitch axis (0-LE, 1-TE), defined along blade span.",
         )
+        self.add_output(
+            "r_thick_yaml_interp",
+            val=np.zeros(n_span),
+            desc="1D array of the relative thickness values defined along blade span.",
+        )
         self.add_output(
             "ref_axis",
             val=np.zeros((n_span, 3)),
@@ -953,6 +967,7 @@ def compute(self, inputs, outputs):
             chord_orig = PchipInterpolator(inputs["s_default"], inputs["chord_yaml"])(nd_span_orig)
             twist_orig = PchipInterpolator(inputs["s_default"], inputs["twist_yaml"])(nd_span_orig)
             pitch_axis_orig = PchipInterpolator(inputs["s_default"], inputs["pitch_axis_yaml"])(nd_span_orig)
+            r_thick_orig = PchipInterpolator(inputs["s_default"], inputs["r_thick_yaml"])(nd_span_orig)
             ref_axis_orig = np.zeros((self.n_span, 3))
             ref_axis_orig[:, 0] = PchipInterpolator(inputs["s_default"], inputs["ref_axis_yaml"][:, 0])(nd_span_orig)
             ref_axis_orig[:, 1] = PchipInterpolator(inputs["s_default"], inputs["ref_axis_yaml"][:, 1])(nd_span_orig)
@@ -980,7 +995,7 @@ def compute(self, inputs, outputs):
             outputs["chord"] = PchipInterpolator(nd_span_orig, chord_orig)(outputs["s"])
             outputs["twist"] = PchipInterpolator(nd_span_orig, twist_orig)(outputs["s"])
             outputs["pitch_axis"] = PchipInterpolator(nd_span_orig, pitch_axis_orig)(outputs["s"])
-
+            outputs["r_thick_yaml_interp"] = PchipInterpolator(nd_span_orig, r_thick_orig)(outputs["s"])
             outputs["ref_axis"][:, 0] = PchipInterpolator(nd_span_orig, ref_axis_orig[:, 0])(outputs["s"])
             outputs["ref_axis"][:, 1] = PchipInterpolator(nd_span_orig, ref_axis_orig[:, 1])(outputs["s"])
             outputs["ref_axis"][:, 2] = PchipInterpolator(nd_span_orig, ref_axis_orig[:, 2])(outputs["s"])
@@ -989,6 +1004,7 @@ def compute(self, inputs, outputs):
             outputs["chord"] = inputs["chord_yaml"]
             outputs["twist"] = inputs["twist_yaml"]
             outputs["pitch_axis"] = inputs["pitch_axis_yaml"]
+            outputs["r_thick_yaml_interp"] = inputs["r_thick_yaml"]
             outputs["ref_axis"] = inputs["ref_axis_yaml"]
 
 
@@ -1033,7 +1049,7 @@ def setup(self):
         # Airfoil properties
         self.add_discrete_input("name", val=n_af * [""], desc="1D array of names of airfoils.")
         self.add_input("ac", val=np.zeros(n_af), desc="1D array of the aerodynamic centers of each airfoil.")
-        self.add_input("r_thick", val=np.zeros(n_af), desc="1D array of the relative thicknesses of each airfoil.")
+        self.add_input("r_thick_discrete", val=np.zeros(n_af), desc="1D array of the relative thicknesses of each airfoil.")
         self.add_input(
             "aoa",
             val=np.zeros(n_aoa),
@@ -1062,6 +1078,11 @@ def setup(self):
             val=np.zeros((n_af, n_xy, 2)),
             desc="3D array of the x and y airfoil coordinates of the n_af airfoils.",
         )
+        self.add_input(
+            "r_thick_yaml",
+            val=np.zeros(n_span),
+            desc="1D array of the relative thicknesses of the blade defined along span.",
+        )
 
         # Polars and coordinates interpolated along span
         self.add_output(
@@ -1111,7 +1132,7 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         for i in range(self.n_af_span):
             for j in range(self.n_af):
                 if self.af_used[i] == discrete_inputs["name"][j]:
-                    r_thick_used[i] = inputs["r_thick"][j]
+                    r_thick_used[i] = inputs["r_thick_discrete"][j]
                     ac_used[i] = inputs["ac"][j]
                     coord_xy_used[i, :, :] = inputs["coord_xy"][j]
                     cl_used[i, :, :, :] = inputs["cl"][j, :, :, :]
@@ -1124,7 +1145,12 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
         spline = PchipInterpolator
         rthick_spline = spline(inputs["af_position"], r_thick_used)
         ac_spline = spline(inputs["af_position"], ac_used)
-        outputs["r_thick_interp"] = rthick_spline(inputs["s"])
+        if np.max(inputs["r_thick_yaml"]) < 1.e-6:
+            rthick_spline = spline(inputs["af_position"], r_thick_used)
+            outputs["r_thick_interp"] = rthick_spline(inputs["s"])
+        else:
+            outputs["r_thick_interp"] = inputs["r_thick_yaml"]
+        ac_spline = spline(inputs["af_position"], ac_used)
         outputs["ac_interp"] = ac_spline(inputs["s"])
 
         # Spanwise interpolation of the profile coordinates with a pchip
@@ -1265,7 +1291,7 @@ def setup(self):
             desc="1D array of the non-dimensional spanwise grid defined along blade axis (0-blade root, 1-blade tip)",
         )
         self.add_input(
-            "r_thick_interp_yaml",
+            "r_thick",
             val=np.zeros(n_span),
             desc="1D array of the relative thicknesses of the blade defined along span.",
         )
@@ -1300,10 +1326,10 @@ 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["rthick"]["n_opt"]))
         self.add_input(
             "r_thick_opt",
-            val=np.ones(aero_shape_opt_options["t/c"]["n_opt"]),
+            val=np.ones(aero_shape_opt_options["rthick"]["n_opt"]),
         )
         self.add_input("s_opt_L_D", val=np.ones(aero_shape_opt_options["L/D"]["n_opt"]))
         self.add_input(
@@ -1371,7 +1397,7 @@ def setup(self):
         self.inn = INN()
 
     def compute(self, inputs, outputs):
-        # Interpolate t/c and L/D from opt grid to full grid
+        # Interpolate rthick and L/D from opt grid to full grid
         spline = PchipInterpolator
         r_thick_spline = spline(inputs["s_opt_r_thick"], inputs["r_thick_opt"])
         r_thick = r_thick_spline(inputs["s"])
@@ -1383,8 +1409,8 @@ def compute(self, inputs, outputs):
         stall_margin = stall_margin_spline(inputs["s"])
 
         # Find indices for start and end of the optimization
-        max_t_c = self.options["rotorse_options"]["inn_af_max_t/c"]
-        min_t_c = self.options["rotorse_options"]["inn_af_min_t/c"]
+        max_t_c = self.options["rotorse_options"]["inn_af_max_rthick"]
+        min_t_c = self.options["rotorse_options"]["inn_af_min_rthick"]
         indices = np.argwhere(np.logical_and(r_thick > min_t_c, r_thick < max_t_c))
         indices = list(np.squeeze(indices))
 
@@ -3254,7 +3280,7 @@ def setup(self):
             desc="Scalar of the rotor radius, defined ignoring prebend and sweep curvatures, and cone and uptilt angles.",
         )
         self.add_input(
-            "r_thick_interp",
+            "r_thick",
             val=np.zeros(n_span),
             desc="1D array of the relative thicknesses of the blade defined along span.",
         )
@@ -3284,7 +3310,7 @@ def compute(self, inputs, outputs):
         cm_corrected = np.zeros((self.n_span, self.n_aoa, self.n_Re, self.n_tab))
         for i in range(self.n_span):
             if (
-                inputs["r_thick_interp"][i] < 0.7 and self.af_correction
+                inputs["r_thick"][i] < 0.7 and self.af_correction
             ):  # Only apply 3D correction to airfoils thinner than 70% to avoid numerical problems at blade root
                 logger.info("3D correction applied to airfoil polars for section " + str(i))
                 for j in range(self.n_Re):

wisdem/glue_code/gc_WT_InitModel.py

@@ -31,7 +31,8 @@ def yaml2openmdao(wt_opt, modeling_options, wt_init, opt_options):
 
     if modeling_options["flags"]["blade"]:
         blade = wt_init["components"]["blade"]
-        wt_opt = assign_blade_values(wt_opt, modeling_options, blade)
+        blade_DV = opt_options['design_variables']['blade']
+        wt_opt = assign_blade_values(wt_opt, modeling_options, blade_DV, blade)
     else:
         blade = {}
 
@@ -106,16 +107,17 @@ def yaml2openmdao(wt_opt, modeling_options, wt_init, opt_options):
     return wt_opt
 
 
-def assign_blade_values(wt_opt, modeling_options, blade):
+def assign_blade_values(wt_opt, modeling_options, blade_DV, blade):
     # Function to assign values to the openmdao group Blade
-    wt_opt = assign_outer_shape_bem_values(wt_opt, modeling_options, blade["outer_shape_bem"])
+    blade_DV_aero = blade_DV['aero_shape']
+    wt_opt = assign_outer_shape_bem_values(wt_opt, modeling_options, blade_DV_aero, blade["outer_shape_bem"])
     wt_opt = assign_internal_structure_2d_fem_values(wt_opt, modeling_options, blade["internal_structure_2d_fem"])
     wt_opt = assign_te_flaps_values(wt_opt, modeling_options, blade)
 
     return wt_opt
 
 
-def assign_outer_shape_bem_values(wt_opt, modeling_options, outer_shape_bem):
+def assign_outer_shape_bem_values(wt_opt, modeling_options, blade_DV_aero, outer_shape_bem):
     # Function to assign values to the openmdao component Blade_Outer_Shape_BEM
 
     nd_span = modeling_options["WISDEM"]["RotorSE"]["nd_span"]
@@ -133,7 +135,16 @@ def assign_outer_shape_bem_values(wt_opt, modeling_options, outer_shape_bem):
     wt_opt["blade.outer_shape_bem.pitch_axis_yaml"] = PchipInterpolator(
         outer_shape_bem["pitch_axis"]["grid"], outer_shape_bem["pitch_axis"]["values"]
     )(nd_span)
-
+    af_opt_flag = blade_DV_aero['af_positions']['flag']
+    if 'rthick' in outer_shape_bem and af_opt_flag == False:
+        # If rthick is defined in input yaml and we are NOT optimizing airfoil positions
+        wt_opt["blade.outer_shape_bem.r_thick_yaml"] = PchipInterpolator(
+            outer_shape_bem["rthick"]["grid"], outer_shape_bem["rthick"]["values"]
+        )(nd_span)
+    elif 'rthick' in outer_shape_bem and af_opt_flag == True:
+        logger.warning('rthick field in input geometry yaml is specified but neglected since you are optimizing airfoil positions')
+    else:
+        logger.warning('rthick field in input geometry yaml not specified. rthick is reconstructed from discrete airfoil positions')
     wt_opt["blade.outer_shape_bem.ref_axis_yaml"][:, 0] = PchipInterpolator(
         outer_shape_bem["reference_axis"]["x"]["grid"], outer_shape_bem["reference_axis"]["x"]["values"]
     )(nd_span)