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():
213def main():
214    wind_speed = 4
215    simulation_time = 20
216    f_calculateDamageAndFatigueLife(wind_speed, simulation_time)