Help with lumping of fatigue data

Hello everyone, I'm coming to you for help with a piece of code I need to write. As part of an assignment I need to perform the lumping of some fatigue data coming from a FE analysis of a rotating shaft. The txt file I've been provided contains the time instants of the analysis, the stess amplitude and the mean stress. I have calculated the equivalent stress with Haigh's formula and since I know the rotating velocity of the shaft I've calculated a vector with the cycles in between two instants. Now as for the lumped histogram I want it to gave the equivalent stress on the y axis and the number of fatigue cycles on the x axis. I have written something but I'm not sure it is implementing correctly the cycles, so I've come here for help. Please tell me if you can find any mistakes and have any suggestions. Thank you very much in advance. I'll paste here the code that I'm using right now:

profile1 = load('profile1.txt');
time_inst = profile1(:,1);
Sa = profile1(:,2);
Sm = profile1(:,3);
omega = 2*pi*900/360;
UTS = 1103; %[MPa]
for jj=1:(length(profile1(:,2))-1)
    ncy(jj) = omega*(time_inst(jj+1)-time_inst(jj));
end
ncy = [0; ncy(:)];
Sa_eq = Sa./(1-(Sm./UTS));
% Lumping
ds = 15; %[MPa]
edges = min(Sa_eq):ds:max(Sa_eq); % Bin edges
[bin_idx, ~] = discretize(Sa_eq, edges);  % bin index for each point
nc = zeros(length(edges)-1,1);           % number of cycles in each bin
for i = 1:length(Sa_eq)
    if ~isnan(bin_idx(i))
        nc(bin_idx(i)) = nc(bin_idx(i)) + ncy(i);
    end
end
centers = edges(1:end-1) + diff(edges)/2;
figure(2)
barh(centers, nc, 'FaceColor', [0.7 0.7 0.7]);
xlabel('N cycles');
ylabel('Stress amplitude [MPa]');
title('Stress-Amplitude Histogram (Horizontal)');
grid on;

[Mathieu NOE]

hello

can you also provide the data file ? if it's big , please zip it first.

use the paperclip button to attach it

[Claudio Gerosa]

There it is, thanks

Answer 1

[Mathieu NOE]

hello again

seems to me you can make the code simpler , no need for for loops

I still have a doubt about the omega units ? and therefore also wonder if ncy does indeed represnt the number of revolutions (cycles) per second of the shaft

as the revolution speed is constant and also the time increment , there are some obvious simplifications (like ncy is always the same value )

profile1 = readmatrix('profile1.txt');

samples = size(profile1,1);
time_inst = profile1(:,1);
dt = mean(diff(time_inst)); % sampling time interval

Sa = profile1(:,2);
Sm = profile1(:,3);
omega = 2*pi*900/360; % units ? 
UTS = 1103; %[MPa]

ncy = omega*dt; % number of rotation (cycles) per second ? 

Sa_eq = Sa./(1-(Sm./UTS));

% Lumping
ds = 15; %[MPa]
edges = floor(min(Sa_eq)):ds:ceil(max(Sa_eq)); % Bin edges (rounded values)
nc = histcounts(Sa_eq, edges); %# get counts and bin locations
nc = nc*ncy; % multiply nc by number of rotations per second

centers = edges(1:end-1) + diff(edges)/2;
figure(2)
barh(centers, nc, 'FaceColor', [0.7 0.7 0.7]);
xlabel('N cycles');
ylabel('Stress amplitude [MPa]');
title('Stress-Amplitude Histogram (Horizontal)');
grid on;

●

[Claudio Gerosa]

First of all thanks for your answer. As for the omega it should be 900 rpm in rad/s but thanks to your answer I just realized I've been using the wrong formula all along, silly me. But apart from that it looks to me like we're getting the same results. So overall my code was correct but overly complicated right?

[Claudio Gerosa]

Also ncy should be a vector that contains the cycles of revolution in between two data of the FE analysis

[Mathieu NOE]

so 900 rpm => 900/60 = 15 rotations per seconds (cycles)

but as you get FE results every dt = 0.1562 s , you have ncy = 900/60*0.1562 = 2.3430 cycles between two successive FE results

what is funny is that the wrong omega your are using is indeed numericaly not so far from the correct one :

wrong omega = 2*pi*900/360 = 15.7080

ncy can be a vector but , unless you have a variying shaft speed , or the FE results are given at unevenly spaced time intervals (which is not the case in your example data ) there is no need to compute the same value for each cycle.

you can double check in your code : ncy is indeed a vector , but of same value. there is no variation .

[Claudio Gerosa]

Yeah I realized my mistake while thinking about your answer. I guess it's part of the process ahahah. Thank you very much.

[Mathieu NOE]

as always, my pleasure !