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

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!

Answer 1

[VBBV]

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>

[Tyl]

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.

Answer 2

[Strider]

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

Answer 3

[AH]

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.