-
Notifications
You must be signed in to change notification settings - Fork 7
/
Copy pathfatigue.py
255 lines (212 loc) · 10.9 KB
/
fatigue.py
1
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
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
133
134
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
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
'''
Created on 04/03/2013
@author: mmpe
'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)
'''
import warnings
import numpy as np
from wetb.fatigue_tools.rainflowcounting import rainflowcount
rainflow_windap = rainflowcount.rainflow_windap
rainflow_astm = rainflowcount.rainflow_astm
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, 1)
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()
with warnings.catch_warnings():
warnings.simplefilter("ignore")
eq_loads = [[((np.nansum(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\n
- if list of tuples, list of (weight, signal), e.g. [(0.1,sig1), (0.8,sig2), (.1,sig3)]\n
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.
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):
weights, ampls, means = np.array([(np.zeros_like(ampl) + weight, ampl, mean) for weight,
signal in signals for ampl, mean in rainflow_func(signal[:]).T], dtype=np.float64).T
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[weights > 0].max()
cycles, ampl_edges, mean_edges = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights)
with warnings.catch_warnings():
warnings.simplefilter("ignore")
ampl_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * ampls)[0]
ampl_bin_mean = np.nanmean(ampl_bin_sum / np.where(cycles, cycles, np.nan), 1)
mean_bin_sum = np.histogram2d(ampls, means, [ampl_bins, mean_bins], weights=weights * means)[0]
mean_bin_mean = np.nanmean(mean_bin_sum / np.where(cycles, cycles, np.nan), 1)
cycles = cycles / 2 # to get full cycles
return cycles, ampl_bin_mean, ampl_edges, mean_bin_mean, mean_edges
def cycle_matrix2(signal, nrb_amp, nrb_mean, rainflow_func=rainflow_windap):
"""
Same as wetb.fatigue_tools.fatigue.cycle_matrix but bin from min_amp to
max_amp instead of 0 to max_amp.
Parameters
----------
Signal : ndarray(n)
1D Raw signal array
nrb_amp : int
Number of bins for the amplitudes
nrb_mean : int
Number of bins for the means
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.
ampl_edges : ndarray, shape(no_bins+1,n)
The amplitude bin edges
mean_edges : ndarray, shape(no_bins+1,n)
The mean bin edges
"""
bins = [nrb_amp, nrb_mean]
ampls, means = rainflow_func(signal)
weights = np.ones_like(ampls)
cycles, ampl_edges, mean_edges = np.histogram2d(ampls, means, bins,
weights=weights)
cycles = cycles / 2 # to get full cycles
return cycles, ampl_edges, mean_edges
def eq_loads_from_Markov(markov_matrix, m_list, neq=1e7):
"""
Calculate the damage equivalent loads from the cycles and amplitudes
of a Markov matrix for different Woehler slopes.
Parameters
----------
markov_matrix : array (Nbins x 2)
Array containing the number of cycles and mean ampltude for each bin.
Can have extra leading dimensions such as Nsensors x Ndimensions
m_list : list (Nwoehlerslopes)
List containing different Woehler slopes.
neq : int or float, optional
Number of equivalent load cycles. The default is 1e7.
Returns
-------
Array (Nwoehlerslopes)
Array containing the equivalent loads for each Woehler slope and
for each extra leading dimension if any
"""
cycles = markov_matrix[..., 0]
amplitudes = markov_matrix[..., 1]
eq_loads = np.array([((np.nansum(cycles*amplitudes**m, axis=-1)/neq)**(1/m))
for m in m_list])
eq_loads = np.transpose(eq_loads, list(range(1, eq_loads.ndim)) + [0])
return eq_loads
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
# equivalent load for default wohler slopes
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))