From 724f70d64c949c867ff9a9ecea6b418cd6454476 Mon Sep 17 00:00:00 2001 From: "Mads M. Pedersen" <mmpe@dtu.dk> Date: Thu, 24 Nov 2016 10:35:40 +0100 Subject: [PATCH] Added bladeData.py + test --- wetb/hawc2/bladeData.py | 142 +++++++++++++++++++++++++++++ wetb/hawc2/tests/test_bladeData.py | 29 ++++++ 2 files changed, 171 insertions(+) create mode 100644 wetb/hawc2/bladeData.py create mode 100644 wetb/hawc2/tests/test_bladeData.py diff --git a/wetb/hawc2/bladeData.py b/wetb/hawc2/bladeData.py new file mode 100644 index 00000000..046e9e0b --- /dev/null +++ b/wetb/hawc2/bladeData.py @@ -0,0 +1,142 @@ +''' +Created on 01/08/2016 + +@author: MMPE +''' +from wetb.hawc2.pc_file import PCFile +from wetb.hawc2.ae_file import AEFile +from wetb.hawc2.htc_file import HTCFile +import os +import numpy as np +from wetb.hawc2.st_file import StFile + +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 H2BladeData(BladeData): + def __init__(self, htc_filename, modelpath): + self.htcfile = htcfile = HTCFile(htc_filename, modelpath) + self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]), + os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])) + s = htcfile.new_htc_structure + blade_name = htcfile.aero.link[2] + mainbodies = [s[k] for k in s.keys() if s[k].name_ == "main_body"] + + blade_main_body = [mb for mb in mainbodies if mb.name[0] == blade_name][0] + self.stFile = StFile(os.path.join(htcfile.modelpath, blade_main_body.timoschenko_input.filename[0])) + self.c2def = np.array([v.values[1:5] for v in blade_main_body.c2_def if v.name_ == "sec"]) + + def plot_geometry(self, plt=None): + if plt is None: + import matplotlib.pyplot as plt + + BladeData.plot_xz_geometry(self, plt=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]) + self.stFile.x_e(z), label='Elastic center') + plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.stFile.x_cg(z), label='Mass center') + plt.plot(z, np.interp(z, self.c2def[:, 2], self.c2def[:, 0]) + self.stFile.x_sh(z), label='Shear center') + + plt.legend() + plt.figure() + BladeData.plot_yz_geometry(self, plt=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]) + 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') + plt.legend() + + plt.show() + + + + +class H2AeroBladeData(BladeData): + def __init__(self, htc_filename, modelpath): + self.htcfile = htcfile = HTCFile(htc_filename, modelpath) + self.pcFile = PCFile(os.path.join(htcfile.modelpath, htcfile.aero.pc_filename[0]), + os.path.join(htcfile.modelpath, htcfile.aero.ae_filename[0])) + self.c2def = np.array([v.values[1:5] for v in htcfile.blade_c2_def if v.name_ == "sec"]) + + def plot_geometry(self, plt=None): + if plt is None: + import matplotlib.pyplot as plt + + BladeData.plot_xz_geometry(self, plt=plt) + z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100) + plt.legend() + plt.figure() + BladeData.plot_yz_geometry(self, plt=plt) + z = np.linspace(self.c2def[0, 2], self.c2def[-1, 2], 100) + plt.legend() + + plt.show() + + diff --git a/wetb/hawc2/tests/test_bladeData.py b/wetb/hawc2/tests/test_bladeData.py new file mode 100644 index 00000000..de3b13f4 --- /dev/null +++ b/wetb/hawc2/tests/test_bladeData.py @@ -0,0 +1,29 @@ +''' +Created on 01/08/2016 + +@author: MMPE +''' +import unittest +import os +from wetb.hawc2.bladeData import H2AeroBladeData, H2BladeData +tfp = os.path.join(os.path.dirname(__file__), 'test_files/') # test file path +class TestBladeData(unittest.TestCase): + + def testBladedata(self): + #bd = H2AeroBladeData(tfp + "h2aero_tb/htc/h2aero.htc", '../') + #bd = H2AeroBladeData(r'C:\mmpe\programming\python\Phd\pitot_tube\tests\test_files\swt36_107\h2a\htc\templates/swt3.6_107.htc', "../../") + #print (bd.pcFile.chord(36)) + #bd = H2BladeData(r"C:\mmpe\HAWC2\models\SWT3.6-107\original_data\SWT4.0-130\dlc12_wsp06_wdir010_s18002.htc", ".") + bd = H2BladeData(r"C:\mmpe\HAWC2\models\SWT3.6-107\htc/stat6/stat6_10.0_0.htc", "../../") + #bd = H2BladeData(r"C:\mmpe\HAWC2\models\NREL 5MW reference wind turbine_play/htc/NREL_5MW_reference_wind_turbine_heli_step.htc", "../") + if 1: + bd.plot_geometry() + + + + + + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'Test.testBladedata'] + unittest.main() -- GitLab