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:
 createspec Spectrum structure constructor gravity returns the constant acceleration of gravity simpson Numerical integration with the Simpson method wspecplot Plot a spectral density error Display message and abort function. gamma Gamma function. hold Hold current graph. ishold Return hold state. linspace Linearly spaced vector. num2str Convert number to string. (Fast version) strcat Concatenate strings.
This function is called by:
 Chapter1 % CHAPTER1 demonstrates some applications of WAFO Chapter2 % CHAPTER2 Modelling random loads and stochastic waves Chapter3 % CHAPTER3 Demonstrates distributions of wave characteristics itmkurs_lab3 Script to computer exercises 3 thspdf2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. wafoinit setup all global variables of the WAFODEMO

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