# 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

##### 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) The parameters are

(2) where is the wind speed at the height 19.4 m. The modal frequency is 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) Here, and mean the wave amplitude, circular frequency, wave number and random phase angle of wave component . The random phase angles are uniformly distributed between 0 and and constant with time. For deep water, and are related by the dispersion relationship

(4) In addition, the wave amplitude can be expressed by

(5) where 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 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);
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. ##### 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 relative to the wind direction, the spectral density is usually expressed

(6) where is defined as

(7) with , and .

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');
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 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. 