B/we mode

WISDEM/docs/inputs/analysis_schema.rst

@@ -90,7 +90,10 @@ Blade twist as a design variable by adding or subtracting radians from the initi
     *Default* = False
 
 :code:`inverse` : Boolean
-    Words TODO?
+    When set to True, the twist is defined inverting the 
+    blade-element momentum equations to achieve a desired margin to stall, 
+    which is defined among the constraints.
+    :code:`flag` and :code:`inverse` cannot be simultaneously be set to True
 
     *Default* = False
 
@@ -172,7 +175,8 @@ Adjust airfoil positions along the blade span.
 
 :code:`af_start` : Integer
     Index of airfoil where the optimization can start shifting airfoil
-    position. The airfoil at blade tip is always locked.
+    position. The airfoil at blade tip is always locked. It is advised 
+    to keep the airfoils close to blade root locked.
 
     *Default* = 4
 

WISDEM/examples/02_reference_turbines/IEA-15-240-RWT.yaml

@@ -347,7 +347,7 @@ components:
         spinner_material: glass_uni
     nacelle:
         drivetrain:
-            uptilt_angle: 0.10471975511965977 # 6 deg
+            uptilt: 0.10471975511965977 # 6 deg
             distance_tt_hub: 5.614
             overhang: 12.0313
             drag_coefficient: 0.5

WISDEM/examples/02_reference_turbines/IEA-3p4-130-RWT.yaml

@@ -280,7 +280,7 @@ components:
     nacelle:
         drivetrain:
             diameter: 3.0
-            uptilt_angle: 0.08726 # 5 deg
+            uptilt: 0.08726 # 5 deg
             distance_tt_hub: 2.
             distance_hub2mb: 1.912
             distance_mb2mb: 0.368

WISDEM/examples/02_reference_turbines/analysis_options.yaml

@@ -144,9 +144,9 @@ constraints:
             flag: False
             margin: 1.4175
         rail_transport:
-            flag: True
+            flag: False
             8_axle: False
-            4_axle: True
+            4_axle: False
         stall:
             flag: False         # Constraint on minimum stall margin
             margin: 0.05233 # Value of minimum stall margin in [rad]

WISDEM/examples/02_reference_turbines/nrel5mw.yaml

@@ -196,7 +196,7 @@ components:
     nacelle:
         drivetrain:
             diameter: 3.0
-            uptilt_angle: 0.08726
+            uptilt: 0.08726
             distance_tt_hub: 2.3
             distance_hub2mb: 1.912
             distance_mb2mb: 0.368

WISDEM/examples/03_blade/blade.yaml

@@ -338,7 +338,7 @@ components:
     nacelle:
         drivetrain:
             diameter: 3.0
-            uptilt_angle: 0.10471975511965977 # 6 deg
+            uptilt: 0.10471975511965977 # 6 deg
             distance_tt_hub: 3.0
             distance_hub2mb: 1.912
             distance_mb2mb: 0.368

WISDEM/examples/05_tower_monopile/modeling_options.yaml

@@ -7,14 +7,16 @@ DriveSE:
     flag: False
 TowerSE:             # Options of TowerSE module
     flag: True
-    nLC: 2           # Number of design load cases 
-    wind: PowerWind  # Wind used 
-    gamma_f: 1.35    # Safety factor for fatigue loads 
-    gamma_m: 1.3     # Safety factor for material properties 
+    nLC: 2           # Number of design load cases
+    wind: PowerWind  # Wind used
+    gamma_f: 1.35    # Safety factor for fatigue loads
+    gamma_m: 1.3     # Safety factor for material properties
     gamma_n: 1.0     # Safety factor for ...
     gamma_b: 1.1     # Safety factor for ...
-    gamma_fatigue: 1.755  # Safety factor for fatigue loads 
-    buckling_length: 30    # Buckling parameter 
+    gamma_fatigue: 1.755  # Safety factor for fatigue loads
+    buckling_length: 30    # Buckling parameter
+    soil_springs: True
+    gravity_foundation: False
     frame3dd:
         shear: True
         geom: True

WISDEM/examples/05_tower_monopile/monopile_direct.py

@@ -39,6 +39,10 @@
 modeling_options["TowerSE"]["buckling_length"] = 30.0
 modeling_options["flags"]["monopile"] = True
 
+# Monopile foundation
+modeling_options["TowerSE"]["soil_springs"] = True
+modeling_options["TowerSE"]["gravity_foundation"] = False
+
 # safety factors
 modeling_options["TowerSE"]["gamma_f"] = 1.35
 modeling_options["TowerSE"]["gamma_m"] = 1.3

