Hgen.m

For the purpose of testing our program frequently with data during the writing process, we wrote an m-file that generates random irregular water level data. 

% Generates a data column of Water level...

N = 20000;
dt=0.1;
t=(0:N-1)*dt;
irr = 0.3*sin(2*pi/2*t)*rand(N,1)+rand(N,1);
w1 = 0.1*sin(2*pi/2*t);
w2 = 0.3*sin(2*pi*t/N).*sin(2*pi/4*t);
trend=2+w1+w2;

H = trend+irr;

This file aims to get a wide spread of wave heights.  Just a sine function with a random noise creates an almost constant wave height.  To get a wave simulation that varies with time, we multiplied two sine functions, one with a high frequency and one with a period that spans the whole data set.

The component that includes a random function multiplied by a sine function was put in to simulate choppiness that varies over time, such as occasional storms or wind gusts.
--------------------------------------------------------------------------------------------------------------------------------------------

Panalyzer.m

Our main program takes water level data, with other known parameters entered by user, calculates the wave heights, determines the 100 year wave height through statistical analysis, calculates the corresponding wave period and wave length based on the lake size, determines if the wave would be breaking, and then it displays the various pressures that would act on a hypothetical wave wall as functions of height from ocean/lake floor.

%========Calculates the Wave heights of a given data set==============

Z = input('enter the water height data (1 column)  ');
dt = input('enter the time interval of the data  ');
d = input('enter the depth at the toe of the structure  ');
n = size(Z);
n = n(1,1);
F = input('enter the fetch length up to the structure  ');

% Eliminate Trend (water level fluctuation)
for i = 1:n;
    x(i,1) = i;
    Trend = polyfit(x,Z,1);
    eta(i,1) = Z(i,1) - Trend(1)*i-Trend(2);
end

%Down crossing analysis
z=0;
for i = 2:n
   if eta(i-1) >= 0
      if eta(i) <= 0
         z=z+1;
         zc(z)=i;
      end
   end
end
nw=z-1;
Tavg=n*dt/nw;
for j = 1:nw
      h(j) = max(eta(zc(j):zc(j+1)))-min(eta(zc(j):zc(j+1)));
      t(j) = j*dt;
end
t = t(1:nw);

figure
plot(t,h)
xlabel('Time (sec)')
ylabel('Wave Height (m)')
title('Wave Heights')

This first portion of the script takes the water level data and determines the wave heights.      It first translates the Z data to "eta" by subtracting the mean water level trend over the time period.  In the down crossing loop, 'z' becomes the number of times that eta crosses zero.  'zc' is an array that records each 'i' where eta crosses zero, so that it can be referenced in the eta array.  The next loop defines the wave height as simply the max eta in each range minus the minimum of the range between 'zc' increments.  Finally, it plots the wave height as a function of time.
This part of the program was modeled after "wavan.m" by Chin Wu.

%==============Wave height statistics============================
minh = min(h);
hmax = max(h);
havg = mean(h);

%Puts wave Heights into b bins
Occ=zeros(50,1);
Height=zeros(50,1);
b = 50;
for i = 1:b;
    Height(i,1) = hmax/b*i;
    for u = 1:nw;
        if h(u) >= Height(i,1);
            Occ(i)=Occ(i)+1;
        end
    end
end

%Computes the probability of the wave height being <= bin height
Cum_Occ=zeros(b,1);
P = zeros (b,1);
Cum_Occ(1)=Occ(1);
for i = 2:b;
    Cum_Occ(i)=Occ(i)+Cum_Occ(i-1);
end
for i = 1:b;
    P(i) = Cum_Occ(i)/(Cum_Occ(b)+1);
end
Z=norminv(P,0,1);
Height=Height(1:b,1);
LH=log(Height);

figure
subplot(2,1,1)
plot(Height,Z)
xlabel('Wave Height')
ylabel('Reduced Variate Z')
title('Normal Distribution')

subplot(2,1,2)
plot(LH,Z)
xlabel('Log of Height')
ylabel('Reduced Variate Z')
title('Log Normal Distribution')

Po=input('Normal or Lognormal fit? enter 1 or 2 respectively  ');

This next segment of the script divides the wave heights into 50 equal sized bins.  The 'Occ' array stands for "occurrence," and it says how many wave heights are greater than or equal to the 'Height' bin size.  'Cum_Occ' stands for "cumulative occurrence;" it gives the cumulative amount of wave heights at or below the defined 'Height' bin.  'P' is the corresponding probability of the wave height falling into each 'Height' bin.  'Z' is the reduced variate corresponding to each probability.  The last part of this section graphs the normal distribution and lognormal distribution of the wave heights.  Since the most straight plot is ideal, the program allows the user to determine which distribution to use at this point.

%=================Calculates 100 year wave height==============

Time=n*dt;
Wpy=31536000*nw/Time;  % Average number of events per year
PH=1-1/(100*Wpy);      % Probability of 100 year event
ZH=norminv(PH,0,1);

