How to determe the Fatigue Damage according Miner's Rule via rainflow counting

54 views (last 30 days)
Tyl
Tyl on 28 Mar 2024
Answered: AH on 19 Apr 2024 at 9:07
Dear Matlab Users,
I have difficulties calculating the total accumulated fatique damage for a set of data (see Attachements). I have determined the histogram based on Rainflow counting, and determined the line for the SN-curve (maximum allowable stress per cycle). The SN-curve has only one slope in this case.
To calculate the total accumulated fatique damage according Miner's Rule, the number of cycles of the SN curve should be divided by the number of cycles of the histogram per stress range, and then be summed up.
My question is: How do I divide the values from the number of cycles of the SN-curve by the the values from the number of cycles of the histogram? How do I make the values of the same size? I want to divide the following (see script for clarification):
N=10.^(L_a1-m1*log(cigma_delta*(T/T_ref)^k)); with cigma_delta= [0:0.05:710]; % N-cycles SN-curve
by
C_7=histogram('BinEdges',edges','BinCounts',sum(hist,2), Orientation='horizontal') % Rainflow counting histogram
Thanks in advance!

Sign in to answer this question.

Accepted Answer

VBBV
VBBV on 28 Mar 2024 at 11:44
Edited: VBBV on 28 Mar 2024 at 12:08
Ran in:
Use linspace function in this line and evaluate the damage
cigma_delta= linspace(1e-3,2*Yield,length(C_7.Values));
clear all, close all, clc
%% Hotspot stress rainflow counting
Data = readtable('Data_hotspotstress2.txt')
Data = 12000×8 table
Var1 Var2 Var3 Var4 Var5 Var6 Var7 Var8 ____ _____ ______ _____ _____ _____ _____ _____ 4 -71.6 -129 -99.5 -35.4 -9.4 -1.5 18.5 3.2 -67.6 -123 -94.8 -35.2 -14.7 -9.7 12.5 2.4 -63.3 -116.4 -89.6 -34.9 -20.3 -18.4 6 2 -63.5 -116.9 -89.8 -35.3 -22.2 -21.2 4.2 2 -67.7 -123.6 -94.7 -36.2 -20.6 -18.7 6.5 1.8 -70.8 -128.6 -98.3 -37 -19.9 -17.8 7.5 1.3 -69.1 -126 -96.2 -37.1 -23.4 -23.1 3.7 0.5 -64.4 -119 -90.8 -36.8 -29.2 -31.9 -2.8 0 -62.2 -115.7 -88.2 -36.8 -32.8 -37.4 -6.8 -0.2 -64.7 -119.8 -91.1 -37.6 -32.6 -37.1 -6.2 -0.3 -68.9 -126.4 -95.9 -38.6 -31.5 -35.4 -4.5 -1 -69.8 -128.1 -96.8 -39.1 -33.5 -38.6 -6.5 -2 -66.8 -123.5 -93.1 -39.1 -39.1 -47.2 -12.8 -2.8 -63.6 -118.8 -89.2 -39.1 -44.6 -55.6 -18.9 -3.3 -64.3 -120.1 -90 -39.7 -46.8 -58.9 -21 -3.5 -67.9 -126 -94.2 -40.8 -46.4 -58.3 -20
cigma_7= table2array(Data(:,[7]));
fs = 20;
t = seconds(0:1/fs:600-1/fs)';
%Max-min stress range
Range_cigma7=abs(max(cigma_7))+abs(min(cigma_7));
% Stress history
figure(1)
plot(t,cigma_7 )
title('Stress history')
xlabel('Time [s]')
ylabel('Stress [MPa]')
%% Rainflow counting: plots number of cycles occurring per stress
%Histogram rainflow counting
figure(2)
[C7,hist,edges] = rainflow(cigma_7','ext');
Warning: Provided extrema do not produce strictly alternating peaks and valleys.
C_7=histogram('BinEdges',edges','BinCounts',sum(hist,2), Orientation='horizontal')
C_7 =
Histogram with properties: Data: [] Values: [1.2775e+03 1123 864.5000 666 533.5000 403.5000 253.5000 207 144 102 58.5000 55 44.5000 32 21 19 10.5000 8 10 8 6.5000 9.5000 3.5000 1.5000 2.5000 5 4 2 1.5000 … ] (1×124 double) NumBins: 124 BinEdges: [0 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 … ] (1×125 double) BinWidth: 1 BinLimits: [0 124] Normalization: 'count' FaceColor: 'auto' EdgeColor: [0 0 0] Use GET to show all properties
set(gca,'xscale','log')
set(gca,'yscale','log')
title('Rainflow counting for \sigma_7')
subtitle(['with a max/min stress difference of ', num2str(Range_cigma7), ' [Mpa]'])
ylim([0 1000])
xlabel('Cycle Counts log(N)')
ylabel('Stress Range log(S)')
%Straight line to visualise behaviour histogram:
hold on;
line([1, 5000], [180, 2], 'LineWidth', 1, 'Color', 'b');
%% SN-curve: plots maximum allowable number of cycles of each stress range
% The total stress fluctuation (i.e. maximum compression and maximum tension) should be considered to be
% transmitted through the welds for fatigue assessments.
% For N< 10^6 cycles and SCF<10
m1= 3; % slope, Hotspot detail parameter
L_a1= 11.764; % log a1
k= 0.25; % Thickness component
Yield= 355 % Yield stress steel
Yield = 355
cigma_delta= linspace(1e-3,2*Yield,length(C_7.Values)); % Stress range in MPa
T= 0.015875; % Thickness through which a crack will most likely grow
T_ref= 0.032; % Reference thickness tubular joints
N=10.^((L_a1-log(cigma_delta*(T/T_ref)^k))/m1) %Number of cycles
N = 1×124
1.0e+06 * 1.9157 0.0025 0.0015 0.0011 0.0009 0.0007 0.0006 0.0006 0.0005 0.0005 0.0004 0.0004 0.0004 0.0003 0.0003 0.0003 0.0003 0.0003 0.0003 0.0003 0.0002 0.0002 0.0002 0.0002 0.0002 0.0002 0.0002 0.0002 0.0002 0.0002
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
% SN-curve hotspot detail
figure(3)
plot(N,cigma_delta)
set(gca,'xscale','log')
set(gca,'yscale','log')
title('SN-curve for \sigma_7')
ylim([0 1000])
xlabel('Cycle Counts log(N)')
ylabel('Stress Range log(S)')
% Combined plot SN-curve and stress histogram
figure(4)
C_7=histogram('BinEdges',edges','BinCounts',sum(hist,2), Orientation='horizontal')
C_7 =
Histogram with properties: Data: [] Values: [1.2775e+03 1123 864.5000 666 533.5000 403.5000 253.5000 207 144 102 58.5000 55 44.5000 32 21 19 10.5000 8 10 8 6.5000 9.5000 3.5000 1.5000 2.5000 5 4 2 1.5000 … ] (1×124 double) NumBins: 124 BinEdges: [0 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 … ] (1×125 double) BinWidth: 1 BinLimits: [0 124] Normalization: 'count' FaceColor: 'auto' EdgeColor: [0 0 0] Use GET to show all properties
hold on;
line([1, 5000], [180, 2], 'LineWidth', 1, 'Color', 'b');
hold on
plot(N,cigma_delta, 'Color', 'r')
set(gca,'xscale','log')
set(gca,'yscale','log')
title('SN-curve for \sigma_7')
xlim([1 10^7])
ylim([1 1000])
xlabel('Cycle Counts log(N)')
ylabel('Stress Range log(S)')
%% Fatigue Damage: Miner's sum
counts=histcounts(C7);
N./C_7.Values
ans = 1×124
1.0e+03 * 1.4996 0.0022 0.0017 0.0016 0.0016 0.0018 0.0025 0.0027 0.0035 0.0045 0.0073 0.0072 0.0083 0.0108 0.0156 0.0164 0.0282 0.0353 0.0270 0.0324 0.0384 0.0253 0.0662 0.1493 0.0867 0.0420 0.0510 0.0990 0.1284 0.0937
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
  1 Comment
Tyl
Tyl on 29 Mar 2024 at 16:53
Edited: Tyl on 29 Mar 2024 at 17:32
Thank you for your answer!
The linspace does make sense, but I do have follow up question. When I change the values for stress range (Yield in cigma_delta) in the SN curve, the curve shifts towards this value (see figure attachment), whereas 2*Yield should actually be the interference stress at N=1. Do you know how I resolve this?
Also, N=10.^((L_a1-log(cigma_delta*(T/T_ref)^k))/m1) should be N=10.^((L_a1-m*log10(cigma_delta*(T/T_ref)^k))). I changed this in the attached file.

Sign in to comment.

More Answers (2)

Strider
Strider on 28 Mar 2024 at 20:47
Ran in:
If you have X and Y coordinates of your SN curve you can intercept the curve for whatever stress values you have.
close all
% make a fake SN curve
X = 1000:1000:10e6;
Y = exp(flip(linspace(1,10,numel(X))));
semilogx(X,Y);
hold on
% intercept curve at a stress that you want
% you can do any amount of binning you want before this.
% wrap a function around this to go through all the bins you have
s = [20000]; % stress I know "Little n" at from rainflow
n = 100;
yq= sort([Y s]);
v = interp1(Y,X,yq,'linear');
idx = ismember(yq,s); % find where it is
N = v(idx); % "Big N"
scatter(N,s) % plot it to check intercept
d_i = n/N
d_i = 9.2400e-04
AH
AH on 19 Apr 2024 at 9:07

You may want to take a look at the example Practical Introduction to Fatigue Analysis Using Rainflow Counting. In this example hysteresis and peak-valley filtering is aaplied to the dat as pre-processing steps before performing the rainflow counting.

Categories

Find more on Vibration Analysis in Help Center and File Exchange

Tags

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!