Skip to content
Snippets Groups Projects
Commit 5adc82b2 authored by Mads M. Pedersen's avatar Mads M. Pedersen
Browse files

restructured bladeinfo.py

parent 7d628546
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
wetb/hawc2/bladeinfo.py
+ 199
− 102
View file @ 5adc82b2
...@@ -14,53 +14,54 @@ from wetb.hawc2.st_file import StFile ...@@ -14,53 +14,54 @@ from wetb.hawc2.st_file import StFile
class BladeInfo(object): # class BladeInfo(object):
def twist(self, radius=None): # def twist(self, radius=None):
"""Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n # """Aerodynamic twist [deg]. Negative when leading edge is twisted towards wind(opposite to normal definition)\n
#
# Parameters
# ---------
# radius : float or array_like, optional
# radius [m] of interest
# If None (default) the twist of all points are returned
#
# Returns
# -------
# twist [deg] : float
# """
# raise NotImplementedError()
#
# def chord(self, radius=None):
# """Aerodynamic chord
#
# Parameters
# ---------
# radius : float or array_like, optional
# radius [m] of interest
# If None (default) the twist of all points are returned
#
# Returns
# -------
# chord [m] : float
# """
# raise NotImplementedError()
#
#
#
# def c2nd(self, radius_nd):
# """Return center line position
#
# Parameters
# ---------
# radius_nd : float or array_like, optional
# non dimensional radius
#
# Returns
# -------
# x,y,z : float
# """
Parameters
---------
radius : float or array_like, optional
radius [m] of interest
If None (default) the twist of all points are returned
Returns
-------
twist [deg] : float
"""
raise NotImplementedError()
def chord(self, radius=None):
"""Aerodynamic chord
Parameters
---------
radius : float or array_like, optional
radius [m] of interest
If None (default) the twist of all points are returned
Returns
-------
chord [m] : float
"""
raise NotImplementedError()
class H2aeroBladeInfo(PCFile, AtTimeFile):
def c2nd(self, radius_nd):
"""Return center line position
Parameters
---------
radius_nd : float or array_like, optional
non dimensional radius
Returns
-------
x,y,z : float
"""
class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
"""Provide HAWC2 info about a blade """Provide HAWC2 info about a blade
From AE file: From AE file:
...@@ -85,45 +86,64 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile): ...@@ -85,45 +86,64 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
""" """
def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None): def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
if isinstance(htcfile, str): if isinstance(htcfile, str):
assert htcfile.lower().endswith('.htc') assert htcfile.lower().endswith('.htc')
htcfile = HTCFile(htcfile) htcfile = HTCFile(htcfile)
self.htcfile = htcfile
blade_name = blade_name or htcfile.aero.link[2] blade_name = blade_name or htcfile.aero.link[2]
s = htcfile.new_htc_structure #s = htcfile.new_htc_structure
at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat")) at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]) pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0]) ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"] #mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0] #self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
st_filename = st_filename or os.path.join(htcfile.modelpath, mainbody_blade.timoschenko_input.filename[0])
if os.path.isfile(pc_filename) and os.path.isfile(ae_filename): if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
PCFile.__init__(self, pc_filename, ae_filename) PCFile.__init__(self, pc_filename, ae_filename)
blade_radius = self.ae_sets[1][-1,0] blade_radius = self.ae_sets[1][-1,0]
if os.path.isfile(st_filename):
StFile.__init__(self, st_filename)
if os.path.isfile(at_time_filename): if os.path.isfile(at_time_filename):
AtTimeFile.__init__(self, at_time_filename, blade_radius) AtTimeFile.__init__(self, at_time_filename, blade_radius)
self.curved_length = self.radius_curved_ac()[-1]
else:
self.c2def = np.array([v.values[1:5] for v in mainbody_blade.c2_def if v.name_ == "sec"]) raise NotImplementedError
#z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
#self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
#self.c2nd = lambda x : interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')(np.max([np.min([x, np.ones_like(x)], 0), np.zeros_like(x) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
x, y, z, twist = self.hawc2_splines()
def interpolate(r):
r = (np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
return np.array([np.interp(r, np.array(z) / z[-1], xyz) for xyz in [x,y,z, twist]]).T
self.c2nd = interpolate #lambda r : interp1d(np.array(z) / z[-1], np.array([x, y, z]).T, axis=0, kind=1)(np.max([np.min([r, np.ones_like(r)], 0), np.zeros_like(r) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
self.hawc2_splines_data = self.hawc2_splines()
@property
def c2def(self):
if not hasattr(self, "_c2def"):
self._c2def = np.array([v.values[1:5] for v in self.htcfile.blade_c2_def if v.name_ == "sec"])
return self._c2def
def c2nd(self, l_nd, curved_length=True):
curve_l_nd, x, y, z, twist = self.hawc2_splines_data
if curved_length:
l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd) + self.c2def[0, 2] / self.c2def[-1, 2]], 0))
return np.array([np.interp(l_nd, curve_l_nd, xyz) for xyz in [x,y,z, twist]]).T
else:
l_nd = (np.max([np.min([l_nd, np.ones_like(l_nd)], 0), np.zeros_like(l_nd)], 0))
return np.array([np.interp(l_nd, z/z[-1], xyz) for xyz in [x,y,z, twist]]).T
def c2(self, l, curved_length=True):
if curved_length:
L = self.curved_length
else:
L = self.c2def[-1,2]
return self.c2nd(l/L, curved_length)
def hawc2_splines(self): 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_l = 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] curve_l_nd = curve_l / curve_l[-1]
def akima(x, y): def akima(x, y):
n = len(x) n = len(x)
...@@ -166,10 +186,11 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile): ...@@ -166,10 +186,11 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
z = np.linspace(0, s[i + 1] - s[i ], 10, endpoint=i >= co.shape[0] - 2) z = np.linspace(0, s[i + 1] - s[i ], 10, endpoint=i >= co.shape[0] - 2)
x.extend(s[i] + z) x.extend(s[i] + z)
y.extend(p(z)) y.extend(p(z))
return y return x,y
x, y, z, twist = [coef2spline(curve_z_nd, akima(curve_z_nd, self.c2def[:, i])) for i in range(4)] x, y, z, twist = [coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, i]))[1] for i in range(4)]
return x, y, z, twist curve_l_nd = coef2spline(curve_l_nd, akima(curve_l_nd, self.c2def[:, 0]))[0]
return curve_l_nd, x, y, z, twist
def xyztwist(self, l=None, curved_length=False): def xyztwist(self, l=None, curved_length=False):
"""Return splined x,y,z and twist """Return splined x,y,z and twist
...@@ -188,44 +209,120 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile): ...@@ -188,44 +209,120 @@ class H2BladeInfo(BladeInfo, PCFile, AtTimeFile):
x,y,z,twist x,y,z,twist
""" """
if l is None: if l is None:
return self.c2def[:, 3] return self.c2def[:, :3]
else: else:
r_nd = np.linspace(0,1,100) r_nd = np.linspace(0,1,1000)
if curved_length: if curved_length:
curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1))) #curved_length = np.cumsum(np.sqrt((np.diff(self.c2nd(np.linspace(0,1,100)),1,0)[:,:3]**2).sum(1)))
assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1]) return self.c2nd(l/self.radius_curved_ac()[-1])
return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1]) # z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
# assert np.all(l>=curved_length[0]) and np.all(l<=curved_length[-1])
# return self.c2nd(r_nd[np.argmin(np.abs(curved_length-l))+1])
else: else:
assert np.all(l>=self.c2def[0,2]) and np.all(l<=self.c2def[-1,2]) assert np.all(l>=self.c2def[0,2]) and np.all(l<=self.c2def[-1,2])
return self.c2nd(l/self.c2def[-1, 2]) return self.c2nd(l/self.c2def[-1, 2])
class H2aeroBladeInfo(H2BladeInfo): class H2BladeInfo(H2aeroBladeInfo):
"""Provide HAWC2 info about a blade
def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
""" From AE file:
at_time_filename: file name of at time file containing twist and chord data - chord(radius=None, set_nr=1):
""" - thickness(radius=None, set_nr=1)
PCFile.__init__(self, pc_filename, ae_filename) - radius_ae(radius=None, set_nr=1)
blade_radius = self.ae_sets[1][-1,0]
AtTimeFile.__init__(self, at_time_filename, blade_radius) From PC file
- CL(radius, alpha, ae_set_nr=1)
assert('twist' in self.attribute_names) - CD(radius, alpha, ae_set_nr=1)
htcfile = HTCFile(htc_filename) - CM(radius, alpha, ae_set_nr=1)
self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
#self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
ac_radii = self.ac_radius()
c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
#self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
def c2nd(self, r_nd):
r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
r_nd_max = np.ones_like(r_nd)
r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
return self._c2nd(r_nd)
From at_time_filename
- attribute_names
- xxx(radius=None, curved_length=None) # xxx for each attribute name
- radius_curved_ac(radius=None) # Curved length of nearest/all aerodynamic calculation points
From ST file
- radius_st(radius=None, mset=1, set=1)
- xxx(radius=None, mset=1, set=1) # xxx for each of r, m, x_cg,y_cg, ri_x, ri_y, xs, ys, E, G, Ix, Iy, K, kx, ky, A, pitch, xe, ye
with template
"""
def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
if isinstance(htcfile, str):
htcfile = HTCFile(htcfile)
s = htcfile.new_htc_structure
# at_time_filename = at_time_filename or ("output_at_time" in htcfile and os.path.join(htcfile.modelpath, htcfile.output_at_time.filename[0] + ".dat"))
# pc_filename = pc_filename or os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0])
# ae_filename = ae_filename or os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])
#
mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
if blade_name is None:
blade_name = htcfile.aero.link[2]
self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
H2aeroBladeInfo.__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):
#
# if isinstance(htcfile, str):
# assert htcfile.lower().endswith('.htc')
# htcfile = HTCFile(htcfile)
#
# blade_name = blade_name or htcfile.aero.link[2]
st_filename = st_filename or os.path.join(htcfile.modelpath, self.mainbody_blade.timoschenko_input.filename[0])
#
# 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 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]
# else:
# raise NotImplementedError
# #z_nd = (np.cos(np.linspace(np.pi, np.pi*2,len(curved_length)-1))+1)/2
# #self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
#
# 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"):
self._c2def = np.array([v.values[1:5] for v in self.mainbody_blade.c2_def if v.name_ == "sec"])
return self._c2def
#
# class H2aeroBladeInfo(H2BladeInfo):
#
# def __init__(self, at_time_filename, ae_filename, pc_filename, htc_filename):
# """
# at_time_filename: file name of at time file containing twist and chord data
# """
# PCFile.__init__(self, pc_filename, ae_filename)
# blade_radius = self.ae_sets[1][-1,0]
# AtTimeFile.__init__(self, at_time_filename, blade_radius)
#
# assert('twist' in self.attribute_names)
# htcfile = HTCFile(htc_filename)
#
#
# self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"])
# #self._c2nd = interp1d(self.c2def[:, 2] / self.c2def[-1, 2], self.c2def[:, :3], axis=0, kind='cubic')
#
# ac_radii = self.radius_curved_ac()
# self.hawc2_splines_data = self.hawc2_splines()
# self.curved_length = np.cumsum(np.sqrt(np.sum(np.diff(self.c2def[:, :3], 1, 0) ** 2, 1)))[-1]
#
# # c2_axis_length = np.r_[0, np.cumsum(np.sqrt((np.diff(self.c2def[:, :3], 1, 0) ** 2).sum(1)))]
# # self._c2nd = interp1d(c2_axis_length / c2_axis_length[-1], self.c2def[:, :3], axis=0, kind='cubic')
# #self._c2nd = interp1d(self.c2def[:,2]/self.c2def[-1,2], self.c2def[:, :3], axis=0, kind='cubic')
# #
# # def c2nd(self, r_nd):
# # r_nd_min = np.zeros_like(r_nd) + self.c2def[0, 2] / self.c2def[-1, 2]
# # r_nd_max = np.ones_like(r_nd)
# # r_nd = np.max([np.min([r_nd, r_nd_max], 0), r_nd_min], 0)
# # return self._c2nd(r_nd)
#
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment