Home > wafo > wdemos > itmkurs > itmkurs_lab3.m

# itmkurs_lab3

## PURPOSE

Script to computer exercises 3

## SYNOPSIS

This is a script file.

## DESCRIPTION

``` ITMKURS_LAB3 Script to computer exercises 3

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Power Spectrum and Rainflow Analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%```

## CROSS-REFERENCE INFORMATION

This function calls:
 cc2dam Calculates the total Palmgren-Miner damage of a cycle count. cmat2dam Calculates the total Palmgren-Miner damage of a cycle matrix. createspec Spectrum structure constructor dat2spec Estimate one-sided spectral density from data. dat2tp Extracts turning points from data, iter Calculates a Markov matrix given a rainflow matrix jonswap Calculates (and plots) a JONSWAP spectral density levels Calculates discrete levels given the parameter matrix. mctpsim Simulates a Markov chain of turning points oscspec Spectral density for a harmonic oscillator pdfplot Plot contents of pdf structures spec2cmat Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) spec2dplus Calculates an upper bound of the damage intensity explicitly spec2mom Calculates spectral moments from spectrum spec2sdat Simulates a Gaussian process and its derivative from spectrum specinterp Interpolation and zero-padding of spectrum torsethaugen Calculates a double peaked (swell + wind) spectrum tp2rfc Finds the rainflow cycles from the sequence of turning points. wspecplot Plot a spectral density
This function is called by:

## SOURCE CODE

```001 %ITMKURS_LAB3 Script to computer exercises 3
002 %
003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
004 % Power Spectrum and Rainflow Analysis
005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
006
007 %
008 %%Power Spectrum and Rainflow Analysis
009 %%Power Spectrum
010 %
011
012 S=jonswap;
013 S=oscspec;
014 S=torsethaugen;
016
017 %Solution Problem1:
018
019 L=spec2mom(S,4);
020 f0=sqrt(L(2)/L(1))/2/pi
021 s=sqrt(L(1))
022 alfa=f0/(sqrt(L(3)/L(2))/2/pi)
023 wspecplot(S);
024 2*pi*f0
025
026 beta=3:0.1:5;
027 gam=5E-9;
028 dpl=gam*spec2dplus(S,beta);
029 Tfpl=(1./dpl)/3600/24  %in days
030 clf, plot(beta,Tfpl)
031
032 %
033 %%Simulation of \$X\$ and the damage intensity
034 %
035
036
037 T=1000/f0 % in seconds
038 gam=5E-9;
039 100*T*gam*spec2dplus(S,4.22) % in percent
040
041
042 max(S.w)/pi
043 dt=1./f0/60
044 S1=specinterp(S,dt);
045 [max(S1.w) pi/dt]
046 clf, wspecplot(S), hold on, wspecplot(S1,1,'r.')
047 clf, plot(beta,Tfpl)
048 hold on
049 for  i=1:5
050    X=spec2sdat(S1,60000);
051    tp=dat2tp(X);
052    rfc=tp2rfc(tp);
053    db=cc2dam(rfc,beta);
054    plot(beta,(1000/f0)*(1./db)*(1/gam)/3600/24,'r.') % Tf in days
055 end
056
057 %
058 %%Theoretical computation of damage intensity
059 %
060
061
062 S=createspec;
063 w=levels([0 4 253]);
064 S.S=2*exp(-2*abs(w-2).^1.4);
065 S.S=S.S+9*exp(-8*(w-0.5).^2);
066 S.w=w;
067 S.note=['S(w)=9*exp(-8(w-0.5)^2)+2*exp(-2|w-2|.^1.4)']
068 wspecplot(S)
069
070 L=spec2mom(S,4);
071 f0=sqrt(L(2)/L(1))/2/pi
072 dpl=gam*spec2dplus(S,beta);
073 Tfpl=(1./dpl)/3600/24  %in days
074 dt=1./f0/60
075 S1=specinterp(S,dt);
076 figure(1)
077 clf, plot(beta,Tfpl)
078 hold on
079 for  i=1:10
080   X=spec2sdat(S1,60000);
081   tp=dat2tp(X);
082   rfc=tp2rfc(tp);
083   db=cc2dam(rfc,beta);
084   plot(beta,(1000/f0)*(1./db)*(1/gam)/3600/24,'r.') % Tf in days
085 end
086
087
088
089 a=sqrt(L(3)/L(2))/2/pi
090 s=sqrt(L(1))
091 help spec2cmat
092 paramu=[-4.5*s 4.5*s 46];
093 nit=2;
094 figure(2), clf
095 [frfc fMm]=spec2cmat(S,[],'rfc',[],paramu,nit);
096 hold on, plot(rfc(:,2),rfc(:,1),'.')
097 clf, pdfplot(fMm)
098 dg=a*cmat2dam(paramu,frfc.f,beta);
099 figure(1), plot(beta,(1./dg)*(1/gam)/3600/24,'g') % Tf in days
100
101
102
103 figure(2), clf
104 [frfc1 fMm1]=spec2cmat(S,[],'rfc',[],paramu,3);
105 dg1=a*cmat2dam(paramu,frfc1.f,beta);
106 figure(1)
107 plot(beta,(1./dg1)*(1/gam)/3600/24,'k') % Tf in days
108
109
110 help mctpsim
111 F{1,1}=fMm.f;
112 F{1,2}=fMm.f';
113 mctp=mctpsim(F,1000);
114 clf,plot(mctp(1:40))
115
116 frfc1=zeros(size(frfc.f));
117 frfc1=triu(frfc.f,5);
118 fMm1=iter(frfc1,fMm.f,20,0.001);
119 clf,contour(fMm.f,20)
120 hold on,  contour(fMm1,20)
121 F{1,1}=fMm1;
122 F{1,2}=fMm1';
123 mctp1=mctpsim(F,1000);
124 clf, subplot(2,1,1)
125 plot(mctp(1:40))
126 subplot(2,1,2)
127 plot(mctp1(1:40))
128
129
130
131
132
133
134
135
136```

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