From 7b7614263675254825017404ce3887a7d8787a42 Mon Sep 17 00:00:00 2001 From: verd_he Date: Thu, 31 Mar 2022 12:31:36 +0200 Subject: [PATCH] Added tool interfaces for the stability analysis use case for alaska/Wind, HAWC2 and HAWCStab2. These interfaces were developed by Nordex (alaska) and LUH (HAWC2, HAWCStab2) --- .../tool_interfaces/HS2ASCIIReader.py | 169 ++++++++ .../tool_interfaces/alaska_interface.py | 309 ++++++++++++++ .../tool_interfaces/hawc2_interface.py | 392 ++++++++++++++++++ .../tool_interfaces/hawcstab2_interface.py | 68 +++ .../tool_interfaces/header_st_file_hawc2.txt | 4 + .../header_st_file_hawc2_6x6.txt | 7 + 6 files changed, 949 insertions(+) create mode 100644 use_cases/stability_analysis/tool_interfaces/HS2ASCIIReader.py create mode 100644 use_cases/stability_analysis/tool_interfaces/alaska_interface.py create mode 100644 use_cases/stability_analysis/tool_interfaces/hawc2_interface.py create mode 100644 use_cases/stability_analysis/tool_interfaces/hawcstab2_interface.py create mode 100644 use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2.txt create mode 100644 use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2_6x6.txt diff --git a/use_cases/stability_analysis/tool_interfaces/HS2ASCIIReader.py b/use_cases/stability_analysis/tool_interfaces/HS2ASCIIReader.py new file mode 100644 index 0000000..6a4601a --- /dev/null +++ b/use_cases/stability_analysis/tool_interfaces/HS2ASCIIReader.py @@ -0,0 +1,169 @@ +import numpy as np + + +class HS2DataContainer(object): + """ data container for HS2 linearization result data """ + + def __init__(self): + super(HS2DataContainer, self).__init__() + + self.names = [] + self.frequency = 0.0 # np array: windspeed , frequency + self.damping = 0.0 # np array: windspeed , damping ratio + self.realpart = 0.0 # np array: windspeed , real part + + +class HS2Data(object): + """ + This is a class for operating with HAWCStab2 linearization result data + + Parameters + ---------- + filenamecmb : string + Path to .cmb HAWCStab2 output file + filenameamp : string + Path to .amp HAWCStab2 output file + filenameop : string + Path to .op/.opt HAWCStab2 output file + """ + def __init__(self, filenamecmb, filenameamp, filenameop): + super(HS2Data, self).__init__() + + self.filenamecmb = filenamecmb + self.filenameamp = filenameamp + self.filenameop = filenameop + self.modes = HS2DataContainer() + + def read_cmb_data(self): + """ + Reads and parses the HS2 result cmb data + + Raises + ------ + OSError + If .cmb file could not be found + + Notes + ----- + There are three blocks for freq./damping/real part. + """ + + # read file + try: + hs2cmd = np.loadtxt(self.filenamecmb,skiprows=1,dtype='float') + + myshape=np.shape(hs2cmd) + if np.shape(myshape)==(1,): + num_windspeeds = 1 + num_modes = int((myshape[0]-1)/3) + else: + num_windspeeds = int(myshape[0]) + num_modes = int((myshape[1]-1)/3) + self.modes.frequency = np.zeros([num_modes,num_windspeeds,2]) + self.modes.damping = np.zeros([num_modes,num_windspeeds,2]) + self.modes.realpart = np.zeros([num_modes,num_windspeeds,2]) + + if np.shape(myshape)==(1,): + for i_mode in range(0,num_modes): + self.modes.frequency[i_mode,:,0] = hs2cmd[0 ] + self.modes.frequency[i_mode,:,1] = hs2cmd[i_mode+ 1] + self.modes.damping [i_mode,:,0] = hs2cmd[0 ] + self.modes.damping [i_mode,:,1] = hs2cmd[i_mode+ 61] + self.modes.realpart [i_mode,:,0] = hs2cmd[0 ] + self.modes.realpart [i_mode,:,1] = hs2cmd[i_mode+121] + else: + for i_mode in range(0,num_modes): + self.modes.frequency[i_mode,:,0] = hs2cmd[:,0 ] + self.modes.frequency[i_mode,:,1] = hs2cmd[:,i_mode+ 1] + self.modes.damping [i_mode,:,0] = hs2cmd[:,0 ] + self.modes.damping [i_mode,:,1] = hs2cmd[:,i_mode+ 61] + self.modes.realpart [i_mode,:,0] = hs2cmd[:,0 ] + self.modes.realpart [i_mode,:,1] = hs2cmd[:,i_mode+121] + + print('INFO: HS2 campbell data loaded successfully:') + print(' - %4i modes'%num_modes) + print(' - %4i wind speeds'%num_windspeeds) + + except OSError: + print('ERROR: HAWCStab2 cmb file %s not found! Abort!'%self.filenamecmb) + raise + + def read_amp_data(self): + """ + Reads and parses the HS2 result amp data + + Raises + ------ + OSError + If .amp file could not be found + + Notes + ----- + 1st column is wind speed, each mode has 30 column + + # Mode number: 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 + # Column num.: 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 + # Wind speed TWR x [m] phase [deg] TWR y [m] phase [deg] TWR yaw [rad] phase [deg] SFT x [m] phase [deg] SFT y [m] phase [deg] SFT tor [rad] phase [deg] Sym edge [m] phase [deg] BW edge [m] phase [deg] FW edge [m] phase [deg] Sym flap [m] phase [deg] BW flap [m] phase [deg] FW flap [m] phase [deg]Sym tors [rad] phase [deg] BW tors [rad] phase [deg] FW tors [rad] phase [deg] + """ + sensor_list = ['TWR SS','TWR FA','TWR yaw','SFT x','SFT y','SFT tor','Sym edge','BW edge','FW edge','Sym flap','BW flap','FW flap','Sym tors','BW tors','FW tors'] + num_sensors = len(sensor_list) + + # read file + try: + hs2amp = np.loadtxt(self.filenameamp,skiprows=5,dtype='float') + + myshape=np.shape(hs2amp) + if np.shape(myshape)==(1,): + num_windspeeds = 1 + num_modes = int((myshape[0]-1)/num_sensors/2) + else: + num_windspeeds = int(myshape[0]) + num_modes = int((myshape[1]-1)/num_sensors/2) + amp_data = np.zeros([num_windspeeds,num_sensors,num_modes]) + + if np.shape(myshape)==(1,): + for i in range(0,num_modes): + for j in range(0,num_sensors): + amp_index = i*num_sensors*2+1+2*j + amp_data[:,j,i] = hs2amp[amp_index] + else: + for i in range(0,num_modes): + for j in range(0,num_sensors): + amp_index = i*num_sensors*2+1+2*j + amp_data[:,j,i] = hs2amp[:,amp_index] + + # determine dominant DOF per mode + for i in range(0,num_modes): + mean_DOF = np.mean(amp_data[:,:,i],axis=0) + self.modes.names.append(sensor_list[np.argmax(mean_DOF)]) + + # override first tower mode + if self.modes.names[2] == sensor_list[0]: + self.modes.names[1] = sensor_list[1] + + print('INFO: HS2 amplitude data loaded successfully:') + print(' - %4i modes'%num_modes) + print(' - %4i wind speeds'%num_windspeeds) + except OSError: + print('ERROR: HAWCStab2 amp file %s not found! Abort!'%self.filenameamp) + raise + + def read_op_data(self): + """ + Reads the operational data from HS2 + + Raises + ------ + OSError + If .op/.opt file could not be found + """ + try: + self.HS2_opdata = np.loadtxt(self.filenameop,skiprows=1,dtype='float') + except OSError: + print('ERROR: HAWCStab2 op file %s not found! Abort!'%self.filenameamp) + raise + + def read_data(self): + self.read_cmb_data() + self.read_amp_data() + self.read_op_data() \ No newline at end of file diff --git a/use_cases/stability_analysis/tool_interfaces/alaska_interface.py b/use_cases/stability_analysis/tool_interfaces/alaska_interface.py new file mode 100644 index 0000000..703963b --- /dev/null +++ b/use_cases/stability_analysis/tool_interfaces/alaska_interface.py @@ -0,0 +1,309 @@ +""" +@author: Sarah Mueller +@date: 17.03.2021 +""" + +# global lib +import sys +import os +import xml.dom.minidom +import shutil +import numpy as np +import subprocess +from copy import deepcopy +from scipy import interpolate + +# local lib +sys.path.append(r"D:\12_Projekte\QuexUS\03_UQ\coding") +from compressed_int_loader import CompressedIntLoader +from tool_interfaces.simulation_model_interface import SimulationModel, SimulationError + + +class AlaskaModel(SimulationModel): + """ + alaska/Wind model interface + + Parameters + ---------- + iteration_run_directory : string + Path to directory of this framework iteration + config : dict + User-specified settings for this tool + + Attributes + ---------- + GEBT_template + fef_template + alaska_exe + run_directory + alaska_model + alaska_temp + """ + def __init__(self, iteration_run_directory, config): + # super(AlaskaModel, self).__init__() # super class constructor is called + self.GEBT_template = config['path_GEBT_template'] + self.fef_template = config['path_fef_template'] + self.alaska_exe = config['path_alaska_exe'] + self.run_directory = iteration_run_directory + self.alaska_model = config['path_alaska_model'] + self.alaska_temp = config['path_alaska_templates'] + + def create_simulation(self, blade_data): + """ + Setup modified alaska model + + Parameters + ---------- + blade_data : dict + Dictionary with structural blade data + + Notes + ----- + blade data: alaska specific stiffness and mass matrices + blade_data['alaska_6x6']['stiffness_matrix'][blade section id][row][column] + blade_data['alaska_6x6']['mass_matrix'][blade section id][row][column] + """ + + # copy alaska model + src_ala = os.path.join(self.alaska_temp, 'alaska_model_UQ_temp') + dest_ala = os.path.join(self.run_directory, 'alaska_model') + shutil.copytree(src_ala, dest_ala) + shutil.copyfile(os.path.join(self.alaska_temp, 'Includes.dat'), os.path.join(self.run_directory, 'Includes.dat')) + + # write parameters + src_param = os.path.join(self.alaska_temp, 'Parameters_UQ_temp.dat') + dest_param = os.path.join(self.run_directory, 'alaska_model', 'Parameters', 'Parameters.dat') + shutil.copyfile(src_param, dest_param) + with open(dest_param,'r') as f: + f_param = f.read() + with open(dest_param,'w') as f: + f.write(f_param % os.path.join(self.run_directory, 'alaska_model')) + + # write load cases + src_LC = os.path.join(self.alaska_temp, 'LC_UQ_temp.xml') + dest_LC = os.path.join(self.run_directory, 'alaska_model', 'LoadCases', 'LC_UQ.xml') + shutil.copyfile(src_LC, dest_LC) + with open(dest_LC, 'r') as f: + f_LC = f.read() + with open(dest_LC, 'w') as f: + f.write(f_LC % os.path.join(self.run_directory, 'alaska_model')) + + # write GEBT xml + #self.write_GEBT_xml(self.GEBT_template, blade_data) + # write fef + self.write_fef(self.fef_template, blade_data) + + def write_GEBT_xml(self, file, blade_data): + """ + function that writes GEBT.xml file + """ + # parsing the xml file + xml_doc = xml.dom.minidom.parse(file) + # get the xml data tree + xml_tree = xml_doc.documentElement + + # change stiffness and mass matrix + for section_id, cross_section in enumerate(xml_tree.getElementsByTagName("CROSS_SECTION")): + # write stiffness data into the xml + k_txt = np.array2string(blade_data['alaska_6x6_element_ref']['stiffness_matrix'][section_id]) + k_txt = k_txt.replace("[", "") + k_txt = k_txt.replace("]]", "]") + k_txt = k_txt.replace("]", "\n") + cross_section.getElementsByTagName("FLEXIBILITY")[0].firstChild.replaceWholeText("\n " + k_txt) + # write mass data into the xml + m_txt = np.array2string(blade_data['alaska_6x6_element_ref']['mass_matrix'][section_id]) + m_txt = m_txt.replace("[", "") + m_txt = m_txt.replace("]]", "]") + m_txt = m_txt.replace("]", "\n") + cross_section.getElementsByTagName("MASS")[0].firstChild.replaceWholeText("\n " + m_txt) + + # write xml + with open(file.replace(".xml", "_UQ_temp.xml"), "w")as f: + xml_string = xml_tree.toprettyxml() + f.write(self.remove_blank_lines_xml(xml_string)) + + def write_fef(self, file, blade_data): + """ + function that writes .fef file + """ + # parsing the xml file + xml_doc = xml.dom.minidom.parse(file) + # get the xml data tree + xml_tree = xml_doc.documentElement + + # define permutation vector + perm_vec = np.array([0, 4, 5, 3, 1, 2]) + stiffmat_perm = np.empty([6, 6]) + massmat_perm = np.empty([6, 6]) + + # get indices for the upper-triangle + i_upper = np.triu_indices(6, k=0) + # get number of blade sections + i_blade_section = len(blade_data['alaska_6x6_element_ref']['stiffness_matrix']) + # initialize string + matrices_str = "" + # calculate averaged matrices (alaska setting: BladeDescriptionType = mean) + for i in range(0, i_blade_section - 1): + stiffmat_ave = (blade_data['alaska_6x6_element_ref']['stiffness_matrix'][i] + blade_data['alaska_6x6_element_ref']['stiffness_matrix'][i + 1]) / 2 + massmat_ave = (blade_data['alaska_6x6_element_ref']['mass_matrix'][i] + blade_data['alaska_6x6_element_ref']['mass_matrix'][i + 1]) / 2 + # apply permutation (not needed for mass matrix) + for j in range(0, 6): + for k in range(0, 6): + stiffmat_perm[j, k] = deepcopy(stiffmat_ave[perm_vec[j], perm_vec[k]]) # use a deepcopy + massmat_perm[j, k] = deepcopy(massmat_ave[perm_vec[j], perm_vec[k]]) # use a deepcopy + # get upper triangle and convert to string + stiffmat_fef = np.array2string(stiffmat_perm[i_upper]) # use permutated matrix + massmat_fef = np.array2string(massmat_ave[i_upper]) # use averaged matrix + # format to one line + matrices_str = (matrices_str + "\n" + "M" + np.str(i+1) + " " + massmat_fef + "\n" + "C" + np.str(i+1) + " " + stiffmat_fef) + + # remove bracket + matrices_str = matrices_str.replace("[", "") + matrices_str = matrices_str.replace("]", "") + # change matrices + xml_tree.getElementsByTagName("TABLE")[6].firstChild.data = ("\n " + matrices_str + "\n") + + # write fef + head, tail = os.path.split(file) + with open(file.replace(".fef", "_UQ_temp.fef"), "w")as f: + xml_string = xml_tree.toprettyxml() + f.write(self.remove_blank_lines_xml(xml_string)) + + src = os.path.join(head, 'FEBlade_UQ_temp.fef') + dest = os.path.join(self.run_directory, 'alaska_model', 'SubModels', 'Blade_%i', 'FEBlade.fef') + for blade_id in range(1, 4): + shutil.copyfile(src, dest % blade_id) + # check if fef file is not empty + for blade_id in range(1, 4): + if os.stat(dest % blade_id).st_size == 0: + print('error-fef: empty fef in blade ') + print(blade_id) + print(self.run_directory) + shutil.copyfile(src, dest % blade_id) + + def remove_blank_lines_xml(self, xml_str): + """ Remove blank lines from xml string """ + new_str = "" + for line in xml_str.split("\n"): + if not line.strip() == "\t" and not line.strip() == "": + new_str += line + "\n" + return new_str + + def run_simulation(self): + """ + Execute alaska simulation + """ + # define command line + command = self.alaska_exe+' '+os.path.join(self.run_directory, 'alaska_model', 'LoadCases', 'LC_UQ.xml') + # execute command + subprocess.run(command, cwd=self.run_directory) + + def extract_results(self): + """ + Extract the useful, desired signals which can be used for damping determination. + + Returns + ------- + result_dict : dict + Collected simulation results. Key names have to match with var_names specified in the postprocessor + settings in the config + """ + # get time series data + #c = CompressedIntLoader(self.alaska_int) + c = CompressedIntLoader(os.path.join(self.run_directory, 'alaska_model', 'Results', 'IWT-7.5-164_UQ_ref_12ms_ref_HS2.int')) + # get sensor names + names = [x.name for x in c.sensors()] + # hard coded radial sensor positions (geometry information only in blade description, not in int file) + radpos = np.array([2, 3.14290000000000, 4.20070000000000, 5.00120000000000, 6.10070000000000, 7.10080000000000, + 8.20070000000000, 8.90060000000000, 10.1002000000000, 10.9000000000000, 12.2000000000000, 13, + 13.9990000000000, 14.8990000000000, 15.9990000000000, 17.0520000000000, 17.9990000000000, + 18.9460000000000, 20.1040000000000, 21.1570000000000, 21.9990000000000, 23.1040000000000, + 23.9469000000000, 24.9989000000000, 25.8989000000000, 26.9989000000000, 27.9468000000000, + 29.0997000000000, 29.9996000000000, 30.9985000000000, 31.8984000000000, 32.9993000000000, + 33.9990000000000, 34.9988000000000, 36.0985000000000, 36.9983000000000, 37.9979000000000, + 38.9615000000000, 40.0320000000000, 40.9965000000000, 41.9958000000000, 42.9591000000000, + 44.0292000000000, 44.9934000000000, 45.9394000000000, 46.9912000000000, 47.9369000000000, + 48.9885000000000, 49.9339000000000, 50.9851000000000, 51.9831000000000, 52.9810000000000, + 53.9785000000000, 54.9761000000000, 55.9203000000000, 56.9701000000000, 57.9397000000000, + 58.9093000000000, 59.9858000000000, 61.0543000000000, 61.9499000000000, 62.9450000000000, + 63.9035000000000, 64.9694000000000, 65.9271000000000, 66.9205000000000, 67.9130000000000, + 68.8528000000000, 69.8967000000000, 70.9406000000000, 71.8782000000000, 72.8159000000000, + 73.9640000000000, 75.0327000000000, 75.8871000000000, 76.8474000000000, 77.7005000000000, + 78.7401000000000, 79.5692000000000, 80.5638000000000, 80.7605000000000, 81.2520000000000, + 81.7235900000000]) + + result_dict = dict() + result_dict['time'] = self.get_time(c) + result_dict['radpos'] = radpos + result_dict['deflectionInPlane_b1'] = -self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['UYR1' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['deflectionInPlane_b2'] = -self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['UYR2' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['deflectionInPlane_b3'] = -self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['UYR3' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['deflectionOutOfPlane_b1'] = self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['UZR1' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['deflectionOutOfPlane_b2'] = self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['UZR2' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['deflectionOutOfPlane_b3'] = self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['UZR3' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['torsion_b1'] = self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['RXV1' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['torsion_b2'] = self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['RXV2' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['torsion_b3'] = self.get_sensor_data_over_blade_position(c, [names.index(i) for i in + ['RXV3' + str(i) for i in + range(1, len(radpos))]], radpos) + result_dict['towertop_fa'] = self.get_sensor_data(c, names.index('GZ')) + result_dict['towertop_ss'] = self.get_sensor_data(c, names.index('GY')) + + return result_dict + + def get_time(self, c): + """ get time vector """ + time = np.arange(c.timeSerie(0).tstart, len(c.timeSerie(0).data)*c.timeSerie(0).tstep, c.timeSerie(0).tstep) + return time + + def get_sensor_data(self, c, idx): + sensordata = c.timeSerie(idx).data * c.timeSerie(idx).scale + return sensordata + + def get_sensor_data_over_blade_position(self, c, idx, radpos): + # get sensordata + sensordata = np.array([c.timeSerie(i).data for i in idx]) + sensordata = sensordata.astype('float64') # change dtype from int32 to float64 + scale = np.array([c.timeSerie(i).scale for i in idx]) + # scale sensordata + for i in range(len(scale)): + sensordata[i] = sensordata[i] * scale[i] + # extrapolate sensordata to last radial position + sensordata_extrap = np.zeros((len(sensordata) + 1, len(sensordata[0]))) + sensordata_extrap[0:-1, :] = sensordata + for i in range(len(sensordata[0])): + f = interpolate.interp1d(radpos[0:-1], sensordata[:, i], fill_value='extrapolate') + sensordata_extrap[-1, i] = f(radpos[-1]) # use interpolation function returned by 'interp1d' + return sensordata_extrap + + +# debug config +# 1) set run_type to "reference" (unchanged blade data) +# 2) set run_type to "test" (with blade data changes) +# +# execute: +# use "python run_analysis.py -c config/alaska/" + +if __name__ == "__main__": + + # raise NotImplementedError('No stand-alone running possible yet') + alaska_example = AlaskaModel() + alaska_example.create_simulation('') + alaska_example.run_simulation() + alaska_example.extract_results() \ No newline at end of file diff --git a/use_cases/stability_analysis/tool_interfaces/hawc2_interface.py b/use_cases/stability_analysis/tool_interfaces/hawc2_interface.py new file mode 100644 index 0000000..3853946 --- /dev/null +++ b/use_cases/stability_analysis/tool_interfaces/hawc2_interface.py @@ -0,0 +1,392 @@ +""" +@author: Otto Braun +@date: 05.01.2021 +""" + +import numpy as np +import os +from os import path +import glob +import shutil +from subprocess import call +import logging + +from tool_interfaces.simulation_model_interface import SimulationModel, SimulationError + + +class HAWC2Model(SimulationModel): + """ + HAWC2 model interface + + Parameters + ---------- + iteration_run_directory : string + Path to directory of this framework iteration + config : dict + User-specified settings for this tool + + Attributes + ---------- + config : dict + User-specified settings for this tool + iteration_run_directory : string + Path to directory of this framework iteration + """ + def __init__(self, iteration_run_directory, config): + self.iteration_run_directory = iteration_run_directory + self.config = config + + def create_simulation(self, blade_data): + """ + Setup modified HAWC2 model + + Parameters + ---------- + blade_data : dict + Dictionary with structural blade data + """ + logger = logging.getLogger('quexus.uq.model.run.hawc2model.create_simulation') + + self.htc_master_file_path = self.config['htc_master_file_path'] + self.master_file_path=self.config['master_file_path'] + + # copy the master model data + source_path = self.master_file_path+'\ModelData' + destination_path = self.iteration_run_directory+'\ModelData' + if path.exists(destination_path): + shutil.rmtree(destination_path) + shutil.copytree(source_path, destination_path) + + source_path = self.master_file_path+'\htc_Snippets' + destination_path = self.iteration_run_directory+'\htc_Snippets' + if path.exists(destination_path): + shutil.rmtree(destination_path) + shutil.copytree(source_path, destination_path) + + source_path = self.master_file_path+'\OutputSensors' + destination_path = self.iteration_run_directory+'\OutputSensors' + if path.exists(destination_path): + shutil.rmtree(destination_path) + shutil.copytree(source_path, destination_path) + + source_path = self.master_file_path+'\Wind' + destination_path = self.iteration_run_directory+'\Wind' + if path.exists(destination_path): + shutil.rmtree(destination_path) + shutil.copytree(source_path, destination_path) + + source_path = self.master_file_path+'\Control' + destination_path = self.iteration_run_directory+'\Control' + if path.exists(destination_path): + shutil.rmtree(destination_path) + shutil.copytree(source_path, destination_path) + + # copy batchfiles from master model + for file in glob.glob(self.master_file_path+'/*.bat'): + shutil.copy(file, self.iteration_run_directory) + + with open(self.htc_master_file_path, 'r') as htc_master_file: + htc_content = htc_master_file.readlines() + self.write_st(blade_data) + self.modify_htc(htc_content) + logger.info('Hawc2 Simulation created sucessfully') + + def write_st(self, blade_data): + """ + Write blade_data to HAWC2 .st file + + Parameters + ---------- + blade_data : dict + Dictionary with structural blade data + """ + logger = logging.getLogger('quexus.uq.model.run.hawc2model.create_simulation.write_st') + + if self.config['fpm_bool'] == '0': + + st_file= np.c_[blade_data['h2_beam']['arc'], + blade_data['h2_beam']['m_pm'], + blade_data['h2_beam']['x_cg'], + blade_data['h2_beam']['y_cg'], + blade_data['h2_beam']['ri_x'], + blade_data['h2_beam']['ri_y'], + blade_data['h2_beam']['x_sc'], + blade_data['h2_beam']['y_sc'], + blade_data['h2_beam']['E'], + blade_data['h2_beam']['G'], + blade_data['h2_beam']['I_x'], + blade_data['h2_beam']['I_y'], + blade_data['h2_beam']['I_T'], + blade_data['h2_beam']['k_x'], + blade_data['h2_beam']['k_y'], + blade_data['h2_beam']['A'], + blade_data['h2_beam']['theta_s'], + blade_data['h2_beam']['x_ec'], + blade_data['h2_beam']['y_ec']] + + header = open(os.path.join(os.getcwd(), 'tool_interfaces', 'header_st_file_hawc2.txt'), 'r') + header_content = header.read() + header.close() + np.savetxt(os.path.join(self.iteration_run_directory, 'st_file.inp'), st_file, newline="\n ", + delimiter=" ", header=header_content, comments='', fmt='%1.12e') + + elif self.config['fpm_bool'] == '1': + + st_file= np.c_[blade_data['h2_FPM']['arc'], + blade_data['h2_FPM']['m_pm'], + blade_data['h2_FPM']['x_cg'], + blade_data['h2_FPM']['y_cg'], + blade_data['h2_FPM']['ri_x'], + blade_data['h2_FPM']['ri_y'], + blade_data['h2_FPM']['theta_s'], + blade_data['h2_FPM']['x_ec'], + blade_data['h2_FPM']['y_ec'], + blade_data['h2_FPM']['K_11'], + blade_data['h2_FPM']['K_12'], + blade_data['h2_FPM']['K_13'], + blade_data['h2_FPM']['K_14'], + blade_data['h2_FPM']['K_15'], + blade_data['h2_FPM']['K_16'], + blade_data['h2_FPM']['K_22'], + blade_data['h2_FPM']['K_23'], + blade_data['h2_FPM']['K_24'], + blade_data['h2_FPM']['K_25'], + blade_data['h2_FPM']['K_26'], + blade_data['h2_FPM']['K_33'], + blade_data['h2_FPM']['K_34'], + blade_data['h2_FPM']['K_35'], + blade_data['h2_FPM']['K_36'], + blade_data['h2_FPM']['K_44'], + blade_data['h2_FPM']['K_45'], + blade_data['h2_FPM']['K_46'], + blade_data['h2_FPM']['K_55'], + blade_data['h2_FPM']['K_56'], + blade_data['h2_FPM']['K_66']] + + header_6x6 = open(os.path.join(os.getcwd(),'tool_interfaces','header_st_file_hawc2_6x6.txt'), 'r') + header_6x6_content = header_6x6.read() + header_6x6.close() + np.savetxt(os.path.join(self.iteration_run_directory, 'st_file.inp'), st_file, newline="\n ", + delimiter=" ", header=header_6x6_content, comments='', fmt='%1.12e') + + else: + logger.exception('fpm_bool not assigned correctly, options are 0 and 1') + raise ValueError('fpm_bool not assigned correctly, options are 0 and 1') + + def modify_htc(self, htc_content): + """ Modify htc file """ + keywords = self.config["keywords"] + for i, key in enumerate(keywords): + key = key.replace('absoulte_path+', self.config["master_file_path"]) + key = key.replace('\\', '/') + keywords[i] = key + + replacements = self.config["replacements"] + if self.config['fpm_bool'] == '1': + for i, key in enumerate(self.config["fpm_keywords"]): + keywords.append(self.config["fpm_keywords"][i]) + replacements.append(self.config["fpm_replacements"][i]) + + for i, key in enumerate(replacements): + key = key.replace('absoulte_path+', self.config["master_file_path"]) + key = key.replace('\\', '/') + if key == './st_file.inp;' and self.config['fpm_bool'] == '1': + key = './st_file.inp;' + key = key+'\n\t\t\tFPM 1 ;' + replacements[i] = key + + with open(os.path.join(self.iteration_run_directory, "HAWC2.htc"), 'w') as output: + for idx, lines in enumerate(htc_content): + for i in range(len(keywords)): + if keywords[i] in lines: + lines = lines.replace(keywords[i], replacements[i]) + htc_content[idx] = lines + output.write(str(lines)) + + def run_simulation(self): + """ + Execute HAWC2 simulation + """ + exe_homedir = self.config['exe_path'] + exe_name = self.config['exe_name'] + call(exe_homedir+exe_name+' HAWC2.htc', cwd=self.iteration_run_directory) + + def extract_results(self): + """ + Extract the useful, desired signals which can be used for damping determination. + + Returns + ------- + result_dict : dict + Collected simulation results. Key names have to match with var_names specified in the postprocessor + settings in the config + + Raises + ------ + FileNotFoundError + If a .sel file could not be found in the iteration_run_directory + + Notes + ----- + Takes the .sel file in the iteration_run_directory to create sensor_key, + these keys are use to look up time series in the .dat-file in the iteration run_directory based on indices + extracted from the .sel file + + resulting dictionary result dict: + result_dict[sensor_key]=# time steps + """ + + logger = logging.getLogger('quexus.uq.model.run.hawc2model.extract_results') + try: + idents = [] + sensor_key_list_file = glob.glob(self.iteration_run_directory + "/**/*.sel", recursive=True) + + with open(sensor_key_list_file[0]) as hawc_keys: + hawc_keys_content = hawc_keys.readlines() + hawc_keys_content = hawc_keys_content[12:-1] + for line in hawc_keys_content: + ident = line.split(maxsplit=1)[1].replace('\t', ' ') + while ' ' in ident: + ident = ident.replace(' ', ' ') + idents.append(ident) + + time_key = 'Time s Time\n' + tower_ss_key = "Dx coo: tower m elastic Deflection tower Dx Mbdy:tower s= 116.00[m] s/S= 1.00 coo: tower center:default\n" + tower_fa_key = "Dy coo: tower m elastic Deflection tower Dy Mbdy:tower s= 116.00[m] s/S= 1.00 coo: tower center:default\n" + + torsion_keys = [] + InPlane_keys = [] + OutOfPlane_keys = [] + for i in[1, 2, 3]: + torsion_keys.append('Rz coo: blade' + str(i) + ' deg elastic Rotation blade' + str(i) + ' Rz Mbdy:blade'+ str(i)) + InPlane_keys.append('Dx coo: blade' + str(i) + ' m elastic Deflection blade' + str(i) + ' Dx Mbdy:blade' + str(i)) + OutOfPlane_keys.append('Dy coo: blade' + str(i) + ' m elastic Deflection blade' + str(i) + ' Dy Mbdy:blade' + str(i)) + + torsion_all_key = dict() + InPlane_all_key = dict() + OutOfPlane_all_key = dict() + + for key in torsion_keys: + torsion_all_key[key] = [] + + for key in InPlane_keys: + InPlane_all_key[key] = [] + + for key in OutOfPlane_keys: + OutOfPlane_all_key[key] = [] + + for lines in idents: + for key in torsion_keys: + if key in lines: + torsion_all_key[key].append(lines) + for key in InPlane_keys: + if key in lines: + InPlane_all_key[key].append(lines) + for key in OutOfPlane_all_key: + if key in lines: + OutOfPlane_all_key[key].append(lines) + + radpos = [] + for line in OutOfPlane_all_key[key]: + ident = line.split('s=', maxsplit=1)[1] + ident = float(ident.split('[m]', maxsplit=1)[0]) + radpos.append(ident) + + sensor_data_file = glob.glob(self.iteration_run_directory + "/**/*.dat", recursive=True) + if len(sensor_data_file) > 1: + raise ValueError("More than one .dat file was found in " + self.iteration_run_directory) + if len(sensor_data_file) < 1: + raise ValueError("No .dat file was found in " + self.iteration_run_directory) + + hawc2data = np.loadtxt(sensor_data_file[0]) + + result_dict = dict() + result_dict['time'] = hawc2data[:, idents.index(time_key)] + result_dict['radpos'] = np.transpose(np.asarray(radpos)) + result_dict['towertop_fa'] = hawc2data[:, idents.index(tower_fa_key)] + result_dict['towertop_ss'] = hawc2data[:, idents.index(tower_ss_key)] + + for i, key in enumerate(torsion_all_key): + result_dict['torsion_b'+str(i+1)] = [] + for single_key in torsion_all_key[key]: + result_dict['torsion_b'+str(i+1)].append(hawc2data[:, idents.index(single_key)]) + result_dict['torsion_b'+str(i+1)] = np.transpose(np.asarray(result_dict['torsion_b'+str(i+1)])) + + for i, key in enumerate(InPlane_all_key): + result_dict['deflectionInPlane_b'+str(i+1)] = [] + for single_key in InPlane_all_key[key]: + result_dict['deflectionInPlane_b'+str(i+1)].append(hawc2data[:, idents.index(single_key)]) + result_dict['deflectionInPlane_b'+str(i+1)] = np.transpose(np.asarray(result_dict['deflectionInPlane_b'+str(i+1)])) + + for i, key in enumerate(OutOfPlane_all_key): + result_dict['deflectionOutOfPlane_b'+str(i+1)] = [] + for single_key in OutOfPlane_all_key[key]: + result_dict['deflectionOutOfPlane_b'+str(i+1)].append(hawc2data[:, idents.index(single_key)]) + result_dict['deflectionOutOfPlane_b'+str(i+1)] = np.transpose(np.asarray(result_dict['deflectionOutOfPlane_b'+str(i+1)])) + + return result_dict + + except FileNotFoundError(): + logger.warning('No .sel file found in ' + self.iteration_run_directory) + return None + + def parse_h2_st(self, h2_st_file): + """ + Read a HAWC2 .st file + + Parameters + ---------- + h2_st_file : string + Path to HAWC2 structural blade data file + + Returns + ------- + h2_st : dict + HAWC2 blade data + + Notes + ----- + HAWC2 structural parameters: + arc + m_pm + x_cg + y_cg + ri_x + ri_y + x_sc + y_sc + E + G + I_x + I_y + I_T + k_x + k_y + A + theta_s + x_ec + y_ec + """ + h2_st = dict() + table_struct_h2 = np.loadtxt(h2_st_file, skiprows=8) + h2_st['arc'] = table_struct_h2[:, 0] + h2_st['m_pm'] = table_struct_h2[:, 1] + h2_st['x_cg'] = table_struct_h2[:, 2] + h2_st['y_cg'] = table_struct_h2[:, 3] + h2_st['ri_x'] = table_struct_h2[:, 4] + h2_st['ri_y'] = table_struct_h2[:, 5] + h2_st['x_sc'] = table_struct_h2[:, 6] + h2_st['y_sc'] = table_struct_h2[:, 7] + h2_st['E'] = table_struct_h2[:, 8] + h2_st['G'] = table_struct_h2[:, 9] + h2_st['I_x'] = table_struct_h2[:, 10] + h2_st['I_y'] = table_struct_h2[:, 11] + h2_st['I_T'] = table_struct_h2[:, 12] + h2_st['k_x'] = table_struct_h2[:, 13] + h2_st['k_y'] = table_struct_h2[:, 14] + h2_st['A'] = table_struct_h2[:, 15] + h2_st['theta_s'] = table_struct_h2[:, 16] + h2_st['x_ec'] = table_struct_h2[:, 17] + h2_st['y_ec'] = table_struct_h2[:, 18] + + return h2_st \ No newline at end of file diff --git a/use_cases/stability_analysis/tool_interfaces/hawcstab2_interface.py b/use_cases/stability_analysis/tool_interfaces/hawcstab2_interface.py new file mode 100644 index 0000000..0bff2ad --- /dev/null +++ b/use_cases/stability_analysis/tool_interfaces/hawcstab2_interface.py @@ -0,0 +1,68 @@ +""" +@author: Otto Braun +@date: 05.01.2021 +""" + +import os + +from tool_interfaces.HS2ASCIIReader import HS2Data +from tool_interfaces.simulation_model_interface import SimulationModel, SimulationError +from tool_interfaces.hawc2_interface import HAWC2Model + + +class HAWCStab2Model(HAWC2Model): + """ + HAWCStab2 model interface + + Parameters + ---------- + iteration_run_directory : string + Path to directory of this framework iteration + config : dict + User-specified settings for this tool + + Attributes + ---------- + config : dict + User-specified settings for this tool + iteration_run_directory : string + Path to directory of this framework iteration + """ + + def extract_results(self): + """ + Extract the useful, desired signals which can be used for damping determination. + + Returns + ------- + result_dict : dict + Collected simulation results. Key names have to match with var_names specified in the postprocessor + settings in the config + """ + try: + self.cmb = os.path.join(self.iteration_run_directory, "HAWC2.cmb") + self.amp= os.path.join(self.iteration_run_directory, "HAWC2.amp") + self.opt= os.path.join(self.master_file_path, "Control", "IWT_7.5-164_Rev-2.5.2_HS2_coarse_stability.opt") + HS2_results = HS2Data(self.cmb, self.amp, self.opt) + HS2_results.read_data() + + # extracting damping of 'critical' mode + # for a HS2 run with multiple operating points, the HS2_results.modes.damping matrix has the shape: + # # of modes x # of wind speeds x 2 + + # find index for windspeed == 12 m/s + desired_wind_speed = 12 # m/s + idx_12ms = HS2_results.modes.damping[0, :, 0].tolist().index(desired_wind_speed) + + # find lowest damping for this wind speed + # NOTE! there are some -100 entries, so only damping values up to -10 % are assumed to be feasible + all_damping_values = HS2_results.modes.damping[:, idx_12ms, 1] + accepted_damping_values = all_damping_values[(all_damping_values > -10)] + + result_dict = dict() + result_dict['damping_criticalmode'] = min(accepted_damping_values) + print(result_dict['damping_criticalmode']) + return result_dict + + except: + print('no HS2 result data were found') \ No newline at end of file diff --git a/use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2.txt b/use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2.txt new file mode 100644 index 0000000..cf05759 --- /dev/null +++ b/use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2.txt @@ -0,0 +1,4 @@ +#1 Main data set for the blade + Radius m x_m y_m r_ix r_iy x_s y_s E G I_x I_y K k_x k_y A theta_s x_e y_e + [m] [kg/m] [m] [m] [m] [m] [m] [m] [N/m^2] [N/m^2] [m^4] [m^4] [m^4/rad] [-] [-] [m^2] [deg] [m] [m] +$1 83 sections flexible blade properties modified by spline method diff --git a/use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2_6x6.txt b/use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2_6x6.txt new file mode 100644 index 0000000..dcdaf98 --- /dev/null +++ b/use_cases/stability_analysis/tool_interfaces/header_st_file_hawc2_6x6.txt @@ -0,0 +1,7 @@ +1 number of sets, Nset +----------------- +#1 BECAS-Result-Data for HAWC2/HAWCStab2 +================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================ + arc [01] m [02] x_cg [03] y_cg [04] ri_x [05] ri_y [06] theta [07] x_e [08] y_e [09] K_11 [10] K_12 [11] K_13 [12] K_14 [13] K_15 [14] K_16 [15] K_22 [16] K_23 [17] K_24 [18] K_25 [19] K_26 [20] K_33 [21] K_34 [22] K_35 [23] K_36 [24] K_44 [25] K_45 [26] K_46 [27] K_55 [28] K_56 [29] K_66 [30] +================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================ +$1 83 Flexible \ No newline at end of file