WISDEM/examples/05_tower_monopile/tower_direct.py

@@ -30,6 +30,10 @@
 modeling_options["TowerSE"]["buckling_length"] = 30.0
 modeling_options["flags"]["monopile"] = False
 
+# Monopile foundation only
+modeling_options["TowerSE"]["soil_springs"] = False
+modeling_options["TowerSE"]["gravity_foundation"] = False
+
 # safety factors
 modeling_options["TowerSE"]["gamma_f"] = 1.35
 modeling_options["TowerSE"]["gamma_m"] = 1.3
@@ -98,10 +102,6 @@
 prob["tower_outfitting_factor"] = 1.07
 prob["yaw"] = 0.0
 
-# offshore specific
-prob["G_soil"] = 140e6
-prob["nu_soil"] = 0.4
-
 # material properties
 prob["E_mat"] = 210e9 * np.ones((n_materials, 3))
 prob["G_mat"] = 80.8e9 * np.ones((n_materials, 3))

WISDEM/examples/09_weis_floating/nrel5mw-semi_oc4.yaml

@@ -197,7 +197,7 @@ components:
     nacelle:
         drivetrain:
             diameter: 3.0
-            uptilt_angle: 0.08726
+            uptilt: 0.08726
             distance_tt_hub: 2.3
             distance_hub2mb: 1.912
             distance_mb2mb: 0.368

WISDEM/examples/09_weis_floating/nrel5mw-spar_oc3.yaml

@@ -197,7 +197,7 @@ components:
     nacelle:
         drivetrain:
             diameter: 3.0
-            uptilt_angle: 0.08726
+            uptilt: 0.08726
             distance_tt_hub: 2.3
             distance_hub2mb: 1.912
             distance_mb2mb: 0.368

WISDEM/wisdem/ccblade/ccblade_component.py

@@ -3,6 +3,7 @@
 import numpy as np
 import wisdem.ccblade._bem as _bem
 from wisdem.commonse.csystem import DirectionVector
+from scipy.interpolate import PchipInterpolator
 
 
 cosd = lambda x: np.cos(np.deg2rad(x))
@@ -532,6 +533,8 @@ def compute(self, inputs, outputs, discrete_inputs, discrete_outputs):
                 )
 
         if self.options["opt_options"]["design_variables"]["blade"]["aero_shape"]["twist"]["inverse"]:
+            if self.options["opt_options"]["design_variables"]["blade"]["aero_shape"]["twist"]["flag"]:
+                raise Exception("Twist cannot be simultaneously optimized and set to be defined inverting the BEM equations. Please check your analysis options yaml.")
             # Find cl and cd for max efficiency
             cl = np.zeros(self.n_span)
             cd = np.zeros(self.n_span)

WISDEM/wisdem/commonse/environment.py

@@ -9,15 +9,10 @@
 
 from __future__ import print_function
 
-import sys
-import math
-
 import numpy as np
 import openmdao.api as om
 from scipy.optimize import brentq
-
-from .constants import gravity
-from .utilities import hstack, vstack
+from wisdem.commonse.constants import gravity
 
 # TODO CHECK
 
@@ -379,10 +374,10 @@ def compute(self, inputs, outputs):
         h = inputs["Hsig_wave"]
 
         # circular frequency
-        omega = 2.0 * math.pi / inputs["Tsig_wave"]
+        omega = 2.0 * np.pi / inputs["Tsig_wave"]
 
         # compute wave number from dispersion relationship
-        k = brentq(lambda k: omega ** 2 - gravity * k * math.tanh(d * k), 0, 1e3 * omega ** 2 / gravity)
+        k = brentq(lambda k: omega ** 2 - gravity * k * np.tanh(d * k), 0, 1e3 * omega ** 2 / gravity)
         self.k = k
         outputs["phase_speed"] = omega / k
 
@@ -427,7 +422,7 @@ def compute_partials(self, inputs, J):
         z = inputs["z"]
         d = inputs["z_surface"] - z_floor
         h = inputs["Hsig_wave"]
-        omega = 2.0 * math.pi / inputs["Tsig_wave"]
+        omega = 2.0 * np.pi / inputs["Tsig_wave"]
         k = self.k
         z_rel = z - inputs["z_surface"]
 
