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

조회 수: 30 (최근 30일)
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!

채택된 답변

VBBV
VBBV 2024년 3월 28일
편집: VBBV 2024년 3월 28일
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
Tyl
Tyl 2024년 3월 29일
편집: Tyl 2024년 3월 29일
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.

댓글을 달려면 로그인하십시오.

추가 답변 (2개)

Strider
Strider 2024년 3월 28일
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 2024년 4월 19일

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.

카테고리

Help CenterFile Exchange에서 Vibration Analysis에 대해 자세히 알아보기

제품


릴리스

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by