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)

The parameters are

(2)

where is the wind speed at the height 19.4 m. The modal frequency is

(3)

In developed sea, significant wave height increase as the square of wind speed that . The parameters A and B can also be calculated by

(4)

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)

Here, and are the constant difference between successive frequencies, wave number and random phase angle of wave component. The random phase angles are uniformly distributed between 0 and 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 relative to the wind direction, the spectral density is usually expressed

(6)

where is defined as

(7)

with , and . 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;
```