From e4cb5de23495cd4fd1b2831aa9bb0804ecc007bb Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 1 Nov 2017 12:09:00 +0100
Subject: [PATCH] more fixes
---
wetb/hawc2/{bladeinfo.py => blade.py} | 36 ++++-
wetb/hawc2/mainbody.py | 217 ++++++++++++++------------
wetb/hawc2/tests/test_blade.py | 48 ++++++
wetb/hawc2/tests/test_main_body.py | 7 +-
4 files changed, 194 insertions(+), 114 deletions(-)
rename wetb/hawc2/{bladeinfo.py => blade.py} (90%)
create mode 100644 wetb/hawc2/tests/test_blade.py
diff --git a/wetb/hawc2/bladeinfo.py b/wetb/hawc2/blade.py
similarity index 90%
rename from wetb/hawc2/bladeinfo.py
rename to wetb/hawc2/blade.py
index 91f3d7b..2e323a5 100644
--- a/wetb/hawc2/bladeinfo.py
+++ b/wetb/hawc2/blade.py
@@ -12,6 +12,7 @@ import os
from wetb.utils.geometry import rad
from wetb.hawc2.st_file import StFile
from wetb.hawc2.ae_file import AEFile
+from wetb.hawc2.mainbody import MainBody
@@ -62,7 +63,7 @@ from wetb.hawc2.ae_file import AEFile
# """
-class H2aeroBladeInfo(PCFile, AEFile, AtTimeFile):
+class H2aeroBlade(PCFile, AEFile, AtTimeFile):
"""Provide HAWC2 info about a blade
From AE file:
@@ -302,7 +303,25 @@ class H2aeroBladeInfo(PCFile, AEFile, AtTimeFile):
return self._Cxxx(radius, alpha, 3, ae_set_nr)
-class H2BladeInfo(H2aeroBladeInfo):
+ def plot_xz_geometry(self, plt):
+
+ z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+ plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
+ plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.chord(z) / 2, label='Leading edge')
+ plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.chord(z) / 2, label="Trailing edge")
+ curve_l_nd, x, y, z, twist = self.hawc2_splines()
+ plt.plot(z, x, label='Hawc2spline')
+
+ def plot_yz_geometry(self, plt):
+
+ z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+ plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
+ plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.thickness(z) / 100 * self.chord(z) / 2, label='Suction side')
+ plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.thickness(z) / 100 * self.chord(z) / 2, label="Pressure side")
+ curve_l_nd, x, y, z, twist = self.hawc2_splines()
+ plt.plot(z, y, label='Hawc2spline')
+
+class H2Blade(H2aeroBlade, MainBody):
"""Provide HAWC2 info about a blade
From AE file:
@@ -340,7 +359,10 @@ class H2BladeInfo(H2aeroBladeInfo):
blade_name = htcfile.aero.link[2]
self.mainbody_blade = htcfile.new_htc_structure.get_subsection_by_name(blade_name)
st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
- H2aeroBladeInfo.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
+ MainBody.__init__(self, htcfile, blade_name)
+ elif st_filename and os.path.isfile(st_filename):
+ StFile.__init__(self, st_filename)
+ H2aeroBlade.__init__(self, htcfile, ae_filename=ae_filename, pc_filename=pc_filename, at_time_filename=at_time_filename, blade_name=blade_name)
# def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
#
@@ -353,9 +375,8 @@ class H2BladeInfo(H2aeroBladeInfo):
# if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
# PCFile.__init__(self, pc_filename, ae_filename)
# blade_radius = self.ae_sets[1][-1,0]
-#
- if st_filename and os.path.isfile(st_filename):
- StFile.__init__(self, st_filename)
+#
+
# if os.path.isfile(at_time_filename):
# AtTimeFile.__init__(self, at_time_filename, blade_radius)
# self.curved_length = self.radius_curved_ac()[-1]
@@ -367,6 +388,7 @@ class H2BladeInfo(H2aeroBladeInfo):
# self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"])
#
# self.hawc2_splines_data = self.hawc2_splines()
+
@property
def c2def(self):
if not hasattr(self, "_c2def"):
@@ -374,7 +396,7 @@ class H2BladeInfo(H2aeroBladeInfo):
return self._c2def
#
-# class H2aeroBladeInfo(H2BladeInfo):
+# class H2aeroBlade(H2Blade):
#
# def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
# """
diff --git a/wetb/hawc2/mainbody.py b/wetb/hawc2/mainbody.py
index 0744294..c74a15c 100644
--- a/wetb/hawc2/mainbody.py
+++ b/wetb/hawc2/mainbody.py
@@ -12,14 +12,23 @@ from wetb.hawc2.st_file import StFile
class MainBody():
- def __init__(self, htc_filename, modelpath, body_name):
- self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+ def __init__(self, htcfile, body_name):
+ if isinstance(htcfile, str):
+ htcfile = HTCFile(htcfile)
+ self.htcfile = htcfile
s = htcfile.new_htc_structure
main_bodies = {s[k].name[0]:s[k] for k in s.keys() if s[k].name_ == "main_body"}
self.main_body = main_bodies[body_name]
+ self.main_body = s.get_subsection_by_name(body_name)
self.stFile = StFile(os.path.join(htcfile.modelpath, self.main_body.timoschenko_input.filename[0]))
- self.c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
+ #self.c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
self.concentrated_mass = [cm.values for cm in self.main_body if cm.name_.startswith('concentrated_mass')]
+
+ @property
+ def c2def(self):
+ if not hasattr(self, "_c2def"):
+ self._c2def = np.array([v.values[1:5] for v in self.main_body.c2_def if v.name_ == "sec"])
+ return self._c2def
def plot_xz_geometry(self, plt=None):
if plt is None:
@@ -43,7 +52,7 @@ class MainBody():
plt.xlabel("z")
plt.ylabel("y")
z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
- plt.plot(self.c2def[:, 2], self.c2def[:, 1], label='Center line')
+ plt.plot(self.c2def[:, 2], self.c2def[:, 1], '.-', label='Center line')
plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_e(z), label='Elastic center')
plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_cg(z), label='Mass center')
plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.stFile.y_sh(z), label='Shear center')
@@ -52,106 +61,106 @@ class MainBody():
plt.legend()
-
-
-class BladeData(object):
- def plot_xz_geometry(self, plt):
-
- z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
- plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
- plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.pcFile.chord(z) / 2, label='Leading edge')
- plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.pcFile.chord(z) / 2, label="Trailing edge")
- x, y, z = self.hawc2_splines()
- #plt.plot(z, x, label='Hawc2spline')
-
- def plot_yz_geometry(self, plt):
-
- z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
- plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
- plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label='Suction side')
- plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label="Pressure side")
- x, y, z = self.hawc2_splines()
- #plt.plot(z, y, label='Hawc2spline')
-
- def hawc2_splines(self):
- curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
- curve_z_nd = curve_z / curve_z[-1]
-
- def akima(x, y):
- n = len(x)
- var = np.zeros((n + 3))
- z = np.zeros((n))
- co = np.zeros((n, 4))
- for i in range(n - 1):
- var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
- var[n + 1] = 2 * var[n] - var[n - 1]
- var[n + 2] = 2 * var[n + 1] - var[n]
- var[1] = 2 * var[2] - var[3]
- var[0] = 2 * var[1] - var[2]
-
- for i in range(n):
- wi1 = abs(var[i + 3] - var[i + 2])
- wi = abs(var[i + 1] - var[i])
- if (wi1 + wi) == 0:
- z[i] = (var[i + 2] + var[i + 1]) / 2
- else:
- z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
-
- for i in range(n - 1):
- dx = x[i + 1] - x[i]
- a = (z[i + 1] - z[i]) * dx
- b = y[i + 1] - y[i] - z[i] * dx
- co[i, 0] = y[i]
- co[i, 1] = z[i]
- co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
- co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
- co[n - 1, 0] = y[n - 1]
- co[n - 1, 1] = z[n - 1]
- co[n - 1, 2] = 0
- co[n - 1, 3] = 0
- return co
-
- def coef2spline(s, co):
- x, y = [], []
- for i, c in enumerate(co.tolist()[:-1]):
- p = np.poly1d(c[::-1])
- z = np.linspace(0, s[i + 1] - s[i ], 10)
- x.extend(s[i] + z)
- y.extend(p(z))
- return y
-
- x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
- return x, y, z
-
-
-
-class Blade(MainBody, BladeData):
- def __init__(self, htc_filename, modelpath=None, blade_number=1):
-
- self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
-
- blade_name = [link[2] for link in htcfile.aero if link.name_.startswith('link') and link[0]==blade_number][0]
- MainBody.__init__(self, htc_filename, modelpath, blade_name)
- self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
- os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
-
- def plot_xz_geometry(self, plt=None):
- if plt is None:
- import matplotlib.pyplot as plt
- plt.figure()
-
- MainBody.plot_xz_geometry(self, plt)
- BladeData.plot_xz_geometry(self, plt=plt)
- plt.legend()
-
- def plot_geometry_yz(self, plt=None):
- if plt is None:
- import matplotlib.pyplot as plt
- plt.figure()
-
- BladeData.plot_yz_geometry(self, plt=plt)
- MainBody.plot_yz_geometry(self, plt)
-
+#
+#
+# class BladeData(object):
+# def plot_xz_geometry(self, plt):
+#
+# z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+# plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]), label='Center line')
+# plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.pcFile.chord(z) / 2, label='Leading edge')
+# plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) - self.pcFile.chord(z) / 2, label="Trailing edge")
+# x, y, z = self.hawc2_splines()
+# #plt.plot(z, x, label='Hawc2spline')
+#
+# def plot_yz_geometry(self, plt):
+#
+# z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100)
+# plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]), label='Center line')
+# plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) + self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label='Suction side')
+# plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 1]) - self.pcFile.thickness(z) / 100 * self.pcFile.chord(z) / 2, label="Pressure side")
+# x, y, z = self.hawc2_splines()
+# #plt.plot(z, y, label='Hawc2spline')
+#
+# def hawc2_splines(self):
+# curve_z = np.r_[0, np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))]
+# curve_z_nd = curve_z / curve_z[-1]
+#
+# def akima(x, y):
+# n = len(x)
+# var = np.zeros((n + 3))
+# z = np.zeros((n))
+# co = np.zeros((n, 4))
+# for i in range(n - 1):
+# var[i + 2] = (y[i + 1] - y[i]) / (x[i + 1] - x[i])
+# var[n + 1] = 2 * var[n] - var[n - 1]
+# var[n + 2] = 2 * var[n + 1] - var[n]
+# var[1] = 2 * var[2] - var[3]
+# var[0] = 2 * var[1] - var[2]
+#
+# for i in range(n):
+# wi1 = abs(var[i + 3] - var[i + 2])
+# wi = abs(var[i + 1] - var[i])
+# if (wi1 + wi) == 0:
+# z[i] = (var[i + 2] + var[i + 1]) / 2
+# else:
+# z[i] = (wi1 * var[i + 1] + wi * var[i + 2]) / (wi1 + wi)
+#
+# for i in range(n - 1):
+# dx = x[i + 1] - x[i]
+# a = (z[i + 1] - z[i]) * dx
+# b = y[i + 1] - y[i] - z[i] * dx
+# co[i, 0] = y[i]
+# co[i, 1] = z[i]
+# co[i, 2] = (3 * var[i + 2] - 2 * z[i] - z[i + 1]) / dx
+# co[i, 3] = (z[i] + z[i + 1] - 2 * var[i + 2]) / dx ** 2
+# co[n - 1, 0] = y[n - 1]
+# co[n - 1, 1] = z[n - 1]
+# co[n - 1, 2] = 0
+# co[n - 1, 3] = 0
+# return co
+#
+# def coef2spline(s, co):
+# x, y = [], []
+# for i, c in enumerate(co.tolist()[:-1]):
+# p = np.poly1d(c[::-1])
+# z = np.linspace(0, s[i + 1] - s[i ], 10)
+# x.extend(s[i] + z)
+# y.extend(p(z))
+# return y
+#
+# x, y, z = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(3)]
+# return x, y, z
+
+
+
+# class Blade(MainBody, BladeData):
+# def __init__(self, htc_filename, modelpath=None, blade_number=1):
+#
+# self.htcfile = htcfile = HTCFile(htc_filename, modelpath)
+#
+# blade_name = [link[2] for link in htcfile.aero if link.name_.startswith('link') and link[0]==blade_number][0]
+# MainBody.__init__(self, htc_filename, modelpath, blade_name)
+# self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]),
+# os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]))
+#
+# def plot_xz_geometry(self, plt=None):
+# if plt is None:
+# import matplotlib.pyplot as plt
+# plt.figure()
+#
+# MainBody.plot_xz_geometry(self, plt)
+# BladeData.plot_xz_geometry(self, plt=plt)
+# plt.legend()
+#
+# def plot_geometry_yz(self, plt=None):
+# if plt is None:
+# import matplotlib.pyplot as plt
+# plt.figure()
+#
+# BladeData.plot_yz_geometry(self, plt=plt)
+# MainBody.plot_yz_geometry(self, plt)
+#
diff --git a/wetb/hawc2/tests/test_blade.py b/wetb/hawc2/tests/test_blade.py
new file mode 100644
index 0000000..f5956b9
--- /dev/null
+++ b/wetb/hawc2/tests/test_blade.py
@@ -0,0 +1,48 @@
+'''
+Created on 3. maj 2017
+
+@author: mmpe
+'''
+import os
+import unittest
+
+from wetb.hawc2.blade import H2Blade
+import numpy as np
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
+class Test(unittest.TestCase):
+
+
+ def testBladeInfo(self):
+ bi = H2Blade(tfp + "simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc")
+ if 0:
+ import matplotlib.pyplot as plt
+ print (dir(bi))
+ #print (bi.radius_s())
+ plt.plot(bi.radius_s(), bi.twist())
+ plt.plot(bi.c2def[:,2], bi.c2def[:,3])
+ x = np.linspace(0,1,1000)
+ plt.plot(bi.blade_radius*x, bi.c2nd(x)[:,3])
+ plt.show()
+
+
+ def testBladeInfo_AE(self):
+ bi = H2Blade(None, tfp + "NREL_5MW_ae.txt", tfp + "NREL_5MW_pc.txt")
+ self.assertEqual(bi.thickness(32), 23.78048780487805)
+ self.assertEqual(bi.chord(32), 3.673)
+ self.assertEqual(bi.pc_set_nr(32), 1)
+
+ def testBladeInfo_PC_AE(self):
+ bi = H2Blade(None, tfp + "NREL_5MW_ae.txt", tfp + "NREL_5MW_pc.txt")
+ self.assertEqual(bi.CL(36, 10), 1.358)
+ self.assertEqual(bi.CD(36, 10), 0.0255)
+ self.assertEqual(bi.CM(36, 10), -0.1103)
+
+# def test_curved_length2radius(self):
+# bi = H2BladeInfo(tfp + "simulation_setup/DTU10MWRef6.0/htc/DTU_10MW_RWT.htc")
+# np.testing.assert_array_equal(bi.xyztwist(86.4979, curved_length=True), bi.xyztwist(86.3655))
+
+
+if __name__ == "__main__":
+ #import sys;sys.argv = ['', 'Test.testBladeInfo']
+ unittest.main()
\ No newline at end of file
diff --git a/wetb/hawc2/tests/test_main_body.py b/wetb/hawc2/tests/test_main_body.py
index 4ddfec2..04b7a93 100644
--- a/wetb/hawc2/tests/test_main_body.py
+++ b/wetb/hawc2/tests/test_main_body.py
@@ -5,13 +5,14 @@ Created on 01/08/2016
'''
import unittest
import os
-from wetb.hawc2.mainbody import MainBody, Blade
+from wetb.hawc2.mainbody import MainBody
+from wetb.hawc2.blade import H2Blade
tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
class TestMainBody(unittest.TestCase):
def test_MainBody(self):
- mainbody = MainBody(tfp+"htcfiles/test.htc", "../", "towertop")
+ mainbody = MainBody(tfp+"htcfiles/test.htc", "blade1")
if 0:
import matplotlib.pyplot as plt
plt.figure()
@@ -21,7 +22,7 @@ class TestMainBody(unittest.TestCase):
plt.show()
def test_Blade(self):
- blade = Blade(tfp+"htcfiles/test.htc", "../")
+ blade = H2Blade(tfp+"htcfiles/test.htc")
if 0:
import matplotlib.pyplot as plt
plt.figure()
--
GitLab