TORSETHAUGEN Calculates a double peaked (swell + wind) spectrum CALL: [S, Ss, Sw] = torsethaugen(w,data,plotflag); [S, Ss, Sw] = torsethaugen(wc,data,plotflag); S, Ss, Sw = spectrum struct for total, swell and wind, respectively. w = angular frequency (default linspace(0,wc,257)) wc = angular cutoff frequency (default 33/Tp) data = [Hm0 Tp A] Hm0 = Significant wave height (default 7) Tp = 2*pi/wp, primary peak period (deafult 11) A = alpha, normalization factor, (default -1) A<0 : A calculated by integration so that int S dw =Hm0^2/16 A==0 : A = (1+f1(N,M)*log(gammai)^f2(N,M))/gammai, original parametric normalization plotflag = 0, do not plot the spectrum (default). 1, plot the spectrum. For zero values, NaN's or parameters not specified in DATA the default values are used. Model: S(w)=Ss(w)+Sw(w) where Ss and Sw are modified JONSWAP spectrums for swell and wind peak, respectively. The energy is divided between the two peaks according to empirical parameters, which peak that is primary depends on parameters. The empirical parameters are found for classes of Hm0 and Tp, originating from a dataset consisting of 20 000 spectra divided into 146 different classes of Hm0 and Tp. (Data measured at the Statfjord field in the North Sea in a period from 1980 to 1989.) The range of the measured Hm0 and Tp for the dataset are from 0.5 to 11 meters and from 3.5 to 19 sec, respectively. See Torsethaugen (1996). Preliminary comparisons with spectra from other areas indicate that some of the empirical parameters are dependent on geographical location. Thus the model must be used with care for other areas than the North Sea and sea states outside the area where measured data are available. Example: [S,Sw,Ss]=torsethaugen([],[6 8],1) See also jonswap, pmspec.
Spectrum structure constructor | |
returns the constant acceleration of gravity | |
Numerical integration with the Simpson method | |
Plot a spectral density | |
Display message and abort function. | |
Gamma function. | |
Hold current graph. | |
Return hold state. | |
Linearly spaced vector. | |
Convert number to string. (Fast version) | |
Concatenate strings. |
% CHAPTER1 demonstrates some applications of WAFO | |
% CHAPTER2 Modelling random loads and stochastic waves | |
% CHAPTER3 Demonstrates distributions of wave characteristics | |
Script to computer exercises 3 | |
Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. | |
setup all global variables of the WAFODEMO |
0001 function [S, Sw, Ss]=torsethaugen(w1,sdata,plotflag) 0002 % TORSETHAUGEN Calculates a double peaked (swell + wind) spectrum 0003 % 0004 % CALL: [S, Ss, Sw] = torsethaugen(w,data,plotflag); 0005 % [S, Ss, Sw] = torsethaugen(wc,data,plotflag); 0006 % 0007 % S, Ss, Sw = spectrum struct for total, swell and wind, respectively. 0008 % w = angular frequency (default linspace(0,wc,257)) 0009 % wc = angular cutoff frequency (default 33/Tp) 0010 % data = [Hm0 Tp A] 0011 % Hm0 = Significant wave height (default 7) 0012 % Tp = 2*pi/wp, primary peak period (deafult 11) 0013 % A = alpha, normalization factor, (default -1) 0014 % A<0 : A calculated by integration so that int S dw =Hm0^2/16 0015 % A==0 : A = (1+f1(N,M)*log(gammai)^f2(N,M))/gammai, original 0016 % parametric normalization 0017 % plotflag = 0, do not plot the spectrum (default). 0018 % 1, plot the spectrum. 0019 % 0020 % For zero values, NaN's or parameters not specified in DATA the 0021 % default values are used. 0022 % 0023 % Model: S(w)=Ss(w)+Sw(w) 0024 % where Ss and Sw are modified JONSWAP spectrums for 0025 % swell and wind peak, respectively. 0026 % The energy is divided between the two peaks according 0027 % to empirical parameters, which peak that is primary depends on parameters. 0028 % The empirical parameters are found for classes of Hm0 and Tp, 0029 % originating from a dataset consisting of 20 000 spectra divided 0030 % into 146 different classes of Hm0 and Tp. (Data measured at the 0031 % Statfjord field in the North Sea in a period from 1980 to 1989.) 0032 % The range of the measured Hm0 and Tp for the dataset 0033 % are from 0.5 to 11 meters and from 3.5 to 19 sec, respectively. 0034 % See Torsethaugen (1996). 0035 % 0036 % Preliminary comparisons with spectra from other areas indicate that 0037 % some of the empirical parameters are dependent on geographical location. 0038 % Thus the model must be used with care for other areas than the 0039 % North Sea and sea states outside the area where measured data 0040 % are available. 0041 % 0042 % Example: [S,Sw,Ss]=torsethaugen([],[6 8],1) 0043 % 0044 % See also jonswap, pmspec. 0045 0046 % References 0047 % Torsethaugen, K. (1996) 0048 % Model for a doubly peaked wave spectrum 0049 % Report No. STF22 A96204. SINTEF Civil and Environm. Engineering, Trondheim 0050 % 0051 % Torsethaugen, K. (1994) 0052 % 'Model for a doubly peaked spectrum. Lifetime and fatigue strength 0053 % estimation implications.' 0054 % International Workshop on Floating Structures in Coastal zone, 0055 % Hiroshima, November 1994. 0056 % 0057 % Torsethaugen, K. (1993) 0058 % 'A two peak wave spectral model.' 0059 % In proceedings OMAE, Glasgow 0060 0061 0062 % Tested on Matlab 5.3 0063 % History: 0064 % Revised pab april 2005 0065 % revised pab 12.12.2000 0066 % - fixed a bug: Hpw = 0 or Hps = 0 is now handled properly 0067 % - sdata may now contain NaN's 0068 % revised pab 16.02.2000 0069 % -added ih, see at the end 0070 % revised pab 20.12.1999 0071 % -added norm 0072 % -changed default value of Hm0 to 7 (equal to jonswap) 0073 % revised pab 01.12.1999 0074 % - updated the reference 0075 % revised es 23.09.1999 0076 % - updated documentation 0077 % - added note 0078 % by pab 23.06.1999 0079 0080 monitor=0; 0081 % default values 0082 Hm0 = 7; 0083 Tp = 11; 0084 Aj = -1; 0085 data2 = [Hm0, Tp, Aj]; 0086 nd2 = length(data2); 0087 if nargin<3|isempty(plotflag), plotflag = 0; end 0088 if (nargin>1) & ~isempty(sdata), 0089 nd = length(sdata); 0090 ind = find(~isnan(sdata(1:min(nd,nd2)))); 0091 if any(ind), % replace default values with those from input data 0092 data2(ind) = sdata(ind); 0093 end 0094 end 0095 if (nd2>0) & (data2(1)>0), Hm0 = data2(1);end 0096 if (nd2>1) & (data2(2)>0), Tp = data2(2);end 0097 if (nd2>2) & (data2(3)==0), Aj = data2(3);end 0098 0099 w = []; 0100 if nargin<1|isempty(w1), 0101 wc = 33/Tp; 0102 elseif length(w1)==1, 0103 wc = w1; 0104 else 0105 w = w1 ; 0106 end 0107 nw = 257; 0108 0109 if isempty(w), w = linspace(0,wc,nw).'; end 0110 0111 w=w(:); 0112 0113 n = length(w); 0114 S = createspec('freq'); 0115 S.S = zeros(n,1); 0116 S.w = w; 0117 S.norm = 0; % not normalized spectra 0118 Sw = S; 0119 Ss = S; 0120 0121 if Hm0>11| Hm0>max((Tp/3.6).^2, (Tp-2)*12/11) 0122 disp('Warning: Hm0 is outside the valid range') 0123 disp('The validity of the spectral density is questionable') 0124 end 0125 if Tp>20|Tp<3 0126 disp('Warning: Tp is outside the valid range') 0127 disp('The validity of the spectral density is questionable') 0128 end 0129 wp = 2*pi/Tp; 0130 0131 0132 %initializing 0133 0134 %sigma used in the peak enhancement factor 0135 sa = 0.07; % sigma_a for w/wp<1 0136 sb = 0.09; % sigma_b for w/wp>1 0137 0138 g = gravity; % m/s^2 0139 0140 % The parameter values below are found comparing the 0141 % model to average measured spectra for the Statfjord Field 0142 % in the Northern North Sea. 0143 Af = 6.6; %m^(-1/3)*sec 0144 AL = 2; %sec/sqrt(m) 0145 Au = 25; %sec 0146 KG = 35; 0147 KG0 = 3.5; 0148 KG1 = 1; % m 0149 r = 0.857; % 6/7 0150 K0 = 0.5; %1/sqrt(m) 0151 K00 = 3.2; 0152 0153 M0 = 4; 0154 B1 = 2; %sec 0155 B2 = 0.7; 0156 B3 = 3.0; %m 0157 S0 = 0.08; %m^2*s 0158 S1 = 3; %m 0159 0160 % Preliminary comparisons with spectra from other areas indicate that 0161 % the parameters on the line below can be dependent on geographical location 0162 A10 = 0.7; A1 = 0.5; A20 = 0.6; A2 = 0.3; A3 = 6; 0163 0164 Tf = Af*(Hm0)^(1/3); 0165 Tl = AL*sqrt(Hm0); % lower limit 0166 Tu = Au; % upper limit 0167 0168 %Non-dimensional scales 0169 % New call pab April 2005 0170 El = min(max((Tf-Tp)/(Tf-Tl),0),1); %wind sea 0171 Eu = min(max((Tp-Tf)/(Tu-Tf),0),1); %Swell 0172 0173 %El = ((Tf-Tp)/(Tf-Tl)), %wind sea 0174 %Eu = ((Tp-Tf)/(Tu-Tf)), %Swell 0175 0176 0177 if Tp<Tf, % Wind dominated seas 0178 % Primary peak (wind dominated) 0179 Nw = K0*sqrt(Hm0)+K00; % high frequency exponent 0180 Mw = M0; % spectral width exponent 0181 Rpw = min((1-A10)*exp(-(El/A1)^2)+A10,1); 0182 Hpw = Rpw*Hm0; % significant waveheight wind 0183 Tpw = Tp; % primary peak period 0184 % peak enhancement factor 0185 gammaw = KG*(1+KG0*exp(-Hm0/KG1))*(2*pi/g*Rpw*Hm0/(Tp^2))^r; 0186 gammaw = max(gammaw,1); 0187 % Secondary peak (swell) 0188 Ns = Nw; % high frequency exponent 0189 Ms = Mw; % spectral width exponent 0190 Rps = sqrt(1-Rpw^2); 0191 Hps = Rps*Hm0; % significant waveheight swell 0192 Tps = Tf+B1; 0193 gammas = 1; 0194 0195 0196 if Rps > 0.1 0197 disp(' Spectrum for Wind dominated sea') 0198 else 0199 disp(' Spectrum for pure wind sea') 0200 end 0201 else %swell dominated seas 0202 0203 % Primary peak (swell) 0204 Ns = K0*sqrt(Hm0)+K00; %high frequency exponent 0205 Ms = M0; %spectral width exponent 0206 Rps = min((1-A20)*exp(-(Eu/A2)^2)+A20,1); 0207 Hps = Rps*Hm0; % significant waveheight swell 0208 Tps = Tp; % primary peak period 0209 % peak enhancement factor 0210 gammas = KG*(1+KG0*exp(-Hm0/KG1))*(2*pi/g*Hm0/(Tf^2))^r*(1+A3*Eu); 0211 gammas = max(gammas,1); 0212 0213 % Secondary peak (wind) 0214 Nw = Ns; % high frequency exponent 0215 Mw = M0*(1-B2*exp(-Hm0/B3)); % spectral width exponent 0216 Rpw = sqrt(1-Rps^2); 0217 Hpw = Rpw*Hm0; % significant waveheight wind 0218 G0w = Mw/((Nw/Mw)^(-(Nw-1)/Mw)*gamma((Nw-1)/Mw )); %normalizing factor 0219 if Hpw>0, 0220 Tpw = (16*S0*(1-exp(-Hm0/S1))*(0.4)^Nw/( G0w*Hpw^2))^(-1/(Nw-1)); 0221 else 0222 Tpw = inf; 0223 end 0224 %Tpw = max(Tpw,2.5) 0225 gammaw = 1; 0226 0227 if Rpw > 0.1 0228 disp(' Spectrum for swell dominated sea') 0229 else 0230 disp(' Spectrum for pure swell sea') 0231 end 0232 end 0233 if (3.6*sqrt(Hm0)<= Tp & Tp<=5*sqrt(Hm0)) 0234 disp(' Jonswap range') 0235 end 0236 if monitor, 0237 disp(['Hm0 = ' num2str(Hm0)]) 0238 disp(['Ns, Ms = ' num2str([Ns Ms]) ' Nw, Mw = ' num2str([Nw Mw])]) 0239 disp(['gammas = ' num2str(gammas) ' gammaw = ' num2str(gammaw)]) 0240 disp(['Rps = ' num2str(Rps) ' Rpw = ' num2str(Rpw)]) 0241 disp(['Hps = ' num2str(Hps) ' Hpw = ' num2str(Hpw)]) 0242 disp(['Tps = ' num2str(Tps) ' Tpw = ' num2str(Tpw)]) 0243 end 0244 0245 %G0s=Ms/((Ns/Ms)^(-(Ns-1)/Ms)*gamma((Ns-1)/Ms )); %normalizing factor 0246 0247 % Wind part 0248 0249 Sw = localJonswap(w,Hpw,Tpw,gammaw,Nw,Mw,Aj); 0250 % Swell part 0251 Ss = localJonswap(w,Hps,Tps,gammas,Ns,Ms,Aj); 0252 0253 0254 0255 S.S = Ss.S+Sw.S; % total spectrum 0256 S.note=strcat('Torsethaugen, Hm0=',num2str(Hm0),', Tp=',num2str(Tp)); 0257 if max(Ss.S)<max(Sw.S) 0258 Ss.note=strcat('Swell, secondary peak of Torsethaugen, Hm0=',num2str(Hm0),', Tp=',num2str(Tp)); 0259 Sw.note=strcat('Wind, primary peak of Torsethaugen, Hm0=',num2str(Hm0),', Tp=',num2str(Tp)); 0260 else 0261 Ss.note=strcat('Swell, primary peak of Torsethaugen, Hm0=',num2str(Hm0),', Tp=',num2str(Tp)); 0262 Sw.note=strcat('Wind, secondary peak of Torsethaugen, Hm0=',num2str(Hm0),', Tp=',num2str(Tp)); 0263 end 0264 if plotflag 0265 wspecplot(S,plotflag,'b') 0266 ih = ishold; 0267 hold on 0268 wspecplot(Ss,plotflag,'b--') 0269 wspecplot(Sw,plotflag,'b--') 0270 if ~ih, hold off, end 0271 end 0272 0273 function S = localJonswap(w,Hm0,Tp,gammai,N,M,A) 0274 %local JONSWAP formulation 0275 0276 n1 = length(w); 0277 S = createspec('freq'); 0278 S.S = zeros(n1,1); 0279 S.w = w; 0280 0281 if (Hm0<=0) 0282 return 0283 end 0284 0285 %sigma used in the peak enhancement factor 0286 sa = 0.07; % sigma_a for wn<1 0287 sb = 0.09; % sigma_b for wn>1 0288 0289 wp = 2*pi/Tp; 0290 wn = w/wp; 0291 0292 B = (N-1)/M; 0293 G0 = (N/M).^(B)*(M/gamma(B)); % Normalizing factor related to Pierson-Moskovitz form 0294 G1 = (Hm0/4)^2/wp*G0; 0295 0296 % for w>wp 0297 k = find(wn>1); 0298 if any(k) 0299 S.S(k)=G1./(wn(k).^N).*(gammai.^(exp(-(wn(k)-1).^2 ... 0300 /(2*sb^2)))).*exp(-N/M*(wn(k)).^(-M)); 0301 end 0302 % for 0<w<=wp 0303 k = find(0<wn & wn<=1); 0304 if any(k) 0305 S.S(k)=G1./(wn(k).^N).*(gammai.^(exp(-(wn(k)-1).^2 ... 0306 /(2*sa^2)))).*exp(-N/M*(wn(k)).^(-M)); 0307 end 0308 0309 if (gammai==1) 0310 A = 1; 0311 end 0312 if A<0, % normalizing by integration 0313 area = simpson(w,S.S); 0314 if area>0 0315 A = (Hm0/4)^2./area;% make sure m0=Hm0^2/16=int S(w)dw 0316 else 0317 A = 0; 0318 end 0319 elseif A==0,% original normalization 0320 % NOTE: that Hm0^2/16 generally is not equal to intS(w)dw 0321 % with this definition of A if sa or sb are changed from the 0322 % default values 0323 if 3<=N & N<=50 & 2<=M & M <=9.5 & 1<= gammai & gammai<=20 0324 f1NM = 4.1*(N-2*M^0.28+5.3)^(-1.45*M^0.1+0.96); 0325 f2NM = (2.2*M^(-3.3) + 0.57)*N^(-0.58*M^0.37+0.53)-1.04*M^(-1.9)+0.94; 0326 A = (1+ f1NM*log(gammai)^f2NM)/gammai; 0327 elseif N == 5 & M == 4, 0328 A = (1-0.287*log(gammai)); 0329 else 0330 error('Not knowing the normalization because N, M or peakedness parameter was out of bounds!') 0331 end 0332 end 0333 0334 S.S = S.S*A; 0335 0336 0337
Comments or corrections to the WAFO group