Newer
Older
# bea1 angle_speed rpm shaft_nacelle angle speed
elif self.ch_details[ch, 0].startswith('bea'):
output_type = self.ch_details[ch, 0].split(' ')[1]
# there is no label option for the bearing output
# and tag it
tag = 'bearing-%s-%s-%s' % (bearing_name, output_type, units)
# save all info in the dict
channelinfo = {}
channelinfo['bearing_name'] = bearing_name
channelinfo['output_type'] = output_type
channelinfo['units'] = units
channelinfo['chi'] = ch
# -----------------------------------------------------------------
# AS DEFINED IN: ch_aero
# AERO CL, CD, CM, VREL, ALFA, LIFT, DRAG, etc
# Cl, R= 0.5 deg Cl of blade 1 at radius 0.49
# Azi 1 deg Azimuth of blade 1
# NOTE THAT RADIUS FROM ch_details[ch, 0] REFERS TO THE RADIUS
# YOU ASKED FOR, AND ch_details[ch, 2] IS WHAT YOU GET, which is
# still based on a mean radius (deflections change the game)
elif self.ch_details[ch, 0].split(',')[0] in ch_aero:
dscr_list = self.ch_details[ch, 2].split(' ')
dscr_list = misc.remove_items(dscr_list, '')
blade_nr = self.ch_details[ch, 2].split('blade ')[1].split()[0]
# sometimes the units for aero sensors are wrong!
# there is no label option
# radius what you get
# radius = dscr_list[-1]
# radius what you asked for
tmp = self.ch_details[ch, 0].split('R=')
radius = misc.remove_items(tmp, '')[-1].strip()
# save all info in the dict
channelinfo = {}
channelinfo['sensortype'] = sensortype
channelinfo['radius'] = float(radius)
channelinfo['blade_nr'] = int(blade_nr)
channelinfo['units'] = units
channelinfo['chi'] = ch
# -----------------------------------------------------------------
# for the induction grid over the rotor
# a_grid, azi 0.00 r 1.74
elif self.ch_details[ch, 0].split(',')[0] in ch_aerogrid:
items = self.ch_details[ch, 0].split(',')
sensortype = items[0]
items2 = items[1].split(' ')
items2 = misc.remove_items(items2, '')
azi = items2[1]
# radius what you asked for
# and tag it
tag = '%s-azi-%s-r-%s' % (sensortype,azi,radius)
# save all info in the dict
channelinfo = {}
channelinfo['sensortype'] = sensortype
channelinfo['radius'] = float(radius)
channelinfo['azimuth'] = float(azi)
channelinfo['units'] = units
channelinfo['chi'] = ch
# -----------------------------------------------------------------
# INDUCTION AT THE BLADE
# 0: Induc. Vz, rpco, R= 1.4
# 1: m/s
# 2: Induced wsp Vz of blade 1 at radius 1.37, RP. coo.
# Induc. Vx, locco, R= 1.4
# Induced wsp Vx of blade 1 at radius 1.37, local ae coo.
# Induc. Vy, blco, R= 1.4
# Induced wsp Vy of blade 1 at radius 1.37, local bl coo.
# Induc. Vz, glco, R= 1.4
# Induced wsp Vz of blade 1 at radius 1.37, global coo.
# Induc. Vx, rpco, R= 8.4
# Induced wsp Vx of blade 1 at radius 8.43, RP. coo.
elif self.ch_details[ch, 0].strip()[:5] == 'Induc':
items = self.ch_details[ch, 2].split(' ')
items = misc.remove_items(items, '')
coord = self.ch_details[ch, 2].split(', ')[1].strip()
# radius what you get
# radius = float(items[8].replace(',', ''))
# radius what you asked for
tmp = self.ch_details[ch, 0].split(' ')
radius = float(misc.remove_items(tmp, '')[-1])
# and tag it
rpl = (coord, blade_nr, component, radius)
tag = 'induc-%s-blade-%1i-%s-r-%03.01f' % rpl
# save all info in the dict
channelinfo = {}
channelinfo['blade_nr'] = blade_nr
channelinfo['sensortype'] = 'induction'
channelinfo['radius'] = radius
channelinfo['coord'] = coord
channelinfo['component'] = component
channelinfo['units'] = units
channelinfo['chi'] = ch
# -----------------------------------------------------------------
# MORE AERO SENSORS
# Ae intfrc Fx, rpco, R= 0.0
# Aero int. force Fx of blade 1 at radius 0.00, RP coo.
# Ae secfrc Fy, R= 25.0
# Aero force Fy of blade 1 at radius 24.11
# Ae pos x, glco, R= 88.2
# Aero position x of blade 1 at radius 88.17, global coo.
elif self.ch_details[ch, 0].strip()[:2] == 'Ae':
units = self.ch_details[ch, 1]
items = self.ch_details[ch, 2].split(' ')
items = misc.remove_items(items, '')
# find blade number
tmp = self.ch_details[ch, 2].split('blade ')[1].strip()
blade_nr = int(tmp.split(' ')[0])
tmp = self.ch_details[ch, 2].split('radius ')[1].strip()
tmp = tmp.split(',')
# radius what you get
# radius = float(tmp[0])
# radius what you asked for
tmp = self.ch_details[ch, 0].split(' ')
radius = float(misc.remove_items(tmp, '')[-1])
if len(tmp) > 1:
coord = tmp[1].strip()
else:
coord = 'aero'
items = self.ch_details[ch, 0].split(' ')
sensortype = items[1]
component = items[2].replace(',', '')
# save all info in the dict
channelinfo = {}
channelinfo['blade_nr'] = blade_nr
channelinfo['sensortype'] = sensortype
channelinfo['radius'] = radius
channelinfo['coord'] = coord
channelinfo['component'] = component
channelinfo['units'] = units
channelinfo['chi'] = ch
rpl = (coord, blade_nr, sensortype, component, radius)
tag = 'aero-%s-blade-%1i-%s-%s-r-%03.01f' % rpl
# TODO: wind speed
# some spaces have been trimmed here
# WSP gl. coo.,Vy m/s
# // Free wind speed Vy, gl. coo, of gl. pos 0.00, 0.00, -2.31
# WSP gl. coo.,Vdir_hor deg
# Free wind speed Vdir_hor, gl. coo, of gl. pos 0.00, 0.00, -2.31
# -----------------------------------------------------------------
# WATER SURFACE gl. coo, at gl. coo, x,y= 0.00, 0.00
elif self.ch_details[ch, 2].startswith('Water'):
units = self.ch_details[ch, 1]
# but remove the comma
x = items[-2][:-1]
y = items[-1]
# and tag it
tag = 'watersurface-global-%s-%s' % (x, y)
# save all info in the dict
channelinfo = {}
channelinfo['coord'] = 'global'
channelinfo['pos'] = (float(x), float(y))
channelinfo['units'] = units
channelinfo['chi'] = ch
# -----------------------------------------------------------------
# WIND SPEED
# WSP gl. coo.,Vx
# Free wind speed Vx, gl. coo, of gl. pos 0.00, 0.00, -6.00 LABEL
elif self.ch_details[ch, 0].startswith('WSP gl.'):
units = self.ch_details[ch, 1]
direction = self.ch_details[ch, 0].split(',')[1]
tmp = self.ch_details[ch, 2].split('pos')[1]
x, y, z = tmp.split(',')
x, y, z = x.strip(), y.strip(), z.strip()
tmp = z.split(' ')
sensortag = ''
if len(tmp) == 2:
z, sensortag = tmp
elif len(tmp) == 1:
z = tmp[0]
# and tag it
tag = 'windspeed-global-%s-%s-%s-%s' % (direction, x, y, z)
# save all info in the dict
channelinfo = {}
channelinfo['coord'] = 'global'
channelinfo['pos'] = (x, y, z)
channelinfo['units'] = units
channelinfo['chi'] = ch
channelinfo['sensortag'] = sensortag
# FIXME: direction is the same as component, right?
channelinfo['direction'] = direction
channelinfo['sensortype'] = 'wsp-global'
# WIND SPEED AT BLADE
# 0: WSP Vx, glco, R= 61.5
# 2: Wind speed Vx of blade 1 at radius 61.52, global coo.
elif self.ch_details[ch, 0].startswith('WSP V'):
units = self.ch_details[ch, 1].strip()
tmp = self.ch_details[ch, 0].split(' ')[1].strip()
direction = tmp.replace(',', '')
blade_nr = self.ch_details[ch, 2].split('blade')[1].strip()[:2]
coord = self.ch_details[ch, 2].split(',')[1].strip()
blade_nr = blade_nr.strip()
# radius what you get
# radius = self.ch_details[ch, 2].split('radius')[1].split(',')[0]
# radius = radius.strip()
# radius what you asked for
tmp = self.ch_details[ch, 0].split(' ')
radius = misc.remove_items(tmp, '')[-1].strip()
# and tag it
rpl = (direction, blade_nr, radius, coord)
tag = 'wsp-blade-%s-%s-%s-%s' % rpl
# save all info in the dict
channelinfo = {}
channelinfo['coord'] = coord
# FIXME: direction is the same as component, right?
channelinfo['direction'] = direction
channelinfo['blade_nr'] = int(blade_nr)
channelinfo['radius'] = float(radius)
channelinfo['units'] = units
channelinfo['chi'] = ch
channelinfo['sensortype'] = 'wsp-blade'
# FLAP ANGLE
# 2: Flap angle for blade 3 flap number 1
elif self.ch_details[ch, 0][:7] == 'setbeta':
units = self.ch_details[ch, 1].strip()
blade_nr = self.ch_details[ch, 2].split('blade')[1].strip()
blade_nr = blade_nr.split(' ')[0].strip()
blade_nr = blade_nr.strip()
# and tag it
tag = 'setbeta-bladenr-%s-flapnr-%s' % (blade_nr, flap_nr)
# save all info in the dict
channelinfo = {}
channelinfo['flap_nr'] = int(flap_nr)
channelinfo['blade_nr'] = int(blade_nr)
channelinfo['units'] = units
channelinfo['chi'] = ch
# harmonic channel output
# Harmonic
# Harmonic sinus function
elif self.ch_details[ch, 0][:7] == 'Harmoni':
func_name = ' '.join(self.ch_details[ch, 1].split(' ')[1:])
channelinfo = {}
channelinfo['output_type'] = func_name
channelinfo['sensortype'] = 'harmonic'
channelinfo['chi'] = ch
base = self.ch_details[ch,2].strip().lower().replace(' ', '_')
tag = base + '_0'
if tag in self.ch_dict:
tag_nr = int(tag.split('_')[-1]) + 1
tag = base + '_%i' % tag_nr
# -----------------------------------------------------------------
# If all this fails, just combine channel name and description
tag = '-'.join(self.ch_details[ch,:3].tolist())
channelinfo = {}
channelinfo['chi'] = ch
channelinfo['units'] = self.ch_details[ch, 1].strip()
# -----------------------------------------------------------------
# add a v_XXX tag in case the channel already exists
if tag in self.ch_dict:
jj = 1
while True:
tag_new = tag + '_v%i' % jj
if tag_new in self.ch_dict:
jj += 1
else:
tag = tag_new
break
self.ch_dict[tag] = copy.copy(channelinfo)
# -----------------------------------------------------------------
# save in for DataFrame format
cols_ch = set(channelinfo.keys())
for col in cols_ch:
df_dict[col].append(channelinfo[col])
# the remainder columns we have not had yet. Fill in blank
for col in (self.cols - cols_ch):
df_dict[col].append('')
df_dict['unique_ch_name'].append(tag)
self.ch_df = pd.DataFrame(df_dict)
self.ch_df.set_index('chi', inplace=True)
def _ch_dict2df(self):
"""
Create a DataFrame version of the ch_dict, and the chi columns is
set as the index
"""
# identify all the different columns
cols = set()
for ch_name, channelinfo in self.ch_dict.items():
cols.update(set(channelinfo.keys()))
df_dict['unique_ch_name'] = []
for ch_name, channelinfo in self.ch_dict.items():
cols_ch = set(channelinfo.keys())
for col in cols_ch:
df_dict[col].append(channelinfo[col])
# the remainder columns we have not had yet. Fill in blank
for col in (cols - cols_ch):
df_dict[col].append('')
df_dict['unique_ch_name'].append(ch_name)
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
self.ch_df = pd.DataFrame(df_dict)
self.ch_df.set_index('chi', inplace=True)
def _data_window(self, nr_rev=None, time=None):
"""
Based on a time interval, create a proper slice object
======================================================
The window will start at zero and ends with the covered time range
of the time input.
Paramters
---------
nr_rev : int, default=None
NOT IMPLEMENTED YET
time : list, default=None
time = [time start, time stop]
Returns
-------
slice_
window
zoomtype
time_range
time_range = [0, time[1]]
"""
# -------------------------------------------------
# determine zome range if necesary
# -------------------------------------------------
time_range = None
if nr_rev:
raise NotImplementedError
# input is a number of revolutions, get RPM and sample rate to
# calculate the required range
# TODO: automatich detection of RPM channel!
time_range = nr_rev/(self.rpm_mean/60.)
# convert to indices instead of seconds
i_range = int(self.Freq*time_range)
window = [0, time_range]
# in case the first datapoint is not at 0 seconds
slice_ = np.r_[i_zero:i_range+i_zero]
zoomtype = '_nrrev_' + format(nr_rev, '1.0f') + 'rev'
elif time.any():
time_range = time[1] - time[0]
i_start = int(time[0]*self.Freq)
i_end = int(time[1]*self.Freq)
slice_ = np.r_[i_start:i_end]
window = [time[0], time[1]]
return slice_, window, zoomtype, time_range
def sig2df(self):
"""Convert sig to dataframe with unique channel names as column names.
"""
# channels that are not part of the naming scheme are not included
df = pd.DataFrame(self.sig[:,self.ch_df.index],
columns=self.ch_df['unique_ch_name'])
return df
# TODO: general signal method, this is not HAWC2 specific, move out
stats = {}
# calculate the statistics values:
stats['max'] = sig[i0:i1, :].max(axis=0)
stats['min'] = sig[i0:i1, :].min(axis=0)
stats['mean'] = sig[i0:i1, :].mean(axis=0)
stats['std'] = sig[i0:i1, :].std(axis=0)
stats['range'] = stats['max'] - stats['min']
stats['absmax'] = np.absolute(sig[i0:i1, :]).max(axis=0)
stats['rms'] = np.sqrt(np.mean(sig[i0:i1, :]*sig[i0:i1, :], axis=0))
stats['int'] = integrate.trapz(sig[i0:i1, :], x=sig[i0:i1, 0], axis=0)
1433
1434
1435
1436
1437
1438
1439
1440
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450
1451
1452
1453
1454
1455
1456
1457
1458
1459
1460
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
def statsdel_df(self, i0=0, i1=None, statchans='all', delchans='all',
m=[3, 4, 6, 8, 10, 12], neq=None, no_bins=46):
"""Calculate statistics and equivalent loads for the current loaded
signal.
Parameters
----------
i0 : int, default=0
i1 : int, default=None
channels : list, default='all'
all channels are selected if set to 'all', otherwise define a list
using the unique channel defintions.
neq : int, default=1
no_bins : int, default=46
Return
------
statsdel : pd.DataFrame
Pandas DataFrame with the statistical parameters and the different
fatigue coefficients as columns, and channels as rows. As index the
unique channel name is used.
"""
stats = ['max', 'min', 'mean', 'std', 'range', 'absmax', 'rms', 'int']
if statchans == 'all':
statchans = self.ch_df['unique_ch_name'].tolist()
statchis = self.ch_df['unique_ch_name'].index.values
else:
sel = self.ch_df['unique_ch_name']
statchis = self.ch_df[sel.isin(statchans)].index.values
if delchans == 'all':
delchans = self.ch_df['unique_ch_name'].tolist()
delchis = self.ch_df.index.values
else:
sel = self.ch_df['unique_ch_name']
delchis = self.ch_df[sel.isin(delchans)].index.values
# delchans has to be a subset of statchans!
if len(set(delchans) - set(statchans)) > 0:
raise ValueError('delchans has to be a subset of statchans')
tmp = np.ndarray((len(statchans), len(stats+m)))
tmp[:,:] = np.nan
m_cols = ['m=%i' % m_ for m_ in m]
statsdel = pd.DataFrame(tmp, columns=stats+m_cols)
statsdel.index = statchans
datasel = self.sig[i0:i1,statchis]
time = self.sig[i0:i1,0]
statsdel['max'] = datasel.max(axis=0)
statsdel['min'] = datasel.min(axis=0)
statsdel['mean'] = datasel.mean(axis=0)
statsdel['std'] = datasel.std(axis=0)
statsdel['range'] = statsdel['max'] - statsdel['min']
statsdel['absmax'] = np.abs(datasel).max(axis=0)
statsdel['rms'] = np.sqrt(np.mean(datasel*datasel, axis=0))
statsdel['int'] = integrate.trapz(datasel, x=time, axis=0)
statsdel['intabs'] = integrate.trapz(np.abs(datasel), x=time, axis=0)
if neq is None:
neq = self.sig[-1,0] - self.sig[0,0]
for chi, chan in zip(delchis, delchans):
signal = self.sig[i0:i1,chi]
eq = self.calc_fatigue(signal, no_bins=no_bins, neq=neq, m=m)
statsdel.loc[chan][m_cols] = eq
return statsdel
# TODO: general signal method, this is not HAWC2 specific, move out
def calc_fatigue(self, signal, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1):
"""
Parameters
----------
signal: 1D array
One dimentional array containing the signal.
no_bins: int
Number of bins for the binning of the amplitudes.
m: list
Values of the slope of the SN curve.
neq: int
Number of equivalent cycles
Returns
-------
eq: list
Damage equivalent loads for each m value.
return eq_load(signal, no_bins=no_bins, m=m, neq=neq)[0]
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
def cycle_matrix(self, signal, no_bins=46):
"""Cycle/Markov matrix.
Convenience function for wetb.fatigue_tools.fatigue.cycle_matrix2
Parameters
----------
signal: 1D array
One dimentional array containing the signal.
no_bins: int
Number of bins for the binning of the amplitudes.
Returns
-------
cycles : ndarray, shape(ampl_bins, mean_bins)
A bi-dimensional histogram of load cycles(full cycles). Amplitudes
are histogrammed along the first dimension and mean values are
histogrammed along the second dimension.
ampl_edges : ndarray, shape(no_bins+1,n)
The amplitude bin edges
mean_edges : ndarray, shape(no_bins+1,n)
The mean bin edges
"""
return cycle_matrix2(signal, no_bins)
def blade_deflection(self):
"""
"""
# select all the y deflection channels
db = misc.DictDB(self.ch_dict)
db.search({'sensortype': 'state pos', 'component': 'z'})
# sort the keys and save the mean values to an array/list
chiz, zvals = [], []
for key in sorted(db.dict_sel.keys()):
zvals.append(-self.sig[:, db.dict_sel[key]['chi']].mean())
chiz.append(db.dict_sel[key]['chi'])
# sort the keys and save the mean values to an array/list
chiy, yvals = [], []
for key in sorted(db.dict_sel.keys()):
yvals.append(self.sig[:, db.dict_sel[key]['chi']].mean())
chiy.append(db.dict_sel[key]['chi'])
return np.array(zvals), np.array(yvals)
David Verelst
committed
def save_chan_names(self, fname):
"""Save unique channel names to text file.
"""
channels = self.ch_df.ch_name.values
channels.sort()
np.savetxt(fname, channels, fmt='%-100s')
def save_channel_info(self, fname):
"""Save all channel info: unique naming + HAWC2 description from *.sel.
"""
p1 = self.ch_df.copy()
# but ignore the units column, we already have that
p2 = pd.DataFrame(self.ch_details,
columns=['Description1', 'units', 'Description2'])
# merge on the index
tmp = pd.merge(p1, p2, right_index=True, how='outer', left_index=True)
tmp.to_excel(fname)
# for a fixed-with text format instead of csv
# header = ''.join(['%100s' % k for k in tmp.columns])
# header = ' windspeed' + header
# np.savetxt(fname, tmp.to_records(), header=header,
# fmt='% 01.06e ')
return tmp
David Verelst
committed
def load_chan_names(self, fname):
dtype = np.dtype('U100')
return np.genfromtxt(fname, dtype=dtype, delimiter=';').tolist()
def save_csv(self, fname, fmt='%.18e', delimiter=','):
"""
Save to csv and use the unified channel names as columns
"""
map_sorting = {}
# first, sort on channel index
for ch_key, ch in self.ch_dict.items():
map_sorting[ch['chi']] = ch_key
header = []
# not all channels might be present...iterate again over map_sorting
for chi in map_sorting:
try:
sensortag = self.ch_dict[map_sorting[chi]]['sensortag']
header.append(map_sorting[chi] + ' // ' + sensortag)
except:
header.append(map_sorting[chi])
# and save
print('saving...', end='')
np.savetxt(fname, self.sig[:, list(map_sorting.keys())], fmt=fmt,
delimiter=delimiter, header=delimiter.join(header))
print(fname)
def save_df(self, fname):
"""
Save the HAWC2 data and sel file in a DataFrame that contains all the
data, and all the channel information (the one from the sel file
and the parsed from this function)
"""
self.sig
self.ch_details
self.ch_dict
def ReadOutputAtTime(fname):
"""Distributed blade loading as generated by the HAWC2 output_at_time
command. From HAWC2 12.3-beta and onwards, there are 7 header columns,
earlier version only have 3.
Parameters
----------
fname : str
header_lnr : int, default=3
Line number of the header (column names) (1-based counting).
# data = pd.read_fwf(fname, skiprows=3, header=None)
# pd.read_table(fname, sep=' ', skiprows=3)
# data.index.names = cols
# because the formatting is really weird, we need to sanatize it a bit
with opent(fname, 'r') as f:
# read the header from line 3
for k in range(7):
line = f.readline()
if line[0:12].lower().replace('#', '').strip() == 'radius_s':
header_lnr = k + 1
break
header = line.replace('\r', '').replace('\n', '')
cols = [k.strip().replace(' ', '_') for k in header.split('#')[1:]]
data = np.loadtxt(fname, skiprows=header_lnr)
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
return pd.DataFrame(data, columns=cols)
def ReadEigenBody(fname, debug=False):
"""
Read HAWC2 body eigenalysis result file
=======================================
Parameters
----------
file_path : str
file_name : str
Returns
-------
results : DataFrame
Columns: body, Fd_hz, Fn_hz, log_decr_pct
"""
# Body data for body number : 3 with the name :nacelle
# Results: fd [Hz] fn [Hz] log.decr [%]
# Mode nr: 1: 1.45388E-21 1.74896E-03 6.28319E+02
FILE = opent(fname)
lines = FILE.readlines()
FILE.close()
df_dict = {'Fd_hz': [], 'Fn_hz': [], 'log_decr_pct': [], 'body': []}
for i, line in enumerate(lines):
if debug: print('line nr: %5i' % i)
# identify for which body we will read the data
if line[:25] == 'Body data for body number':
body = line.split(':')[2].rstrip().lstrip()
# remove any annoying characters
if debug: print('modes for body: %s' % body)
# identify mode number and read the eigenfrequencies
elif line[:8] == 'Mode nr:':
linelist = line.replace('\n', '').replace('\r', '').split(':')
# modenr = linelist[1].rstrip().lstrip()
# text after Mode nr can be empty
try:
eigenmodes = linelist[2].rstrip().lstrip().split(' ')
except IndexError:
eigenmodes = ['0', '0', '0']
if debug: print(eigenmodes)
# in case we have more than 3, remove all the empty ones
# this can happen when there are NaN values
if not len(eigenmodes) == 3:
eigenmodes = linelist[2].rstrip().lstrip().split(' ')
eigmod = []
for k in eigenmodes:
if len(k) > 1:
eigmod.append(k)
else:
eigmod = eigenmodes
# remove any trailing spaces for each element
for k in range(len(eigmod)):
1747
1748
1749
1750
1751
1752
1753
1754
1755
1756
1757
1758
1759
1760
1761
1762
1763
1764
1765
1766
1767
1768
1769
1770
1771
1772
1773
1774
1775
1776
1777
1778
1779
1780
1781
1782
1783
1784
1785
1786
1787
1788
1789
1790
1791
1792
1793
1794
1795
df_dict['body'].append(body)
df_dict['Fd_hz'].append(eigmod[0])
df_dict['Fn_hz'].append(eigmod[1])
df_dict['log_decr_pct'].append(eigmod[2])
return pd.DataFrame(df_dict)
def ReadEigenStructure(file_path, file_name, debug=False, max_modes=500):
"""
Read HAWC2 structure eigenalysis result file
============================================
The file looks as follows:
#0 Version ID : HAWC2MB 11.3
#1 ___________________________________________________________________
#2 Structure eigenanalysis output
#3 ___________________________________________________________________
#4 Time : 13:46:59
#5 Date : 28:11.2012
#6 ___________________________________________________________________
#7 Results: fd [Hz] fn [Hz] log.decr [%]
#8 Mode nr: 1: 3.58673E+00 3.58688E+00 5.81231E+00
#...
#302 Mode nr:294: 0.00000E+00 6.72419E+09 6.28319E+02
Parameters
----------
file_path : str
file_name : str
debug : boolean, default=False
max_modes : int
Stop evaluating the result after max_modes number of modes have been
identified
Returns
-------
modes_arr : ndarray(3,n)
An ndarray(3,n) holding Fd, Fn [Hz] and the logarithmic damping
decrement [%] for n different structural eigenmodes
"""
# 0 Version ID : HAWC2MB 11.3
# 1 ___________________________________________________________________
# 2 Structure eigenanalysis output
# 3 ___________________________________________________________________
# 4 Time : 13:46:59
# 5 Date : 28:11.2012
# 6 ___________________________________________________________________
# 7 Results: fd [Hz] fn [Hz] log.decr [%]
# 8 Mode nr: 1: 3.58673E+00 3.58688E+00 5.81231E+00
# Mode nr:294: 0.00000E+00 6.72419E+09 6.28319E+02
FILE = opent(os.path.join(file_path, file_name))
lines = FILE.readlines()
FILE.close()
header_lines = 8
# we now the number of modes by having the number of lines
nrofmodes = len(lines) - header_lines
for i, line in enumerate(lines):
if i > max_modes:
# cut off the unused rest
break
# ignore the header
if i < header_lines:
continue
# split up mode nr from the rest
parts = line.split(':')
# get fd, fn and damping, but remove all empty items on the list
modes_arr[:, i-header_lines]=misc.remove_items(parts[2].split(' '), '')
return modes_arr
"""
"""
def __init__(self):
pass
def __call__(self, z_h, r_blade_tip, a_phi=None, shear_exp=None, nr_hor=3,
nr_vert=20, h_ME=500.0, io=None, wdir=None):
1846
1847
1848
1849
1850
1851
1852
1853
1854
1855
1856
1857
1858
1859
1860
1861
1862
1863
1864
1865
1866
1867
1868
1869
1870
"""
Parameters
----------
z_h : float
Hub height
r_blade_tip : float
Blade tip radius
a_phi : float, default=None
:math:`a_{\\varphi} \\approx 0.5` parameter for the modified
Ekman veer distribution. Values vary between -1.2 and 0.5.
shear_exp : float, default=None
nr_vert : int, default=3
nr_hor : int, default=20
h_ME : float, default=500
Modified Ekman parameter. Take roughly 500 for off shore sites,
1000 for on shore sites.
io : str or io buffer, default=None
When specified, the HAWC2 user defined shear input file will be
written.
wdir : float, default=None
A constant veer angle, or yaw angle. Equivalent to setting the
wind direction. Angle in degrees.
Returns
-------
uu, vv, ww, xx, zz
"""
x, z = self.create_coords(z_h, r_blade_tip, nr_vert=nr_vert,
nr_hor=nr_hor)
if a_phi is not None:
phi_rad = WindProfiles.veer_ekman_mod(z, z_h, h_ME=h_ME, a_phi=a_phi)
assert len(phi_rad) == nr_vert
else:
nr_vert = len(z)
phi_rad = np.zeros((nr_vert,))
# add any yaw error on top of
if wdir is not None:
# because wdir cw positive, and phi veer ccw positive
phi_rad -= (wdir*np.pi/180.0)
u, v, w = self.decompose_veer(phi_rad, nr_hor)
# when no shear is defined
if shear_exp is None:
uu = u
vv = v
ww = w
else:
# scale the shear on top of the veer
shear = WindProfiles.powerlaw(z, z_h, shear_exp)
uu = u*shear[:,np.newaxis]
vv = v*shear[:,np.newaxis]
ww = w*shear[:,np.newaxis]
# and write to a file
if isinstance(io, str):
with open(io, 'wb') as fid:
fid = self.write(fid, uu, vv, ww, x, z)
self.fid =None
elif io is not None:
io = self.write(io, uu, vv, ww, x, z)
self.fid = io
return uu, vv, ww, x, z
def create_coords(self, z_h, r_blade_tip, nr_vert=3, nr_hor=20):
"""
Utility to create the coordinates of the wind field based on hub heigth
and blade length. Add 15% to r_blade_tip to make sure horizontal edges
are defined wide enough.
"""
# take 15% extra space after the blade tip
z = np.linspace(0, z_h + r_blade_tip*1.15, nr_vert)
# along the horizontal, coordinates with 0 at the rotor center
x = np.linspace(-r_blade_tip*1.15, r_blade_tip*1.15, nr_hor)
return x, z
def deltaphi2aphi(self, d_phi, z_h, r_blade_tip, h_ME=500.0):
"""For a given `\\Delta \\varphi` over the rotor diameter, estimate
the corresponding `a_{\\varphi}`.
Parameters
----------
`\\Delta \\varphi` : ndarray or float
Veer angle difference over the rotor plane from lowest to highest
blade tip position.
z_h : float
Hub height in meters.
r_blade_tip : float
Blade tip radius/length.
h_ME : float, default=500.0
Modified Ekman parameter. For on shore,
:math:`h_{ME} \\approx 1000`, for off-shore,
:math:`h_{ME} \\approx 500`
Returns
-------
`a_{\\varphi}` : ndarray or float
"""
t1 = r_blade_tip * 2.0 * np.exp(-z_h/(h_ME))
a_phi = d_phi * np.sqrt(h_ME*z_h) / t1
return a_phi
def deltaphi2aphi_opt(self, deltaphi, z, z_h, r_blade_tip, h_ME):
"""
convert delta_phi over a given interval z to a_phi using
scipy.optimize.fsolve on veer_ekman_mod.
Parameters
----------
deltaphi : float
Desired delta phi in rad over interval z[0] at bottom to z[1] at
the top.
def func(a_phi, z, z_h, h_ME, deltaphi_target):
phis = WindProfiles.veer_ekman_mod(z, z_h, h_ME=h_ME, a_phi=a_phi)
return np.abs(deltaphi_target - (phis[1] - phis[0]))
args = (z, z_h, h_ME, deltaphi)
return sp.optimize.fsolve(func, [0], args=args)[0]
def decompose_veer(self, phi_rad, nr_hor):
"""
Convert a veer angle into u, v, and w components, ready for the
HAWC2 user defined veer input file. nr_vert refers to the number of
vertical grid points.
Paramters
---------
phi_rad : ndarray(nr_vert)
veer angle in radians as function of height