Irregular ocean wave generation codes for MATLAB/Simulink

Contents
  • Software for this post
  • Introduction
  • Wave spectrum (PM Spectrum)
  • Description of irregular wave by MATLAB (one direction)
  • An example
  • Waves in two dimensional plane
  • Conclusion
Software for this post

MATLAB/Simulink (verified by version R2015b).

Introduction

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/Simulink.

 Pierson-Moskowitz Spectrum

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

     \[ \omega_0=\big(\frac{4B}{5} \big)^{0.25}. \]

The MATLAB codes can be expressed as

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 
Description of Waves

Linear theory is usually used to simulate irregular wave. 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.

(3)    \begin{equation*}  \zeta=\sum_{j=1}^N A_j \sin(\omega_j t - k_j x + \epsilon_j). \end{equation*}

Here,  A_j,\ \omega_j,\ k_j and  \epsilon_j mean the wave amplitude, circular frequency, wave number and random phase angle of wave component  j . The random phase angles  \epsilon_j are uniformly distributed between 0 and  2 \pi and constant with time. For deep water,  \omega_j and  k_j are related by the dispersion relationship

(4)    \begin{equation*} \frac{\omega_j^2}{g}=k_j. \end{equation*}

In addition, the wave amplitude can be expressed by

(5)    \begin{equation*} A_j=\sqrt{2 S(\omega_j)\Delta \omega}, \end{equation*}

where  \Delta \omega is a constant difference between successive frequencies.

The MATLAB codes for computing the elevation of irregular waves can be expressed

function [el] = 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
el=sum(A.*sin(w.* t-k.*x+Ph));
end
An Example

Save the above functions in the same directory, and then run the following commands in MATLAB window to test the wave elevation at  x=0 by setting the wind speed at the height of 19.4 m as 12 m/s.

clear;
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);
el=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
    el(i)=waveGen(S,w,x,t(i));
end

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

The wave spectrum and wave elevation are shown in the following figures.

wave spectrum

The above files including a Simulink model can be downloaded from here.

Waves in two Dimensional plane

In the above, we only discussed wave propagating along the x-axis. 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 MATLAB codes can be expressed as

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

Then, the wave generation codes in 2d plane can be written as

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

An example of using the above functions to generate waves is given as follows.

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;

The generated wave at  t=0~s is shown in the following figure.

Conclusion

Irregular wave generation based on PM spectrum was studied and codes for MATLAB were presented. The functions can be simply utilized for the simulation in Simulink using MATLAB function block.

The functions for generating waves in 2d plane can be download here.

Leave a Reply

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