From 734dd08569f849578d61caacf28b62da0536fa2a Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Wed, 1 Nov 2017 11:42:30 +0100
Subject: [PATCH] new pc file without ae file. Old pc-ae file integrated in
bladeinfo
---
wetb/hawc2/bladeinfo.py | 134 ++++++++++++++++++++++++-------
wetb/hawc2/pc_file.py | 34 ++++----
wetb/hawc2/tests/test_pc_file.py | 19 ++---
3 files changed, 128 insertions(+), 59 deletions(-)
diff --git a/wetb/hawc2/bladeinfo.py b/wetb/hawc2/bladeinfo.py
index e89bd9a..91f3d7b 100644
--- a/wetb/hawc2/bladeinfo.py
+++ b/wetb/hawc2/bladeinfo.py
@@ -11,6 +11,7 @@ from scipy.interpolate.interpolate import interp1d
import os
from wetb.utils.geometry import rad
from wetb.hawc2.st_file import StFile
+from wetb.hawc2.ae_file import AEFile
@@ -61,7 +62,7 @@ from wetb.hawc2.st_file import StFile
# """
-class H2aeroBladeInfo(PCFile, AtTimeFile):
+class H2aeroBladeInfo(PCFile, AEFile, AtTimeFile):
"""Provide HAWC2 info about a blade
From AE file:
@@ -86,36 +87,37 @@ class H2aeroBladeInfo(PCFile, AtTimeFile):
"""
- def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
-
- if isinstance(htcfile, str):
- assert htcfile.lower().endswith('.htc')
- htcfile = HTCFile(htcfile)
- self.htcfile = htcfile
- blade_name = blade_name or htcfile.aero.link[2]
- #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])
-
+ def __init__(self, htcfile=None, ae_filename=None, pc_filename=None, at_time_filename=None, blade_name=None):
+ if htcfile:
+ if isinstance(htcfile, str):
+ assert htcfile.lower().endswith('.htc')
+ htcfile = HTCFile(htcfile)
+ self.htcfile = htcfile
+ blade_name = blade_name or htcfile.aero.link[2]
+ 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])
+ self.hawc2_splines_data = self.hawc2_splines()
+
#mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"]
#self.mainbody_blade = [mb for mb in mainbodies if mb.name[0] == blade_name][0]
if os.path.isfile(pc_filename) and os.path.isfile(ae_filename):
- PCFile.__init__(self, pc_filename, ae_filename)
+ AEFile.__init__(self, ae_filename)
+ PCFile.__init__(self, pc_filename)
blade_radius = self.ae_sets[1][-1,0]
- if os.path.isfile(at_time_filename):
+ if at_time_filename and 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
+ self.curved_length = None
#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.hawc2_splines_data = self.hawc2_splines()
+
@property
@@ -222,6 +224,83 @@ class H2aeroBladeInfo(PCFile, AtTimeFile):
assert np.all(l>=self.c2def[0,2]) and np.all(l<=self.c2def[-1,2])
return self.c2nd(l/self.c2def[-1, 2])
+ def _Cxxx(self, radius, alpha, column, ae_set_nr=1):
+ thickness = self.thickness(radius, ae_set_nr)
+ pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+ thicknesses, profiles = self.pc_sets[pc_set_nr]
+ index = np.searchsorted(thicknesses, thickness)
+ if index == 0:
+ index = 1
+
+ Cx0, Cx1 = profiles[index - 1:index + 1]
+ Cx0 = np.interp(alpha, Cx0[:, 0], Cx0[:, column])
+ Cx1 = np.interp(alpha, Cx1[:, 0], Cx1[:, column])
+ th0, th1 = thicknesses[index - 1:index + 1]
+ return Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
+
+ def _CxxxH2(self, radius, alpha, column, ae_set_nr=1):
+ thickness = self.thickness(radius, ae_set_nr)
+ pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+ thicknesses, profiles = self.pc_sets[pc_set_nr]
+ index = np.searchsorted(thicknesses, thickness)
+ if index == 0:
+ index = 1
+
+ Cx0, Cx1 = profiles[index - 1:index + 1]
+
+ Cx0 = np.interp(np.arange(360), Cx0[:,0]+180, Cx0[:,column])
+ Cx1 = np.interp(np.arange(360), Cx1[:,0]+180, Cx1[:,column])
+ #Cx0 = np.interp(alpha, Cx0[:, 0], Cx0[:, column])
+ #Cx1 = np.interp(alpha, Cx1[:, 0], Cx1[:, column])
+ th0, th1 = thicknesses[index - 1:index + 1]
+ cx = Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
+ return np.interp(alpha+180, np.arange(360), cx)
+
+
+
+ def CL(self, radius, alpha, ae_set_nr=1):
+ """Lift coefficient
+
+ Parameters
+ ---------
+ radius : float
+ radius [m]
+ alpha : float
+ Angle of attack [deg]
+ ae_set_nr : int optional
+ Aerdynamic set number, default is 1
+
+ Returns
+ -------
+ Lift coefficient : float
+ """
+ return self._Cxxx(radius, alpha, 1, ae_set_nr)
+
+
+ def CL_H2(self, radius, alpha, ae_set_nr=1):
+ return self._CxxxH2(radius, alpha, 1, ae_set_nr)
+
+ def CD(self, radius, alpha, ae_set_nr=1):
+ """Drag coefficient
+
+ Parameters
+ ---------
+ radius : float
+ radius [m]
+ alpha : float
+ Angle of attack [deg]
+ ae_set_nr : int optional
+ Aerdynamic set number, default is 1
+
+ Returns
+ -------
+ Drag coefficient : float
+ """
+ return self._Cxxx(radius, alpha, 2, ae_set_nr)
+
+ def CM(self, radius, alpha, ae_set_nr=1):
+ return self._Cxxx(radius, alpha, 3, ae_set_nr)
+
class H2BladeInfo(H2aeroBladeInfo):
"""Provide HAWC2 info about a blade
@@ -247,18 +326,20 @@ class H2BladeInfo(H2aeroBladeInfo):
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
+ def __init__(self, htcfile=None, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
+ if htcfile is not 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]
+ 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 = 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)
# def __init__(self, htcfile, ae_filename=None, pc_filename=None, at_time_filename=None, st_filename=None, blade_name=None):
@@ -268,13 +349,12 @@ class H2BladeInfo(H2aeroBladeInfo):
# 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):
+ 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)
diff --git a/wetb/hawc2/pc_file.py b/wetb/hawc2/pc_file.py
index 301a5b8..84b15cc 100644
--- a/wetb/hawc2/pc_file.py
+++ b/wetb/hawc2/pc_file.py
@@ -16,7 +16,7 @@ standard_library.install_aliases()
from wetb.hawc2.ae_file import AEFile
import numpy as np
-class PCFile(AEFile):
+class PCFile(object):
"""Read HAWC2 PC (profile coefficients) file
examples
@@ -37,9 +37,7 @@ class PCFile(AEFile):
>>> pcfile.CM(36,10) # CM at radius=36m and AOA=10deg
-0.1103
"""
- def __init__(self, filename, ae_filename):
- AEFile.__init__(self, ae_filename)
-
+ def __init__(self, filename):
with open (filename) as fid:
lines = fid.readlines()
nsets = int(lines[0].split()[0])
@@ -61,9 +59,7 @@ class PCFile(AEFile):
lptr += n_rows
self.pc_sets[nset] = (np.array(thicknesses), profiles)
- def _Cxxx(self, radius, alpha, column, ae_set_nr=1):
- thickness = self.thickness(radius, ae_set_nr)
- pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+ def _Cxxx(self, thickness, alpha, column, ae_set_nr=1, pc_set_nr=1):
thicknesses, profiles = self.pc_sets[pc_set_nr]
index = np.searchsorted(thicknesses, thickness)
if index == 0:
@@ -75,9 +71,7 @@ class PCFile(AEFile):
th0, th1 = thicknesses[index - 1:index + 1]
return Cx0 + (Cx1 - Cx0) * (thickness - th0) / (th1 - th0)
- def _CxxxH2(self, radius, alpha, column, ae_set_nr=1):
- thickness = self.thickness(radius, ae_set_nr)
- pc_set_nr = self.pc_set_nr(radius, ae_set_nr)
+ def _CxxxH2(self, thickness, alpha, column, ae_set_nr=1,pc_set_nr=1):
thicknesses, profiles = self.pc_sets[pc_set_nr]
index = np.searchsorted(thicknesses, thickness)
if index == 0:
@@ -95,13 +89,13 @@ class PCFile(AEFile):
- def CL(self, radius, alpha, ae_set_nr=1):
+ def CL(self, thickness, alpha, ae_set_nr=1,pc_set_nr=1):
"""Lift coefficient
Parameters
---------
- radius : float
- radius [m]
+ thickness : float
+ thickness [5]
alpha : float
Angle of attack [deg]
ae_set_nr : int optional
@@ -111,13 +105,13 @@ class PCFile(AEFile):
-------
Lift coefficient : float
"""
- return self._Cxxx(radius, alpha, 1, ae_set_nr)
+ return self._Cxxx(thickness, alpha, 1, ae_set_nr, pc_set_nr)
- def CL_H2(self, radius, alpha, ae_set_nr=1):
- return self._CxxxH2(radius, alpha, 1, ae_set_nr)
+ def CL_H2(self, thickness, alpha, ae_set_nr=1, pc_set_nr=1):
+ return self._CxxxH2(thickness, alpha, 1, ae_set_nr, pc_set_nr)
- def CD(self, radius, alpha, ae_set_nr=1):
+ def CD(self, thickness, alpha, ae_set_nr=1, pc_set_nr=1):
"""Drag coefficient
Parameters
@@ -133,10 +127,10 @@ class PCFile(AEFile):
-------
Drag coefficient : float
"""
- return self._Cxxx(radius, alpha, 2, ae_set_nr)
+ return self._Cxxx(thickness, alpha, 2, ae_set_nr, pc_set_nr)
- def CM(self, radius, alpha, ae_set_nr=1):
- return self._Cxxx(radius, alpha, 3, ae_set_nr)
+ def CM(self, thickness, alpha, ae_set_nr=1, pc_set_nr=1):
+ return self._Cxxx(thickness, alpha, 3, ae_set_nr,pc_set_nr)
if __name__ == "__main__":
pc = PCFile(r"C:\mmpe\Projects\inflow\Hawc2aero_setup/data/Hawc_pc.b52", r"C:\mmpe\Projects\inflow\Hawc2aero_setup/data/S36_ae_h2.001")
diff --git a/wetb/hawc2/tests/test_pc_file.py b/wetb/hawc2/tests/test_pc_file.py
index e5a7027..0ac9e5b 100644
--- a/wetb/hawc2/tests/test_pc_file.py
+++ b/wetb/hawc2/tests/test_pc_file.py
@@ -8,6 +8,7 @@ from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
+from wetb.hawc2.ae_file import AEFile
standard_library.install_aliases()
import os
import unittest
@@ -24,19 +25,13 @@ class TestPCFile(unittest.TestCase):
unittest.TestCase.setUp(self)
self.testfilepath = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path
-
def test_PCFile_ae(self):
- pc = PCFile(self.testfilepath + "NREL_5MW_pc.txt", self.testfilepath + "NREL_5MW_ae.txt")
- self.assertEqual(pc.thickness(32), 23.78048780487805)
- self.assertEqual(pc.chord(32), 3.673)
- self.assertEqual(pc.pc_set_nr(32), 1)
-
- def test_PCFile2(self):
- pc = PCFile(self.testfilepath + "NREL_5MW_pc.txt", self.testfilepath + "NREL_5MW_ae.txt")
- self.assertEqual(pc.CL(36, 10), 1.358)
- self.assertEqual(pc.CD(36, 10), 0.0255)
- self.assertEqual(pc.CM(36, 10), -0.1103)
-
+ pc = PCFile(self.testfilepath + "NREL_5MW_pc.txt")
+ ae = AEFile(self.testfilepath + "NREL_5MW_ae.txt")
+ thickness = ae.thickness(36)
+ self.assertEqual(pc.CL(thickness, 10), 1.358)
+ self.assertEqual(pc.CD(thickness, 10), 0.0255)
+ self.assertEqual(pc.CM(thickness, 10), -0.1103)
if __name__ == "__main__":
#import sys;sys.argv = ['', 'Test.testName']
--
GitLab