sensotwin.lifetime_extrapolation.simulation.lib.fatigueAnalysisCCX
1import pandas as pd 2import numpy as np 3import os 4from pathlib import Path 5from . import readVTT 6import vtk 7from vtk.util import numpy_support as vtknp 8import sys 9 10def f_calculateDamageAndFatigueLife(wind_speed, simulation_time, working_dir = './'): 11 ''' 12 Calculates the damage and the fatigue life of every element 13 based on constant uniform wind, shifted linear Goodman diagrams 14 and Miner's sum. 15 16 This follows the calculation method by Germanischer Lloyd for 17 fibre-reinforced polymers under fatigue loads, p. 5-22 f.: 18 Germanischer Lloyd: Rules and Regulations, IV-Non-marine 19 Technology, Part 1-Wind Energy In: Guidelines for the 20 Certification of Wind Turbines. Guidelines for the 21 Certification of Wind Turbines, 2010. 22 23 Parameters: 24 - wind_speed (float): 25 Average wind speed acting on the entire rotor disc in m/s. 26 - simulation_time (float): 27 Total time the rotor blade is exposed to wind_speed in years. 28 - working_dir (string): 29 Working directory to execute file in. 30 31 Returns: 32 - None 33 ''' 34 35 # read in material data from table 36 current_dir = working_dir 37 material_data_df = readVTT.f_readVTT2DF(os.path.join(current_dir, 'IN', 'material_fatigue_data.vtt')) 38 39 # Compute gamma_Ma and gamma_Mb 40 material_data_df["gamma_Ma"] = material_data_df.gamma_M0*material_data_df.C1a*material_data_df.C2a*material_data_df.C3a*material_data_df.C4a; 41 material_data_df["gamma_Mb"] = material_data_df.gamma_M0*material_data_df.C2b*material_data_df.C3b*material_data_df.C4b*material_data_df.C5b; 42 43 # define the number of azimuthal positions 44 AzimuthalPositions = 4; 45 StepSize = 360/AzimuthalPositions; 46 47 # import RPM Data 48 rpm_df = readVTT.f_readVTT2DF(os.path.join(current_dir,'IN', 'rpm_windspeed_0p01.vtt')) 49 50 vwind_values = rpm_df["windspeed_mpers"] 51 rpm_values = rpm_df["rpm"] 52 53 rpm_at_wind_speed_interp = np.interp(wind_speed,vwind_values,rpm_values) 54 total_rotations = int(round((simulation_time*365*24*60)*rpm_at_wind_speed_interp,0)); 55 56 # import stress data 57 58 # for aerodynamics and centrifugal forces 59 file_path = os.path.join(Path(current_dir).parent, 'Static', 'OUT', 'AerodynamicsCentrifugal') 60 all_files = os.listdir(file_path) 61 all_processed_vtu_files = [x for x in all_files if x.endswith('_processed.vtu')] 62 filename = all_processed_vtu_files[0] 63 64 reader = vtk.vtkXMLUnstructuredGridReader() 65 reader.SetFileName(os.path.join(file_path,filename)) 66 reader.Update() 67 input_data_aerocent = reader.GetOutput() 68 69 70 # for gravity forces 71 all_files = os.listdir(os.path.join(Path(current_dir).parent, 'Static', 'OUT')) 72 all_azimuth_folders = [x for x in all_files if x.startswith('Gravity')] 73 74 input_data_gravity = [] 75 gravity_stresses = [] 76 77 for i,folder in enumerate(all_azimuth_folders): 78 file_path = os.path.join(Path(current_dir).parent, 'Static', 'OUT', folder) 79 all_files = os.listdir(file_path) 80 all_processed_vtu_files = [x for x in all_files if x.endswith('_processed.vtu')] 81 filename = all_processed_vtu_files[0] 82 83 reader = vtk.vtkXMLUnstructuredGridReader() 84 reader.SetFileName(os.path.join(file_path,filename)) 85 reader.Update() 86 input_data_gravity.append(reader.GetOutput()) 87 88 gravity_stresses.append([x[0] for x in vtknp.vtk_to_numpy(input_data_gravity[i].GetCellData().GetArray('S'))]) # get first component of stress vector (S11) 89 90 # get LAYIDS for element assignment 91 material_names = list(material_data_df["Material"]) 92 layer_ids = [] 93 94 for material in material_names: 95 current_array = input_data_aerocent.GetCellData().GetArray('LAYID_' + material) 96 if current_array is None: 97 print('WARNING: Material ' + material + ' not found in results. Ignoring it.') 98 else: 99 layer_ids.append([material, vtknp.vtk_to_numpy(current_array)]) 100 101 # compute fatigue for each material 102 miner_sum_materials = [] 103 fatigue_life_materials = [] 104 105 for i,material in enumerate(material_data_df["Material"]): 106 m = float(material_data_df.iloc[i]['m']) 107 UTS = float(material_data_df.iloc[i]['UTS']) 108 UCS = float(material_data_df.iloc[i]['UCS']) 109 gamma_Ma = float(material_data_df.iloc[i]['gamma_Ma']) 110 gamma_Mb = float(material_data_df.iloc[i]['gamma_Mb']) 111 112 # get stresses for that material 113 # find all stresses where the field LAYIDS_<current_material> is not NaN 114 current_aerocent_stresses = [] 115 for j,x in enumerate(vtknp.vtk_to_numpy(input_data_aerocent.GetCellData().GetArray('S'))): 116 if not np.isnan(layer_ids[i][1][j]): 117 current_aerocent_stresses.append(x[0]) 118 else: 119 current_aerocent_stresses.append(np.nan) 120 121 current_total_stresses = [] 122 123 # for every azimuth angle 124 for input_data in input_data_gravity: 125 current_gravity_stresses_azimuth = [] 126 for j,x in enumerate(vtknp.vtk_to_numpy(input_data.GetCellData().GetArray('S'))): 127 if not np.isnan(layer_ids[i][1][j]): 128 current_gravity_stresses_azimuth.append(x[0]) 129 else: 130 current_gravity_stresses_azimuth.append(np.nan) 131 132 current_total_stresses_azimuth = [] 133 134 for j in range(len(current_gravity_stresses_azimuth)): # for every element 135 # compute total stress as superposition of aerodynamical, centrifugal and gravity stresses for every azimuth angle 136 current_total_stresses_azimuth.append(current_aerocent_stresses[j] + current_gravity_stresses_azimuth[j]) 137 138 current_total_stresses.append(current_total_stresses_azimuth) # list of lists [[value_per_element, value_per_element, ...], <-- azimuth_angle_1 139 # [value_per_element, value_per_element, ...], ...] <-- azimuth_angle_2 140 141 current_miner_sum = [] 142 current_fatigue_life = [] 143 144 for j in range(len(current_total_stresses[0])): # for every element 145 sigma_max = max([x[j] for x in current_total_stresses]); 146 sigma_min = min([x[j] for x in current_total_stresses]); 147 sigma_mean = (sigma_max+sigma_min)/2; 148 sigma_amplitude = sigma_max-sigma_mean; 149 150 # compute allowable number and Miner's sum 151 N_allowable = np.power(((UTS + abs(UCS) - abs(2*gamma_Ma*sigma_mean - UTS + abs(UCS)))/(2*gamma_Mb*sigma_amplitude)),m) 152 miner_sum = total_rotations/N_allowable 153 154 current_miner_sum.append(miner_sum) 155 current_fatigue_life.append(simulation_time/miner_sum) 156 157 miner_sum_materials.append(current_miner_sum) 158 fatigue_life_materials.append(current_fatigue_life) 159 160 # write results 161 file_path = os.path.join(Path(current_dir).parent, 'Static', 'OUT', 'AerodynamicsCentrifugal') 162 all_files = os.listdir(file_path) 163 all_processed_vtu_files = [x for x in all_files if x.endswith('_processed.vtu')] 164 filename = all_processed_vtu_files[0] 165 166 reader = vtk.vtkXMLUnstructuredGridReader() 167 reader.SetFileName(os.path.join(file_path,filename)) 168 reader.Update() 169 input_data = reader.GetOutput() 170 171 # delete S, E, U and ERROR cell data array 172 dataset = input_data.GetCellData() 173 dataset.RemoveArray('S') 174 dataset.RemoveArray('S_Mises') 175 dataset.RemoveArray('S_Principal') 176 dataset.RemoveArray('E') 177 dataset.RemoveArray('E_Mises') 178 dataset.RemoveArray('E_Principal') 179 dataset.RemoveArray('U') 180 dataset.RemoveArray('ERROR') 181 dataset.Modified() 182 183 # rename U point data to U_AEROCENT as it only shows the deformation due to aerodynamical and centrifugal forces 184 input_data.GetPointData().GetArray('U').SetName('U_AEROCENT') 185 dataset.Modified() 186 187 # create new array with damage and fatigue life for each material 188 189 for i,material in enumerate(material_data_df["Material"]): 190 array = vtknp.numpy_to_vtk(miner_sum_materials[i]) 191 array.SetName('DFAT_' + material) 192 dataset = input_data.GetCellData() 193 dataset.AddArray(array) 194 dataset.Modified() 195 196 array = vtknp.numpy_to_vtk(fatigue_life_materials[i]) 197 array.SetName('TFAT_' + material) 198 dataset = input_data.GetCellData() 199 dataset.AddArray(array) 200 dataset.Modified() 201 202 # write to new .vtu file 203 if not os.path.exists(os.path.join(current_dir, 'OUT')): 204 os.mkdir(os.path.join(current_dir, 'OUT')) 205 206 writer = vtk.vtkXMLUnstructuredGridWriter() 207 writer.SetFileName(os.path.join(current_dir, 'OUT', filename.strip('_processed.vtu') + '_FATIGUE.vtu')) 208 writer.CompressorType = 'NONE' 209 writer.SetInputData(input_data) 210 writer.Write() 211 212def main(): 213 wind_speed = 4 214 simulation_time = 20 215 f_calculateDamageAndFatigueLife(wind_speed, simulation_time) 216 217if __name__ == '__main__': 218 main()
def
f_calculateDamageAndFatigueLife(wind_speed, simulation_time, working_dir='./'):
11def f_calculateDamageAndFatigueLife(wind_speed, simulation_time, working_dir = './'): 12 ''' 13 Calculates the damage and the fatigue life of every element 14 based on constant uniform wind, shifted linear Goodman diagrams 15 and Miner's sum. 16 17 This follows the calculation method by Germanischer Lloyd for 18 fibre-reinforced polymers under fatigue loads, p. 5-22 f.: 19 Germanischer Lloyd: Rules and Regulations, IV-Non-marine 20 Technology, Part 1-Wind Energy In: Guidelines for the 21 Certification of Wind Turbines. Guidelines for the 22 Certification of Wind Turbines, 2010. 23 24 Parameters: 25 - wind_speed (float): 26 Average wind speed acting on the entire rotor disc in m/s. 27 - simulation_time (float): 28 Total time the rotor blade is exposed to wind_speed in years. 29 - working_dir (string): 30 Working directory to execute file in. 31 32 Returns: 33 - None 34 ''' 35 36 # read in material data from table 37 current_dir = working_dir 38 material_data_df = readVTT.f_readVTT2DF(os.path.join(current_dir, 'IN', 'material_fatigue_data.vtt')) 39 40 # Compute gamma_Ma and gamma_Mb 41 material_data_df["gamma_Ma"] = material_data_df.gamma_M0*material_data_df.C1a*material_data_df.C2a*material_data_df.C3a*material_data_df.C4a; 42 material_data_df["gamma_Mb"] = material_data_df.gamma_M0*material_data_df.C2b*material_data_df.C3b*material_data_df.C4b*material_data_df.C5b; 43 44 # define the number of azimuthal positions 45 AzimuthalPositions = 4; 46 StepSize = 360/AzimuthalPositions; 47 48 # import RPM Data 49 rpm_df = readVTT.f_readVTT2DF(os.path.join(current_dir,'IN', 'rpm_windspeed_0p01.vtt')) 50 51 vwind_values = rpm_df["windspeed_mpers"] 52 rpm_values = rpm_df["rpm"] 53 54 rpm_at_wind_speed_interp = np.interp(wind_speed,vwind_values,rpm_values) 55 total_rotations = int(round((simulation_time*365*24*60)*rpm_at_wind_speed_interp,0)); 56 57 # import stress data 58 59 # for aerodynamics and centrifugal forces 60 file_path = os.path.join(Path(current_dir).parent, 'Static', 'OUT', 'AerodynamicsCentrifugal') 61 all_files = os.listdir(file_path) 62 all_processed_vtu_files = [x for x in all_files if x.endswith('_processed.vtu')] 63 filename = all_processed_vtu_files[0] 64 65 reader = vtk.vtkXMLUnstructuredGridReader() 66 reader.SetFileName(os.path.join(file_path,filename)) 67 reader.Update() 68 input_data_aerocent = reader.GetOutput() 69 70 71 # for gravity forces 72 all_files = os.listdir(os.path.join(Path(current_dir).parent, 'Static', 'OUT')) 73 all_azimuth_folders = [x for x in all_files if x.startswith('Gravity')] 74 75 input_data_gravity = [] 76 gravity_stresses = [] 77 78 for i,folder in enumerate(all_azimuth_folders): 79 file_path = os.path.join(Path(current_dir).parent, 'Static', 'OUT', folder) 80 all_files = os.listdir(file_path) 81 all_processed_vtu_files = [x for x in all_files if x.endswith('_processed.vtu')] 82 filename = all_processed_vtu_files[0] 83 84 reader = vtk.vtkXMLUnstructuredGridReader() 85 reader.SetFileName(os.path.join(file_path,filename)) 86 reader.Update() 87 input_data_gravity.append(reader.GetOutput()) 88 89 gravity_stresses.append([x[0] for x in vtknp.vtk_to_numpy(input_data_gravity[i].GetCellData().GetArray('S'))]) # get first component of stress vector (S11) 90 91 # get LAYIDS for element assignment 92 material_names = list(material_data_df["Material"]) 93 layer_ids = [] 94 95 for material in material_names: 96 current_array = input_data_aerocent.GetCellData().GetArray('LAYID_' + material) 97 if current_array is None: 98 print('WARNING: Material ' + material + ' not found in results. Ignoring it.') 99 else: 100 layer_ids.append([material, vtknp.vtk_to_numpy(current_array)]) 101 102 # compute fatigue for each material 103 miner_sum_materials = [] 104 fatigue_life_materials = [] 105 106 for i,material in enumerate(material_data_df["Material"]): 107 m = float(material_data_df.iloc[i]['m']) 108 UTS = float(material_data_df.iloc[i]['UTS']) 109 UCS = float(material_data_df.iloc[i]['UCS']) 110 gamma_Ma = float(material_data_df.iloc[i]['gamma_Ma']) 111 gamma_Mb = float(material_data_df.iloc[i]['gamma_Mb']) 112 113 # get stresses for that material 114 # find all stresses where the field LAYIDS_<current_material> is not NaN 115 current_aerocent_stresses = [] 116 for j,x in enumerate(vtknp.vtk_to_numpy(input_data_aerocent.GetCellData().GetArray('S'))): 117 if not np.isnan(layer_ids[i][1][j]): 118 current_aerocent_stresses.append(x[0]) 119 else: 120 current_aerocent_stresses.append(np.nan) 121 122 current_total_stresses = [] 123 124 # for every azimuth angle 125 for input_data in input_data_gravity: 126 current_gravity_stresses_azimuth = [] 127 for j,x in enumerate(vtknp.vtk_to_numpy(input_data.GetCellData().GetArray('S'))): 128 if not np.isnan(layer_ids[i][1][j]): 129 current_gravity_stresses_azimuth.append(x[0]) 130 else: 131 current_gravity_stresses_azimuth.append(np.nan) 132 133 current_total_stresses_azimuth = [] 134 135 for j in range(len(current_gravity_stresses_azimuth)): # for every element 136 # compute total stress as superposition of aerodynamical, centrifugal and gravity stresses for every azimuth angle 137 current_total_stresses_azimuth.append(current_aerocent_stresses[j] + current_gravity_stresses_azimuth[j]) 138 139 current_total_stresses.append(current_total_stresses_azimuth) # list of lists [[value_per_element, value_per_element, ...], <-- azimuth_angle_1 140 # [value_per_element, value_per_element, ...], ...] <-- azimuth_angle_2 141 142 current_miner_sum = [] 143 current_fatigue_life = [] 144 145 for j in range(len(current_total_stresses[0])): # for every element 146 sigma_max = max([x[j] for x in current_total_stresses]); 147 sigma_min = min([x[j] for x in current_total_stresses]); 148 sigma_mean = (sigma_max+sigma_min)/2; 149 sigma_amplitude = sigma_max-sigma_mean; 150 151 # compute allowable number and Miner's sum 152 N_allowable = np.power(((UTS + abs(UCS) - abs(2*gamma_Ma*sigma_mean - UTS + abs(UCS)))/(2*gamma_Mb*sigma_amplitude)),m) 153 miner_sum = total_rotations/N_allowable 154 155 current_miner_sum.append(miner_sum) 156 current_fatigue_life.append(simulation_time/miner_sum) 157 158 miner_sum_materials.append(current_miner_sum) 159 fatigue_life_materials.append(current_fatigue_life) 160 161 # write results 162 file_path = os.path.join(Path(current_dir).parent, 'Static', 'OUT', 'AerodynamicsCentrifugal') 163 all_files = os.listdir(file_path) 164 all_processed_vtu_files = [x for x in all_files if x.endswith('_processed.vtu')] 165 filename = all_processed_vtu_files[0] 166 167 reader = vtk.vtkXMLUnstructuredGridReader() 168 reader.SetFileName(os.path.join(file_path,filename)) 169 reader.Update() 170 input_data = reader.GetOutput() 171 172 # delete S, E, U and ERROR cell data array 173 dataset = input_data.GetCellData() 174 dataset.RemoveArray('S') 175 dataset.RemoveArray('S_Mises') 176 dataset.RemoveArray('S_Principal') 177 dataset.RemoveArray('E') 178 dataset.RemoveArray('E_Mises') 179 dataset.RemoveArray('E_Principal') 180 dataset.RemoveArray('U') 181 dataset.RemoveArray('ERROR') 182 dataset.Modified() 183 184 # rename U point data to U_AEROCENT as it only shows the deformation due to aerodynamical and centrifugal forces 185 input_data.GetPointData().GetArray('U').SetName('U_AEROCENT') 186 dataset.Modified() 187 188 # create new array with damage and fatigue life for each material 189 190 for i,material in enumerate(material_data_df["Material"]): 191 array = vtknp.numpy_to_vtk(miner_sum_materials[i]) 192 array.SetName('DFAT_' + material) 193 dataset = input_data.GetCellData() 194 dataset.AddArray(array) 195 dataset.Modified() 196 197 array = vtknp.numpy_to_vtk(fatigue_life_materials[i]) 198 array.SetName('TFAT_' + material) 199 dataset = input_data.GetCellData() 200 dataset.AddArray(array) 201 dataset.Modified() 202 203 # write to new .vtu file 204 if not os.path.exists(os.path.join(current_dir, 'OUT')): 205 os.mkdir(os.path.join(current_dir, 'OUT')) 206 207 writer = vtk.vtkXMLUnstructuredGridWriter() 208 writer.SetFileName(os.path.join(current_dir, 'OUT', filename.strip('_processed.vtu') + '_FATIGUE.vtu')) 209 writer.CompressorType = 'NONE' 210 writer.SetInputData(input_data) 211 writer.Write()
Calculates the damage and the fatigue life of every element based on constant uniform wind, shifted linear Goodman diagrams and Miner's sum.
This follows the calculation method by Germanischer Lloyd for fibre-reinforced polymers under fatigue loads, p. 5-22 f.: Germanischer Lloyd: Rules and Regulations, IV-Non-marine Technology, Part 1-Wind Energy In: Guidelines for the Certification of Wind Turbines. Guidelines for the Certification of Wind Turbines, 2010.
Parameters:
- wind_speed (float): Average wind speed acting on the entire rotor disc in m/s.
- simulation_time (float): Total time the rotor blade is exposed to wind_speed in years.
- working_dir (string): Working directory to execute file in.
Returns:
- None
def
main():