@@ -497,27 +492,31 @@ class TowerSoil(om.ExplicitComponent):
         Spring stiffness (x, theta_x, y, theta_y, z, theta_z)
     """
 
+    def initialize(self):
+        self.options.declare("npts", default=1)
+
     def setup(self):
-        super(TowerSoil, self).setup()
+        npts = self.options["npts"]
 
         # variable
         self.add_input("d0", 1.0, units="m")
-        self.add_input("depth", 1.0, units="m")
+        self.add_input("depth", 0.0, units="m")
 
         # inputeter
         self.add_input("G", 140e6, units="Pa")
         self.add_input("nu", 0.4)
         self.add_input("k_usr", -1 * np.ones(6), units="N/m")
 
-        self.add_output("k", np.zeros(6), units="N/m")
+        self.add_output("z_k", np.zeros(npts), units="N/m")
+        self.add_output("k", np.zeros((npts, 6)), units="N/m")
 
         self.declare_partials("k", ["d0", "depth"])
 
     def compute(self, inputs, outputs):
 
         G = inputs["G"]
         nu = inputs["nu"]
-        h = inputs["depth"]
+        h = np.linspace(inputs["depth"], 0.0, self.options["npts"])
         r0 = 0.5 * inputs["d0"]
 
         # vertical
@@ -533,17 +532,18 @@ def compute(self, inputs, outputs):
         k_thetax = 8.0 * G * r0 ** 3 * eta / (3.0 * (1.0 - nu))
 
         # torsional
-        k_phi = 16.0 * G * r0 ** 3 / 3.0
+        k_phi = 16.0 * G * r0 ** 3 * np.ones(h.size) / 3.0
 
-        outputs["k"] = np.array([k_x, k_thetax, k_x, k_thetax, k_z, k_phi]).flatten()
+        outputs["k"] = np.c_[k_x, k_thetax, k_x, k_thetax, k_z, k_phi]
+        outputs["z_k"] = -h
         ind = np.nonzero(inputs["k_usr"] >= 0.0)[0]
-        outputs["k"][ind] = inputs["k_usr"][ind]
+        outputs["k"][:, ind] = inputs["k_usr"][np.newaxis, ind]
 
     def compute_partials(self, inputs, J):
 
         G = inputs["G"]
         nu = inputs["nu"]
-        h = inputs["depth"]
+        h = np.linspace(inputs["depth"], 0.0, self.options["npts"])
         r0 = 0.5 * inputs["d0"]
 
         # vertical
@@ -574,13 +574,13 @@ def compute_partials(self, inputs, J):
         dkphi_dr0 = 16.0 * G * 3 * r0 ** 2 / 3.0
         dkphi_dh = 0.0
 
-        dk_dr0 = np.array([dkx_dr0, dkthetax_dr0, dkx_dr0, dkthetax_dr0, dkz_dr0, dkphi_dr0])
+        dk_dr0 = np.c_[dkx_dr0, dkthetax_dr0, dkx_dr0, dkthetax_dr0, dkz_dr0, dkphi_dr0]
         # dk_dr0[inputs['rigid']] = 0.0
-        dk_dh = np.array([dkx_dh, dkthetax_dh, dkx_dh, dkthetax_dh, dkz_dh, dkphi_dh])
+        dk_dh = np.c_[dkx_dh, dkthetax_dh, dkx_dh, dkthetax_dh, dkz_dh, dkphi_dh]
         # dk_dh[inputs['rigid']] = 0.0
 
         J["k", "d0"] = 0.5 * dk_dr0
         J["k", "depth"] = dk_dh
         ind = np.nonzero(inputs["k_usr"] >= 0.0)[0]
-        J["k", "d0"][ind] = 0.0
-        J["k", "depth"][ind] = 0.0
+        J["k", "d0"][:, ind] = 0.0
+        J["k", "depth"][:, ind] = 0.0

WISDEM/wisdem/commonse/fileIO.py

@@ -78,6 +78,7 @@ def save_data(fname, prob, npz_file=True, mat_file=True, xls_file=True):
             data["description"].append(var_dict[k][1]["desc"])
         df = pd.DataFrame(data)
         df.to_excel(froot + ".xlsx", index=False)
+        df.to_csv(froot + ".csv", index=False)
 
 
 def load_data(fname, prob):

WISDEM/wisdem/commonse/vertical_cylinder.py

@@ -1,13 +1,11 @@
 import numpy as np
 import openmdao.api as om
-
-from wisdem.commonse import gravity, eps, NFREQ
 import wisdem.commonse.frustum as frustum
-import wisdem.commonse.manufacturing as manufacture
-from wisdem.commonse.utilization_constraints import hoopStressEurocode, hoopStress
 import wisdem.commonse.utilities as util
 import wisdem.pyframe3dd.pyframe3dd as pyframe3dd
-
+import wisdem.commonse.manufacturing as manufacture
+from wisdem.commonse import NFREQ, eps, gravity
+from wisdem.commonse.utilization_constraints import hoopStress, hoopStressEurocode
 
 RIGID = 1e30
 NREFINE = 3
@@ -61,6 +59,7 @@ class CylinderDiscretization(om.ExplicitComponent):
     def initialize(self):
         self.options.declare("nPoints")
         self.options.declare("nRefine", default=NREFINE)
+        self.options.declare("nPin", default=0)
 
     def setup(self):
         nPoints = self.options["nPoints"]
@@ -90,12 +89,14 @@ def compute(self, inputs, outputs):
 
         nRefine = int(np.round(self.options["nRefine"]))
         z_param = float(inputs["foundation_height"]) + np.r_[0.0, np.cumsum(inputs["section_height"].flatten())]
+
         # Have to regine each element one at a time so that we preserve input nodes
         z_full = np.array([])
         for k in range(z_param.size - 1):
             zref = np.linspace(z_param[k], z_param[k + 1], nRefine + 1)
             z_full = np.append(z_full, zref)
         z_full = np.unique(z_full)
+
         outputs["z_full"] = z_full
         outputs["d_full"] = np.interp(z_full, z_param, inputs["diameter"])
         z_section = 0.5 * (z_full[:-1] + z_full[1:])
@@ -355,8 +356,8 @@ class CylinderFrame3DD(om.ExplicitComponent):
         6-degree polynomial coefficients of mode shapes in the x-direction
     y_mode_shapes : numpy array[NFREQ2, 5]
         6-degree polynomial coefficients of mode shapes in the x-direction
-    top_deflection : float, [m]
-        Deflection of cylinder top in yaw-aligned +x direction
+    cylinder_deflection : numpy array[npts], [m]
+        Deflection of cylinder nodes in yaw-aligned +x direction
     Fz_out : numpy array[npts-1], [N]
         Axial foce in vertical z-direction in cylinder structure.
     Vx_out : numpy array[npts-1], [N]
@@ -463,7 +464,7 @@ def setup(self):
         self.add_output("y_mode_shapes", val=np.zeros((NFREQ2, 5)))
         self.add_output("x_mode_freqs", val=np.zeros(NFREQ2))
         self.add_output("y_mode_freqs", val=np.zeros(NFREQ2))
-        self.add_output("top_deflection", val=0.0, units="m")
+        self.add_output("cylinder_deflection", val=np.zeros(npts), units="m")
         self.add_output("Fz_out", val=np.zeros(npts - 1), units="N")
         self.add_output("Vx_out", val=np.zeros(npts - 1), units="N")
         self.add_output("Vy_out", val=np.zeros(npts - 1), units="N")
@@ -623,7 +624,9 @@ def compute(self, inputs, outputs):
         outputs["y_mode_shapes"] = mshapes_y[:NFREQ2, :]
 
         # deflections due to loading (from cylinder top and wind/wave loads)
-        outputs["top_deflection"] = displacements.dx[iCase, n - 1]  # in yaw-aligned direction
+        outputs["cylinder_deflection"] = np.sqrt(
+            displacements.dx[iCase, :] ** 2 + displacements.dy[iCase, :] ** 2
+        )  # in yaw-aligned direction
 
         # shear and bending, one per element (convert from local to global c.s.)
         Fz = forces.Nx[iCase, 1::2]
@@ -635,8 +638,14 @@ def compute(self, inputs, outputs):
         Mxx = -forces.Mzz[iCase, 1::2]
 
         # Record total forces and moments
-        outputs["base_F"] = -1.0 * np.r_[-forces.Vz[iCase, 0], forces.Vy[iCase, 0], forces.Nx[iCase, 0]]
-        outputs["base_M"] = -1.0 * np.r_[-forces.Mzz[iCase, 0], forces.Myy[iCase, 0], forces.Txx[iCase, 0]]
+        base_idx = 2 * int(inputs["kidx"].max())
+        outputs["base_F"] = (
+            -1.0 * np.r_[-forces.Vz[iCase, base_idx], forces.Vy[iCase, base_idx], forces.Nx[iCase, base_idx]]
+        )
+        outputs["base_M"] = (
+            -1.0 * np.r_[-forces.Mzz[iCase, base_idx], forces.Myy[iCase, base_idx], forces.Txx[iCase, base_idx]]
+        )
+
         outputs["Fz_out"] = Fz
         outputs["Vx_out"] = Vx
         outputs["Vy_out"] = Vy