Home > wafo > papers > tutorcom > Chapter4.m

Chapter4

PURPOSE ^

% CHAPTER4 contains the commands used in Chapter 4 of the tutorial

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 % 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;

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003