Home > wafo > spec > torsethaugen.m

torsethaugen

PURPOSE ^

Calculates a double peaked (swell + wind) spectrum

SYNOPSIS ^

[S, Sw, Ss]=torsethaugen(w1,sdata,plotflag)

DESCRIPTION ^

  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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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