if Po == 1
    Coef=polyfit(Height,Z,1);
    HYH=(ZH-Coef(2))/Coef(1);
else
    Coef=polyfit(LH,Z,1);
    HYH=exp((ZH-Coef(2))/Coef(1));
end
HYH

This part of the file determines the 100 year wave height based on the regression of the of the chosen distribution.  'HYH' is the 100 year wave height.

%===================Calculates the Period======================

Uf = (3832031.25*HYH^2/F)^0.5;
Ut = (163255*HYH^4)^(1/5);
if Uf >= Ut  % Fetch limited (higher required wind for fetch analysis)
    U = Uf;
    Fs = 9.81*F/U^2;
else
    U = Ut;  % Limited by a duration of 4 hours
    Fs = (2053.256/U)^(3/2);
end
Tp = U*0.02915*(Fs)^(1/3)

Since the program aims to measure the pressures resulting from a 100 year occurrence wave, it wouldn't work to simply take the period calculated directly from the data.  Assuming that the user would know the approximate size of the body of water, the script does a reverse of the JONSWAP method to determine the required wind speed to generate a wave height equal to the hundred year height already calculated.

This script determines the required wind speed for a fetch-limited scenario and a 4 hour duration-limited scenario.  The scenario that requires the higher wind speed is the limiting factor, so the program choses that scenario and determines the corresponding wave period for the wind speed.

%===============Calculates the Wavelength======================
Lo = 1.56*Tp^2;
L1 = Lo*(tanh(2*pi*d/Lo))^0.5;
L2 = Lo*tanh(2*pi*d/L1);
L = (L1 + L2)/2

This section goes through the standard iteration method of calculating the wave length when period and depth are given.

%==============Determines if wave is breaking==================

Bc = tanh(2*pi*d/L)/7;
slope = HYH/L;
if slope >= Bc
    disp('Wave is breaking; data validity is comprimised')
end

Since the plots of pressure on the wall don't take into account wave breaking conditions, this section determines of the wave would be breaking based on the "Miche" breaking criteria, and displays to the user whether or not the 100 year wave would be breaking at the location (depth) where the hypothetical wave wall would be.

%====Calculates and plots the Pressures acting on seawall=======

ro = 1000; % density of water (kg/m^3)
s = pi*HYH^2/L*coth(2*pi*d/L);   % vertical run-up on the wall
pb = ro*9.81*HYH/cosh(2*pi*d/L)/1000;   % Wave pressure at bottom (KPa)
pc = (pb+ro*9.81*d)*(HYH+s)/(d+HYH+s)/1000;  % Wave pressure at crest (KPa)
pt = ro*9.81*(HYH-s)/1000;   % Wave pressure at trough (KPa)
ps = ro*9.81*d/1000;   % Hydrostatic pressure at bottom (KPa)

G=300;
y=zeros(G,1);
for i = 1:G
    y(i,1)=i/G*(d+HYH+s);
end
E = G*d/(d+HYH+s);
e = floor(E);    % determines point of mean water level
xs = zeros(G,1);
for i = 1:e
    xs(i,1)=ro/1000*9.81*d-i/e*(ro/1000*9.81*d);
end
xw = zeros(G,1);
for i = 1:e
    xw(i,1) = pb+(pc-pb)/e*i;
end
for i = e:G
    xw(i,1) = pc-((i-e)/(G-e))*pc;
end
ge = G*(d+s-HYH)/(d+s+HYH);
g = floor(ge);     % determines the point of water level at trough
xt = zeros(G,1);
for i = 1:g
    xt(i,1) = -pb-(pt-pb)/g*i;
end
for i = g:e
    xt(i,1) = -pt + (i-g)/(e-g)*pt;
end
x1 = xs + xw; % water pressure distribution at crest
x2 = xs + xt; % water pressure distribution at trough

figure
subplot(1,3,1)
plot(xs,y)
xlabel('Pressure (KPa)')
ylabel('Height from bottom')
title('Hydrostatic Pressure')

subplot(1,3,2)
plot(xw,y)
xlabel('Pressure (KPa)')
ylabel('Height from bottom')
title('Pressure from waves at crest')

subplot(1,3,3)
plot(xt,y)
xlabel('Pressure (KPa)')
ylabel('Height from bottom')
title('Pressure from waves at trough')

figure
subplot(1,2,1)
plot(x1,y)
xlabel('Pressure (KPa)')
ylabel('Height from bottom')
title('Pressure Distribution for wave crest')

subplot(1,2,2)
plot(x2,y)
xlabel('Pressure (KPa)')
ylabel('Height from bottom')
title('Pressure Distribution for wave trough')

This sections calculates and plots the pressures as a function of height based on the methods and equations in the Coastal Engineering Manual by the Army Corps of Engineers.  These equations and distributions are displayed in the diagram below.  Each of the parameters are labeled in the script, and the generated graphs are shown in the "Results" section of this website.