Newer
Older
'eq_load' calculate equivalent loads using one of the two rain flow counting methods
'cycle_matrix' calculates a matrix of cycles (binned on amplitude and mean value)
'eq_load_and_cycles' is used to calculate eq_loads of multiple time series (e.g. life time equivalent load)
The methods uses the rainflow counting routines (See documentation in top of methods):
- 'rainflow_windap': (Described in "Recommended Practices for Wind Turbine Testing - 3. Fatigue Loads",
2. edition 1990, Appendix A)
or
- 'rainflow_astm' (based on the c-implementation by Adam Nieslony found at the MATLAB Central File Exchange
http://www.mathworks.com/matlabcentral/fileexchange/3026)
from __future__ import division
from __future__ import print_function
from __future__ import unicode_literals
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from wetb.fatigue_tools.rainflowcounting import rainflowcount
rainflow_windap = rainflowcount.rainflow_windap
rainflow_astm = rainflowcount.rainflow_astm
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
def eq_load(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=1, rainflow_func=rainflow_windap):
"""Equivalent load calculation
Calculate the equivalent loads for a list of Wohler exponent and number of equivalent loads
Parameters
----------
signals : list of tuples or array_like
- if list of tuples: list must have format [(sig1_weight, sig1),(sig2_weight, sig1),...] where\n
- sigx_weight is the weight of signal x\n
- sigx is signal x\n
- if array_like: The signal
no_bins : int, optional
Number of bins in rainflow count histogram
m : int, float or array-like, optional
Wohler exponent (default is [3, 4, 6, 8, 10, 12])
neq : int, float or array-like, optional
The equivalent number of load cycles (default is 1, but normally the time duration in seconds is used)
rainflow_func : {rainflow_windap, rainflow_astm}, optional
The rainflow counting function to use (default is rainflow_windap)
Returns
-------
eq_loads : array-like
List of lists of equivalent loads for the corresponding equivalent number(s) and Wohler exponents
Examples
--------
>>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
>>> eq_load(signal, no_bins=50, neq=[1, 17], m=[3, 4, 6], rainflow_func=rainflow_windap)
[[10.311095426959747, 9.5942535021382174, 9.0789213365013932], # neq = 1, m=[3,4,6]
[4.010099657859783, 4.7249689509841746, 5.6618639965313005]], # neq = 17, m=[3,4,6]
eq_load([(.4, signal), (.6, signal)], no_bins=50, neq=[1, 17], m=[3, 4, 6], rainflow_func=rainflow_windap)
[[10.311095426959747, 9.5942535021382174, 9.0789213365013932], # neq = 1, m=[3,4,6]
[4.010099657859783, 4.7249689509841746, 5.6618639965313005]], # neq = 17, m=[3,4,6]
"""
try:
return eq_load_and_cycles(signals, no_bins, m, neq, rainflow_func)[0]
except TypeError:
return [[np.nan] * len(np.atleast_1d(m))] * len(np.atleast_1d(neq))
def eq_load_and_cycles(signals, no_bins=46, m=[3, 4, 6, 8, 10, 12], neq=[10 ** 6, 10 ** 7, 10 ** 8], rainflow_func=rainflow_windap):
"""Calculate combined fatigue equivalent load
Parameters
----------
signals : list of tuples or array_like
- if list of tuples: list must have format [(sig1_weight, sig1),(sig2_weight, sig1),...] where\n
- sigx_weight is the weight of signal x\n
- sigx is signal x\n
- if array_like: The signal
no_bins : int, optional
Number of bins for rainflow counting
m : int, float or array-like, optional
Wohler exponent (default is [3, 4, 6, 8, 10, 12])
neq : int or array-like, optional
Equivalent number, default is [10^6, 10^7, 10^8]
rainflow_func : {rainflow_windap, rainflow_astm}, optional
The rainflow counting function to use (default is rainflow_windap)
Returns
-------
eq_loads : array-like
List of lists of equivalent loads for the corresponding equivalent number(s) and Wohler exponents
cycles : array_like
2d array with shape = (no_ampl_bins, no_mean_bins)
ampl_bin_mean : array_like
mean amplitude of the bins
ampl_bin_edges
Edges of the amplitude bins
"""
cycles, ampl_bin_mean, ampl_bin_edges, _, _ = cycle_matrix(signals, no_bins, 1, rainflow_func)
if 0: #to be similar to windap
ampl_bin_mean = (ampl_bin_edges[:-1] + ampl_bin_edges[1:]) / 2
cycles, ampl_bin_mean = cycles.flatten(), ampl_bin_mean.flatten()
eq_loads = [[((np.sum(cycles * ampl_bin_mean ** _m) / _neq) ** (1. / _m)) for _m in np.atleast_1d(m)] for _neq in np.atleast_1d(neq)]
return eq_loads, cycles, ampl_bin_mean, ampl_bin_edges
def cycle_matrix(signals, ampl_bins=10, mean_bins=10, rainflow_func=rainflow_windap):
"""Markow load cycle matrix
Calculate the Markow load cycle matrix
Parameters
----------
Signals : array-like or list of tuples
if array-like, the raw signal
if list of tuples, list of (weight, signal), e.g. [(0.1,sig1), (0.8,sig2), (.1,sig3)]
ampl_bins : int or array-like, optional
if int, Number of amplitude value bins (default is 10)
if array-like, the bin edges for amplitude
mean_bins : int or array-like, optional
if int, Number of mean value bins (default is 10)
if array-like, the bin edges for mea
rainflow_func : {rainflow_windap, rainflow_astm}, optional
The rainflow counting function to use (default is rainflow_windap)
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.
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
ampl_bin_mean : ndarray, shape(ampl_bins,)
The average cycle amplitude of the bins
ampl_edges : ndarray, shape(ampl_bins+1,)
The amplitude bin edges
mean_bin_mean : ndarray, shape(ampl_bins,)
The average cycle mean of the bins
mean_edges : ndarray, shape(mean_bins+1,)
The mean bin edges
Examples
--------
>>> signal = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
>>> cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = cycle_matrix(signal)
>>> cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges = cycle_matrix([(.4, signal), (.6,signal)])
"""
if isinstance(signals[0], tuple):
ampls = np.empty((0,), dtype=np.float64)
means = np.empty((0,), dtype=np.float64)
weights = np.empty((0,), dtype=np.float64)
for w, signal in signals:
a, m = rainflow_func(signal[:])
ampls = np.r_[ampls, a]
means = np.r_[means, m]
weights = np.r_[weights, (np.zeros_like(a) + w)]
else:
ampls, means = rainflow_func(signals[:])
weights = np.ones_like(ampls)
if isinstance(ampl_bins, int):
ampl_bins = np.linspace(0, 1, num=ampl_bins + 1) * ampls.max()
cycles, ampl_edges, mean_edges = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights)
ampl_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * ampls)[0]
ampl_bin_mean = np.zeros_like(cycles)
mask = (cycles > 0)
ampl_bin_mean[mask] = ampl_bin_sum[mask] / cycles[mask]
mean_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * means)[0]
mean_bin_mean = np.zeros_like(cycles)
mean_bin_mean[cycles > 0] = mean_bin_sum[cycles > 0] / cycles[cycles > 0]
cycles = cycles / 2 # to get full cycles
return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
if __name__ == "__main__":
signal1 = np.array([-2.0, 0.0, 1.0, 0.0, -3.0, 0.0, 5.0, 0.0, -1.0, 0.0, 3.0, 0.0, -4.0, 0.0, 4.0, 0.0, -2.0])
signal2 = signal1 * 1.1
print (eq_load(signal1, no_bins=50, neq=17, rainflow_func=rainflow_windap))
print (eq_load(signal1, no_bins=50, neq=17, rainflow_func=rainflow_astm))
# Cycle matrix with 4 amplitude bins and 4 mean value bins
print (cycle_matrix(signal1, 4, 4, rainflow_func=rainflow_windap))
print (cycle_matrix(signal1, 4, 4, rainflow_func=rainflow_astm))
# Cycle matrix where signal1 and signal2 contributes with 50% each
print (cycle_matrix([(.5, signal1), (.5, signal2)], 4, 8, rainflow_func=rainflow_astm))