% CHAPTER4 contains the commands used in Chapter 4 of the tutorial CALL: Chapter4 Some of the commands are edited for fast computation. Each set of commands is followed by a 'pause' command. This routine also can print the figures; For printing the figures on directory ../bilder/ edit the file and put printing=1;
Calculates the amplitudes from a cycle count. | |
Calculates the cycle count matrix from a cycle count. | |
Calculates the total Palmgren-Miner damage of a cycle count. | |
Plots a cycle count as a point process in the plane. | |
Calculates the total Palmgren-Miner damage of a cycle matrix. | |
Computes the (Palmgren-Miner) damage matrix from a cycle matrix. | |
Calculates the level crossings from a cycle matrix. | |
Plots a cycle matrix, e.g. a rainflow matrix. | |
The sequence of discretized turning points from a signal. | |
Extracts level-crossing spectrum from data, | |
Extracts turning points from data, | |
Calculates rainflow matrix from discrete turning points. | |
Calculates the fatigue failure time distribution. | |
Calculate transformation, g, proposed by Winterstein | |
Calculates (and plots) a JONSWAP spectral density | |
Simulates process with given irregularity factor and crossing spectrum | |
Calculates discrete levels given the parameter matrix. | |
Calculates the rainflow matrix for a MCTP. | |
Simulates a Markov chain of turning points | |
Makes test matrices for min-max (and max-min) matrices. | |
Plot contents of pdf structures | |
Surface elevation dataset used in WAT version 1.1. | |
Plots SN-data and estimates parameters | |
Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) | |
Simulates a Gaussian process and its derivative from spectrum | |
Estimates the moments of 2'nd order non-linear waves | |
Calculates the number of upcrossings from the turning points. | |
Calculates min2Max and Max2min cycles from a sequence of turning points | |
Finds the rainflow cycles from the sequence of turning points. | |
Prints a caption "made by WAFO" in current figure. | |
Plots data on a Normal distribution paper |
0001 %% CHAPTER4 contains the commands used in Chapter 4 of the tutorial 0002 % 0003 % CALL: Chapter4 0004 % 0005 % Some of the commands are edited for fast computation. 0006 % Each set of commands is followed by a 'pause' command. 0007 % 0008 % This routine also can print the figures; 0009 % For printing the figures on directory ../bilder/ edit the file and put 0010 % printing=1; 0011 0012 % Tested on Matlab 5.3 0013 % History 0014 % Revised pab sept2005 0015 % Added sections -> easier to evaluate using cellmode evaluation. 0016 % revised pab Feb2004 0017 % updated call to lc2sdat 0018 % Created by GL July 13, 2000 0019 % from commands used in Chapter 4 0020 % 0021 0022 %% Chapter 4 Fatigue load analysis and rain-flow cycles 0023 0024 pstate = 'off'; 0025 0026 printing=0; 0027 %set(0,'DefaultAxesFontSize',15') 0028 0029 %% Section 4.3.1 Crossing intensity 0030 load sea.dat; xx_sea=sea; 0031 tp_sea = dat2tp(xx_sea); 0032 lc_sea = tp2lc(tp_sea); 0033 T_sea = xx_sea(end,1)-xx_sea(1,1); 0034 lc_sea(:,2) = lc_sea(:,2)/T_sea; 0035 clf 0036 subplot(221), plot(lc_sea(:,1),lc_sea(:,2)) 0037 title('Crossing intensity, (u, \mu(u))') 0038 subplot(222), semilogx(lc_sea(:,2),lc_sea(:,1)) 0039 title('Crossing intensity, (log \mu(u), u)') 0040 if (printing==1), print -deps ../bilder/fatigue_3.eps 0041 end 0042 wafostamp([],'(ER)') 0043 pause(pstate) 0044 0045 m_sea = mean(xx_sea(:,2)); 0046 f0_sea = interp1(lc_sea(:,1),lc_sea(:,2),m_sea,'linear') 0047 extr_sea = length(tp_sea)/(2*T_sea); 0048 alfa_sea = f0_sea/extr_sea 0049 pause(pstate) 0050 0051 %% Section 4.3.2 Extraction of rainflow cycles 0052 %% Min-max and rainflow cycle plots 0053 RFC_sea=tp2rfc(tp_sea); 0054 mM_sea=tp2mm(tp_sea); 0055 clf 0056 subplot(122), ccplot(mM_sea); 0057 title('min-max cycle count') 0058 subplot(121), ccplot(RFC_sea); 0059 title('Rainflow cycle count') 0060 if (printing==1), print -deps ../bilder/fatigue_4.eps 0061 end 0062 wafostamp([],'(ER)') 0063 pause(pstate) 0064 0065 %% Min-max and rainflow cycle distributions 0066 ampmM_sea=cc2amp(mM_sea); 0067 ampRFC_sea=cc2amp(RFC_sea); 0068 clf 0069 subplot(221), hist(ampmM_sea,25); 0070 title('min-max amplitude distribution') 0071 subplot(222), hist(ampRFC_sea,25); 0072 title('Rainflow amplitude distribution') 0073 if (printing==1), print -deps ../bilder/fatigue_5.eps 0074 end 0075 wafostamp([],'(ER)') 0076 pause(pstate) 0077 0078 %% Section 4.3.3 Simulation of rainflow cycles 0079 %% Simulation of cycles in a Markov model 0080 n=41; param_m=[-1 1 n]; param_D=[1 n n]; 0081 u_markov=levels(param_m); 0082 G_markov=mktestmat(param_m,[-0.2 0.2],0.15,1); 0083 T_markov=5000; 0084 xxD_markov=mctpsim({G_markov []},T_markov); 0085 xx_markov=[(1:T_markov)' u_markov(xxD_markov)']; 0086 clf 0087 plot(xx_markov(1:50,1),xx_markov(1:50,2)) 0088 title('Markov chain of turning points') 0089 if (printing==1), print -deps ../bilder/fatigue_6.eps 0090 end 0091 wafostamp([],'(ER)') 0092 pause(pstate) 0093 0094 0095 %% Rainflow cycles in a transformed Gaussian model 0096 %% Hermite transformed wave data and rainflow filtered turning points, h = 0.2. 0097 me = mean(xx_sea(:,2)); 0098 sa = std(xx_sea(:,2)); 0099 Hm0_sea = 4*sa; 0100 Tp_sea = 1/max(lc_sea(:,2)); 0101 spec = jonswap([],[Hm0_sea Tp_sea]); 0102 0103 [sk, ku] = spec2skew(spec); 0104 spec.tr = hermitetr([],[sa sk ku me]); 0105 param_h = [-1.5 2 51]; 0106 spec_norm = spec; 0107 spec_norm.S = spec_norm.S/sa^2; 0108 xx_herm = spec2sdat(spec_norm,[2^15 1],0.1); 0109 % ????? PJ, JR 11-Apr-2001 0110 % NOTE, in the simulation program spec2sdat 0111 %the spectrum must be normalized to variance 1 0112 % ????? 0113 h = 0.2; 0114 [dtp,u_herm,xx_herm_1]=dat2dtp(param_h,xx_herm,h); 0115 clf 0116 plot(xx_herm(:,1),xx_herm(:,2),'k','LineWidth',2); hold on; 0117 plot(xx_herm_1(:,1),xx_herm_1(:,2),'k--','Linewidth',2); 0118 axis([0 50 -1 1]), hold off; 0119 title('Rainflow filtered wave data') 0120 if (printing==1), print -deps ../bilder/fatigue_7.eps 0121 end 0122 wafostamp([],'(ER)') 0123 pause(pstate) 0124 0125 %% Rainflow cycles and rainflow filtered rainflow cycles in the transformed Gaussian process. 0126 tp_herm=dat2tp(xx_herm); 0127 RFC_herm=tp2rfc(tp_herm); 0128 mM_herm=tp2mm(tp_herm); 0129 h=0.2; 0130 [dtp,u,tp_herm_1]=dat2dtp(param_h,xx_herm,h); 0131 RFC_herm_1 = tp2rfc(tp_herm_1); 0132 clf 0133 subplot(121), ccplot(RFC_herm) 0134 title('h=0') 0135 subplot(122), ccplot(RFC_herm_1) 0136 title('h=0.2') 0137 if (printing==1), print -deps ../bilder/fatigue_8.eps 0138 end 0139 wafostamp([],'(ER)') 0140 pause(pstate) 0141 0142 %% Section 4.3.4 Calculating the rainflow matrix 0143 0144 0145 Grfc_markov=mctp2rfm({G_markov []}); 0146 clf 0147 subplot(121), cmatplot(u_markov,u_markov,G_markov), axis('square') 0148 subplot(122), cmatplot(u_markov,u_markov,Grfc_markov), axis('square') 0149 wafostamp([],'(ER)') 0150 pause(pstate) 0151 0152 %% 0153 clf 0154 cmatplot(u_markov,u_markov,{G_markov Grfc_markov},3) 0155 wafostamp([],'(ER)') 0156 pause(pstate) 0157 0158 %% Min-max-matrix and theoretical rainflow matrix for test Markov sequence. 0159 cmatplot(u_markov,u_markov,{G_markov Grfc_markov},4) 0160 subplot(121), axis('square'), title('min2max transition matrix') 0161 subplot(122), axis('square'), title('Rainflow matrix') 0162 if (printing==1), print -deps ../bilder/fatigue_9.eps 0163 end 0164 wafostamp([],'(ER)') 0165 pause(pstate) 0166 0167 %% Observed and theoretical rainflow matrix for test Markov sequence. 0168 n=length(u_markov); 0169 Frfc_markov=dtp2rfm(xxD_markov,n); 0170 clf 0171 cmatplot(u_markov,u_markov,{Frfc_markov Grfc_markov*T_markov/2},3) 0172 subplot(121), axis('square'), title('Observed rainflow matrix') 0173 subplot(122), axis('square'), title('Theoretical rainflow matrix') 0174 if (printing==1), print -deps ../bilder/fatigue_10.eps 0175 end 0176 wafostamp([],'(ER)') 0177 pause(pstate) 0178 0179 %% Smoothed observed and calculated rainflow matrix for test Markov sequence. 0180 tp_markov=dat2tp(xx_markov); 0181 RFC_markov=tp2rfc(tp_markov); 0182 h=1; 0183 Frfc_markov_smooth=cc2cmat(param_m,RFC_markov,[],1,h); 0184 clf 0185 cmatplot(u_markov,u_markov,{Frfc_markov_smooth Grfc_markov*T_markov/2},4) 0186 subplot(121), axis('square'), title('Smoothed observed rainflow matrix') 0187 subplot(122), axis('square'), title('Theoretical rainflow matrix') 0188 if (printing==1), print -deps ../bilder/fatigue_11.eps 0189 end 0190 wafostamp([],'(ER)') 0191 pause(pstate) 0192 0193 %% Rainflow matrix from spectrum 0194 clf 0195 %GmM3_herm=spec2mmtpdf(spec,[],'Mm',[],[],2); 0196 GmM3_herm=spec2cmat(spec,[],'Mm',[],param_h,2); 0197 pdfplot(GmM3_herm) 0198 wafostamp([],'(ER)') 0199 pause(pstate) 0200 0201 0202 %% Min-max matrix and theoretical rainflow matrix for Hermite-transformed Gaussian waves. 0203 Grfc_herm=mctp2rfm({GmM3_herm.f []}); 0204 u_herm=levels(param_h); 0205 clf 0206 cmatplot(u_herm,u_herm,{GmM3_herm.f Grfc_herm},4) 0207 subplot(121), axis('square'), title('min-max matrix') 0208 subplot(122), axis('square'), title('Theoretical rainflow matrix') 0209 if (printing==1), print -deps ../bilder/fatigue_12.eps 0210 end 0211 wafostamp([],'(ER)') 0212 pause(pstate) 0213 0214 %% 0215 clf 0216 Grfc_direct_herm=spec2cmat(spec,[],'rfc',[],[],2); 0217 subplot(121), pdfplot(GmM3_herm), axis('square'), hold on 0218 subplot(122), pdfplot(Grfc_direct_herm), axis('square'), hold off 0219 if (printing==1), print -deps ../bilder/fig_mmrfcjfr.eps 0220 end 0221 wafostamp([],'(ER)') 0222 pause(pstate) 0223 0224 0225 %% Observed smoothed and theoretical min-max matrix, 0226 %% (and observed smoothed and theoretical rainflow matrix for Hermite-transformed Gaussian waves). 0227 tp_herm=dat2tp(xx_herm); 0228 RFC_herm=tp2rfc(tp_herm); 0229 mM_herm=tp2mm(tp_herm); 0230 h=0.2; 0231 FmM_herm_smooth=cc2cmat(param_h,mM_herm,[],1,h); 0232 Frfc_herm_smooth=cc2cmat(param_h,RFC_herm,[],1,h); 0233 T_herm=xx_herm(end,1)-xx_herm(1,1); 0234 clf 0235 cmatplot(u_herm,u_herm,{FmM_herm_smooth GmM3_herm.f*length(mM_herm) ; ... 0236 Frfc_herm_smooth Grfc_herm*length(RFC_herm)},4) 0237 subplot(221), axis('square'), title('Observed smoothed min-max matrix') 0238 subplot(222), axis('square'), title('Theoretical min-max matrix') 0239 subplot(223), axis('square'), title('Observed smoothed rainflow matrix') 0240 subplot(224), axis('square'), title('Theoretical rainflow matrix') 0241 if (printing==1), print -deps ../bilder/fatigue_13.eps 0242 end 0243 wafostamp([],'(ER)') 0244 pause(pstate) 0245 0246 %% Section 4.3.5 Simulation from crossings and rainflow structure 0247 0248 %% Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) 0249 %% for simulated process with irregularity factor 0.25. 0250 clf 0251 cross_herm=dat2lc(xx_herm); 0252 alpha1=0.25; 0253 alpha2=0.75; 0254 xx_herm_sim1=lc2sdat(cross_herm,500,alpha1); 0255 cross_herm_sim1=dat2lc(xx_herm_sim1); 0256 subplot(211) 0257 plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2))) 0258 hold on 0259 stairs(cross_herm_sim1(:,1),... 0260 cross_herm_sim1(:,2)/max(cross_herm_sim1(:,2))) 0261 hold off 0262 title('Crossing intensity, \alpha = 0.25') 0263 subplot(212) 0264 plot(xx_herm_sim1(:,1),xx_herm_sim1(:,2)) 0265 title('Simulated load, \alpha = 0.25') 0266 if (printing==1), print -deps ../bilder/fatigue_14_25.eps 0267 end 0268 wafostamp([],'(ER)') 0269 0270 %% Crossing spectrum (smooth curve) and obtained spectrum (wiggled curve) 0271 %% for simulated process with irregularity factor 0.75. 0272 xx_herm_sim2=lc2sdat(cross_herm,500,alpha2); 0273 cross_herm_sim2=dat2lc(xx_herm_sim2); 0274 subplot(211) 0275 plot(cross_herm(:,1),cross_herm(:,2)/max(cross_herm(:,2))) 0276 hold on 0277 stairs(cross_herm_sim2(:,1),... 0278 cross_herm_sim2(:,2)/max(cross_herm_sim2(:,2))) 0279 hold off 0280 title('Crossing intensity, \alpha = 0.75') 0281 subplot(212) 0282 plot(xx_herm_sim2(:,1),xx_herm_sim2(:,2)) 0283 title('Simulated load, \alpha = 0.75') 0284 if (printing==1), print -deps ../bilder/fatigue_14_75.eps 0285 end 0286 wafostamp([],'(ER)') 0287 pause(pstate) 0288 0289 %% Section 4.4 Fatigue damage and fatigue life distribution 0290 %% Section 4.4.1 Introduction 0291 beta=3.2; gam=5.5E-10; T_sea=xx_sea(end,1)-xx_sea(1,1); 0292 d_beta=cc2dam(RFC_sea,beta)/T_sea; 0293 time_fail=1/gam/d_beta/3600 %in hours of the specific storm 0294 pause(pstate) 0295 0296 %% Section 4.4.2 Level crossings 0297 %% Crossing intensity as calculated from the Markov matrix (solid curve) and from the observed rainflow matrix (dashed curve). 0298 clf 0299 mu_markov=cmat2lc(param_m,Grfc_markov); 0300 muObs_markov=cmat2lc(param_m,Frfc_markov/(T_markov/2)); 0301 clf 0302 plot(mu_markov(:,1),mu_markov(:,2),muObs_markov(:,1),muObs_markov(:,2),'--') 0303 title('Theoretical and observed crossing intensity ') 0304 if (printing==1), print -deps ../bilder/fatigue_15.eps 0305 end 0306 wafostamp([],'(ER)') 0307 pause(pstate) 0308 0309 %% Section 4.4.3 Damage 0310 %% Distribution of damage from different RFC cycles, from calculated theoretical and from observed rainflow matrix. 0311 beta = 4; 0312 Dam_markov = cmat2dam(param_m,Grfc_markov,beta) 0313 DamObs1_markov = cc2dam(RFC_markov,beta)/(T_markov/2) 0314 DamObs2_markov = cmat2dam(param_m,Frfc_markov,beta)/(T_markov/2) 0315 pause(pstate) 0316 0317 Dmat_markov = cmat2dmat(param_m,Grfc_markov,beta); 0318 DmatObs_markov = cmat2dmat(param_m,Frfc_markov,beta)/(T_markov/2); 0319 clf 0320 subplot(121), cmatplot(u_markov,u_markov,Dmat_markov,4) 0321 title('Theoretical damage matrix') 0322 subplot(122), cmatplot(u_markov,u_markov,DmatObs_markov,4) 0323 title('Observed damage matrix') 0324 if (printing==1), print -deps ../bilder/fatigue_16.eps 0325 end 0326 wafostamp([],'(ER)') 0327 pause(pstate) 0328 0329 0330 %% 0331 %Damplus_markov = lc2dplus(mu_markov,beta) 0332 pause(pstate) 0333 0334 %% Section 4.4.4 Estimation of S-N curve 0335 0336 %% Load SN-data and plot in log-log scale. 0337 SN = load('sn.dat'); 0338 s = SN(:,1); 0339 N = SN(:,2); 0340 clf 0341 loglog(N,s,'o'), axis([0 14e5 10 30]) 0342 %if (printing==1), print -deps ../bilder/fatigue_?.eps end 0343 wafostamp([],'(ER)') 0344 pause(pstate) 0345 0346 0347 %% Check of S-N-model on normal probability paper. 0348 0349 wnormplot(reshape(log(N),8,5)) 0350 if (printing==1), print -deps ../bilder/fatigue_17.eps 0351 end 0352 wafostamp([],'(ER)') 0353 pause(pstate) 0354 0355 %% Estimation of S-N-model on linear scale. 0356 0357 [e0,beta0,s20] = snplot(s,N,12) 0358 title('S-N-data with estimated N(s)','FontSize',20) 0359 set(gca,'FontSize',20) 0360 if (printing==1), print -deps ../bilder/fatigue_18a.eps 0361 end 0362 wafostamp([],'(ER)') 0363 pause(pstate) 0364 0365 %% Estimation of S-N-model on log-log scale. 0366 [e0,beta0,s20] = snplot(s,N,14) 0367 title('S-N-data with estimated N(s)','FontSize',20) 0368 set(gca,'FontSize',20) 0369 if (printing==1), print -deps ../bilder/fatigue_18b.eps 0370 end 0371 wafostamp([],'(ER)') 0372 pause(pstate) 0373 0374 %% Section 4.4.5 From S-N curve to fatigue life distribution 0375 %% Damage intensity as function of $\beta$ 0376 beta = 3:0.1:8; 0377 DRFC = cc2dam(RFC_sea,beta); 0378 dRFC = DRFC/T_sea; 0379 plot(beta,dRFC), axis([3 8 0 0.25]) 0380 title('Damage intensity as function of \beta') 0381 if (printing==1), print -deps ../bilder/fatigue_19.eps 0382 end 0383 wafostamp([],'(ER)') 0384 pause(pstate) 0385 0386 %% Fatigue life distribution with sea load. 0387 dam0 = cc2dam(RFC_sea,beta0)/T_sea; 0388 [t0,F0] = ftf(e0,dam0,s20,0.5,1); 0389 [t1,F1] = ftf(e0,dam0,s20,0,1); 0390 [t2,F2] = ftf(e0,dam0,s20,5,1); 0391 plot(t0,F0,t1,F1,t2,F2) 0392 title('Fatigue life distribution function') 0393 if (printing==1), print -deps ../bilder/fatigue_20.eps 0394 end 0395 wafostamp([],'(ER)') 0396 save Chap4 0397 0398
Comments or corrections to the WAFO group