Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Update wisdem #98

Merged
merged 2 commits into from
Mar 25, 2021
Merged

Update wisdem #98

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
ec72583
Squashed 'WISDEM/' changes from 03a0fc61c..bd87e520a
ptrbortolotti Mar 24, 2021
570f034
Merge commit 'ec725839af9ceb5e3135a9ab04687ada38b6f404' into update_w…
ptrbortolotti Mar 24, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 6 additions & 6 deletions WISDEM/wisdem/glue_code/gc_PoseOptimization.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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":
Expand Down
9 changes: 5 additions & 4 deletions WISDEM/wisdem/glue_code/gc_WT_InitModel.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
2 changes: 1 addition & 1 deletion WISDEM/wisdem/inputs/analysis_schema.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
97 changes: 55 additions & 42 deletions WISDEM/wisdem/rotorse/rotor_structure.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):

Expand Down Expand Up @@ -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(
Expand Down Expand Up @@ -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)
Expand All @@ -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)
# -----------------------------------

Expand Down Expand Up @@ -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]
Expand All @@ -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(
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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
Expand Down Expand Up @@ -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)),
Expand All @@ -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
Expand Down Expand Up @@ -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")
self.connect("frame.root_M", "brs.root_M")
2 changes: 1 addition & 1 deletion WISDEM/wisdem/test/test_gluecode/test_gc_modified_yaml.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)


Expand Down
8 changes: 4 additions & 4 deletions WISDEM/wisdem/test/test_gluecode/test_gluecode.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand All @@ -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):
Expand All @@ -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)


Expand Down