Irregular Waves by MATLAB

Ocean wave model is necessary for numerical studies of marine crafts and structures. The process of wave generation due to wind starts with small wavelets. These short waves continue to grow until they finally break and their energy is dissipated. A developing sea starts with high frequencies creating a spectrum with a peak at a relatively high frequency. A developed sea is created when a storm has lasted for a long time. After the wind has stopped, a low-frequency decaying sea or swell is formed. This study will show you wave generation codes using Pierson-Moskowitz Wave Spectrum with MATLAB.

Pierson and Moskowitz have developed a two-parameter wave spectral formulation for fully developed wind-generated seas from analyses of wave spectra in the North Atlantic Ocean:

(1)    \begin{equation*} S(\omega)=A \omega^{-5} \text{exp}(-B \omega^{-4}). \end{equation*}

The parameters are

(2)    \begin{equation*} A=8.1 \times 10^{-3} g^2, \end{equation*} \begin{equation*} B=0.74\big(\frac{g}{V_{19.4}}\big)^4, \end{equation*}

where  $V_{19.4}$ is the wind speed at the height 19.4 m. The modal frequency is

(3)    \begin{equation*} \omega_0=\big(\frac{4B}{5} \big)^{0.25}. \end{equation*}

In developed sea, significant wave height increase as the square of wind speed that  $H_s=2.06V_{10}^2/g$ . The parameters A and B can also be calculated by

(4)    \begin{equation*} A=\frac{4\pi^3*Hs^2}{Tz^4}, \end{equation*} \begin{equation*} B=\frac{16\pi^3}{Tz^4}, \end{equation*}

The wave elevation of a long-crested irregular sea propagating along positive x-axis can be written as the sum of a large number of wave components, i.e.

(5)    \begin{equation*}  \zeta=\sum_{j=1}^N \sqrt{2S(\omega_j)\Delta \omega} \sin(\omega_j t - k_j x + \epsilon_j). \end{equation*}

Here,  $\Delta \omega,\ k_j$ and  $\epsilon_j$ are the constant difference between successive frequencies, wave number and random phase angle of wave component. The random phase angles  $\epsilon_j$ are uniformly distributed between 0 and  $2 \pi$ and constant with time.

Wave spectrum generation by MATLAB can be expressed based on wind speed

function [S] = waveSpectrum(Vwind,w) % Vwind is the wind speed at height 19.4 m 
% w is the list of angular frequencies, for example: w=0.01:0.01:4. 
    g=9.81; 
    A=8.1e-3*g^2; 
    B=0.74*(g/Vwind)^4; 
    S=A*w.^(-5).*exp(-B./(w.^(4))); 
end 

or significant wave height

function [S] = waveSpectrum(Hs,Tz,w)
     A = 4*pi^3*Hs.^2./(Tz.^4);
     B = 16*pi^3./(Tz.^4);
     S=A*w.^(-5).*exp(-B./(w.^(4)));
end

One Dimensional Waves

The elevation of irregular waves propagating along x-axis can be generated by

function [zeta] = waveGen(S,w,x,t)
% S is the wave spectrum
% w is the list of angular frequencies, for example: w=0.01:0.01:4.
% x is the location of wave from the origin
% t is the current time

persistent Ph                  % random angle phase should be 
   if isempty(Ph)                 % defined only one time
       N=length(w);
       Ph=rand(1,N)*2*pi;         % random angle phases   
   end

   dw=w(2)-w(1);
   g=9.81;
   k=w.^2./g;
   A=sqrt(2*S*dw);             % magnitude of wave
   zeta=sum(A.*sin(w.* t-k.*x+Ph));
end

After saving the above function, the following example can generate the waves

Vwind=12;              % wind speed at 19.4 m
dw=0.01;                % the difference between successive frequencies
w=0.01:dw:4;            % angular frequencies
t=0.1:0.02:500;         % simulation time
x=0;

N=length(t);
zeta=zeros(1,N);

