Tutorial 2.5: Wave simulation

In this tutorial you will learn how to simulate a wave elevation using the wave spectrum and assuming potential theory to describe the water kinematics.

In this example we assume one wave direction (long-crested waves) only.

NOTE: you will require wave_spectrum.py
IMPORTANT: Download the required external functions and put them in the same folder as this notebook
rdeke/ComModHOS_double

Part 1: Calculate harmonic wave component parameters from a wave spectrum.

We start by defining the numerical values of all the parameters that will be needed:

• spectrum_type

  • SpecType =1 - Bretschneitder (p1=A,p2=B)

  • SpecType =2 - Pierson-Moskowitz (p1=Vwind20)

  • SpecType =3 - ITTC-Modified Pierson-Moskowitz (p1=Hs,p2=T0)

  • SpecType =4 - ITTC-Modified Pierson-Moskowitz (p1=Hs,p2=T1)

  • SpecType =5 - ITTC-Modified Pierson-Moskowitz (p1=Hs,p2=Tz)

  • SpecType =6 - JONSWAP (p1=Vwind10,p2=Fetch)

  • SpecType =7 - JONSWAP (p1=Hs,p2=w0,p3=gamma)

  • SpecType =8 - Torsethaugen (p1=Hs,p2=w0)
    

• hs - Significant wave heigh in sea state [m]

• T0 - Spectrum period at peak frequency in spectrum [s]

• omega_peak - Peak frequency in spectrum [rad/s]

• depth - Average water depth, for calculation of wave numbers [m]

• nfreq - Number of frequency components [-]

• omega_max - Cutoff frequency of spectrum = cutoff*omega_mean [rad]

• rand_seed - Random number seed, applies to all random numbers (phase, rand frequency)

from wave_spectrum_generator import wave_spectrum
import numpy as np
import matplotlib.pyplot as plt

# Parameters
spectrum_type = 3
hs = 1.0
T0 = 9
omega_peak = 2*np.pi/T0
depth = 20
nfreq = 20
omega_max = 3*omega_peak
rand_seed = 584
g = 9.81

With the given parameters we need to generate the following harmonic wave component parameters required for the time-simulations.

• zeta_a - Vector of harmonic wave amplitudes

• omega - Vector of harmonic wave frequencies

• phase - Vector of harmonic wave phases (random)

• wavenum - Vector of harmonic wave numbers

# Frequency step
delta_omega = omega_max/nfreq

# Maximum simulation time (before the signal gets repeated)
max_time_sim = 2*np.pi/delta_omega

print("Max simulation time: "+str(max_time_sim))

# Frequency vector, starting at delta_omega
omega_vec = np.arange(delta_omega,omega_max + delta_omega,delta_omega)

Max simulation time: 60.0

Set random generator required to create the same set of random number based on the seed number. If this line is commented or unused, everytime we run the code we will have different random numbers.

np.random.seed(rand_seed)

# Create evenly distributed random phases
phase = np.random.randn(1,nfreq)*2*np.pi

# Generate the spectral densities using the provided function
s_omega = wave_spectrum(spectrum_type, [hs, T0], omega_vec, 0)

Use: S = wave_spectrum(SpecType, Par, W, PlotFlag)

Input:

SpecType - Spectrum type

Par - Spectrum parameters

W - List of wave frequencies [rad/s]

PlotFlag - 1 to plot the spectrum, 0 for no plot

Output:

S - List of wave spectrum values [m^2 s]; evaluated at W[k]

We can plot the spectrum and the frequencies employed:

plt.figure()
plt.plot(omega_vec, s_omega)
plt.plot(omega_vec, s_omega, 'r*')
plt.xlabel("$\omega$ [rad/s]")
plt.ylabel("S [m$^2$s]")
plt.title("Long-crested wave spectrum, "+str(nfreq)+" harmonic components as red stars");

You can try using different number of harmonic components and check the simulation time.

Now we need to calculate the wave amplitude and wave number for every harmonic component.

from scipy.optimize import root_scalar
zeta_a = []
wavenum = []
for ii in np.arange(0, nfreq):
    # Calculate the wave amplitude.
    # The area under the spectrum represents 0.5*zeta_a^2
    # here we assume that this area is simply calculated by the S(w)*Deltaw.
    zeta_a.append(np.sqrt(2*s_omega[ii]*delta_omega))

    # Calculate the wave numbers using the dispersion relation
    # we are going to use the root_scaler function to solve the non-
    # linear equation
    wavenum_0 = omega_vec[ii]**2/g
    def func(k):
        return k*np.tanh(k*depth) - omega_vec[ii]**2/g
    sol = root_scalar(func, method='toms748', bracket=[wavenum_0, 1e10])
    wavenum.append(sol.root)
omega = omega_vec

Now we have defined all our required parameters!

Part 2: Create time domain simulation of wave eleveation

List of time values for simulation Note: depending on number of frequency components (nfreq), the time response will be repeated after \(2\pi/\)delta_omega

t = np.arange(0, 100.1, 0.1)

# Define position
x = 0.0

For each time step the wave elevation is obtained as the superpostiion of the different harmonic components. Please note that it can be done in a vectorial way for all frequencies and is not necessary to use a for loop.

zeta = []
for i in np.arange(0, len(t)):
    zeta.append(np.sum(zeta_a*np.cos(omega*t[i] + phase - np.array(wavenum)*x)))

plt.figure()
plt.plot(t, zeta)
plt.xlabel("t [s]")
plt.ylabel("Surface elevation $zeta$ [m]");

---

Problem: Calculate the heave force and roll moment on a rectangular barge of dimensions 50x50 m. Plot the results in the time domain over an adequate time period.

Hint: You will need to make use of potential theory to calculate the pressure of the water underthe barge. At this point, this step should be straightforward since the wave parameters for the different harmonic components are already calculated in this tutorial. Note that the position where the pressure needs to be determined depends on to the motion of the structure.

---

Some additional parameters:

• Barge dimensions: 50x50x4 m, in deep water

• Barge draft: 0.5 m in seawater \(\rho\) = 1025 \(kg/m^3\), assume uniform weight distribution

• Hydrodynamic added mass (in heave and roll) = 10% and damping = 1% of dry mass

• \(H_s\) = 1.0 m, \(T_0\) = 9 s, g = 9.81 N/kg (Same as before)

# Implement here ...

---

The solution can be found here.