S=waveSpectrum(Vwind,w);
figure;
plot(w,S);
xlabel('frequency [rad/s]');
ylabel('wave spectrum [m^2/s]');
grid on;
for i=1:N
    zeta(i)=waveGen(S,w,x,t(i));
end

figure;
plot(t,zeta);
xlabel('time [s]');
ylabel('wave elevation [m]');
grid on;

Two Dimensional Waves

Wind driven waves would be propagated in all direction relative to the wind. For each angle $ \phi$ relative to the wind direction, the spectral density  $S(\omega,\phi) $ is usually expressed

(6)    \begin{equation*} S(\omega,\phi)=S(\omega)f(\phi) \end{equation*}

where  $f(\phi)$ is defined as

(7)    \begin{equation*} f(\phi)=\left{ \begin{array}{ll} A_n \cos^n(\phi), \ & |\phi| \le \frac{\pi}{2} \ 0, & \text{otherwise} \end{array} \right. \end{equation*}

with  n=2 \ \text{or} \ 4 , and  A_2=2/\pi \ \text{and} \ A_4=8/(3\pi) . The generation of 2d irregular waves in MATLAB can be realized as follows:

function [S] = waveSpectrum2d(Vwind,w,phi,n)
% Vwind is the wind speed at height 19.4 m
% w is the list of angular frequencies, for example: w=0.01:0.01:4.
% phi is the wave propagation angle list relative to the x-axis, e.g., phi=-pi:0.01:pi
% n: n=2 or n=4

g=9.81;
A=8.1e-3*g^2;
B=0.74*(g/Vwind)^4;
S0=A*w.^(-5).*exp(-B./(w.^(-4)));

if(n~=2 && n~=4)
   error('n must be 2 or 4');
end
if(n==2)
   An=2/pi;
elseif(n==4)
   An=8/(3*pi);
end
f=An*cos(phi).^n.*(abs(phi) < pi/2);
if(iscolumn(S0)==0)
   S0=S0';
end
if(isrow(f)==0)
   f=f';
end
S=kron(S0,f);
end
function [el] = waveGen2d(S,w,phi,x,y,t)
% S is the wave spectrum
% w is the list of angular frequencies, for example: w=0.01:0.01:4.
% phi is the wave propagation angle relative to the x-axis
% x, y are the location of wave from the origin
% t is the current time

persistent Ph                  % random angle phase should be 
if isempty(Ph)                 % defined only one time
    M=length(w);
    N=length(phi);
    Ph=rand(M,N)*2*pi;         % random angle phases   
end

dw=w(2)-w(1);
dphi=phi(2)-phi(1);
g=9.81;
k=w.^2./g;
A=sqrt(2*S*dw*dphi);             % magnitude of wave
el=0;
for i=1:length(w)
   el=el+sum(A(i,:).*sin(w(i)*t+k(i)*(x*cos(phi)+y*sin(phi))+Ph(i,:)));
end
end

Example:

clear;
Vwind=12;              % wind speed at 19.4 m
dw=0.01;                % the difference between successive frequencies
w=0.01:dw:4;            % angular frequencies
dphi=0.02;
phi=-pi:dphi:pi;
t=0:0.1:50;         % simulation time
x=-20:0.5:20;
y=-20:0.5:20;
Nx=length(x);
Ny=length(y);

Nt=length(t);
el=zeros(Nx,Ny,Nt);

n=2;
S=waveSpectrum2d(Vwind,w,phi,n);
figure;
surf(w,phi,S');
xlabel('frequency [rad/s]');
ylabel('phi [rad]');
zlabel('spectrum [m/s^2]');
grid on;
for i=1:Nt
    for j=1:Nx
        for k=1:Ny
            el(j,k,i)=waveGen2d(S,w,phi,x(j),y(k),t(i));
        end
    end
end

figure;
surf(x,y,el(:,:,1)');
xlabel('x [m]');
ylabel('y [m]');
zlabel('wave elevation [m]');
grid on;

Leave a Comment

Your email address will not be published. Required fields are marked *