This notebook computes the dynamics and limits of a paddle-style wavemaker. First we define some constants:
In [30]:
# Some constants
h = 2.44 # Water depth in meters
w = 3.66 # Tank width in meters
rho = 1000.0 # Water density
g = 9.81 # Acceleration due to gravity
m_paddle = 1270.0 # Paddle mass in kg, 1270.0 from OMB manual
h_piston = 3.3147 # Piston attachment height in meters
max_halfstroke = 0.16 # Maximum stroke amplitude of piston
min_period = 0.5 # Minimum period limit
max_period = 5.0 # Maximum period limit
Next, we write functions to solve the dispersion relation (Dean and Dalrymple)
$$ \sigma^2 = gk \tanh kh. $$
In [31]:
import numpy as np
from numpy import pi
from __future__ import division, print_function
def wavenumber(rad_freq, depth=2.44, decimals=5, verbose=False):
k = np.arange(0, 30, 10**(-1*decimals))
residuals = np.abs(rad_freq**2 - g*k*np.tanh(k*depth))
k = k[residuals==residuals.min()][0]
if verbose:
print("k =", k, "\nresidual =", residuals.min())
return k
# Test this function
period = 1.0
rad_freq = 2*np.pi/period
depth = 2.44
k = wavenumber(rad_freq, depth, verbose=True)
print("L =", 2*pi/k, "m")
Next, we write the inverse function, to solve for radian frequency
In [32]:
def frequency(wavenumber, depth=2.44, decimals=5, verbose=False):
k = wavenumber
sigma = np.arange(0, 10, 10**(-1*decimals))
residuals = np.abs(sigma**2 - g*k*np.tanh(k*depth))
sigma = sigma[residuals==residuals.min()][0]
if verbose:
print("sigma =", sigma, "\nresidual =", residuals.min())
return sigma
# Now see if this function is consistent with the wavenumber function
sigma_calculated = frequency(k, depth, verbose=True)
print("sigma_nominal =", rad_freq)
print("error =", (sigma_calculated - rad_freq)/rad_freq)
Dean and Dalrymple provide a formula (Eq. 6.25) for the ratio of wave height $H$ to piston stroke $S$ (at the water surface)
$$ \frac{H}{S} = 4 \left( \frac{\sinh k_p h}{k_p h} \right) \frac{k_p h \sinh k_p h -\cosh k_p h + 1} {\sinh 2 k_p h + 2 k_p h}, $$where $k$ is the wavenumber (the $p$ subscript meaning progressive waves), and $h$ is the water depth. With this formula we can write functions for converting stroke amplitude to wave height and vice versa:
In [33]:
def stroke_amp_to_height(piston_stroke_amp, period, paddle_height=3.3147, depth=2.44):
sigma = 2*pi/period
k = wavenumber(sigma, depth)
kh = k*depth
surface_stroke = 2*piston_stroke_amp*depth/paddle_height # Note the multiplication by 2
height = surface_stroke*(4*(np.sinh(kh)/(kh))*(kh*np.sinh(kh) \
- np.cosh(kh) + 1)/(np.sinh(2*kh) + 2*kh))
return height
def height_to_stroke_amp(height, period, paddle_height=3.3147, depth=2.44):
"""Note that this returns the piston stroke amplitude, not the stroke at the water surface."""
sigma = 2*pi/period
k = wavenumber(sigma, depth)
kh = k*depth
surface_stroke = height/(4*(np.sinh(kh)/(kh))*(kh*np.sinh(kh) \
- np.cosh(kh) + 1)/(np.sinh(2*kh) + 2*kh))
piston_stroke = surface_stroke*paddle_height/depth/2
return piston_stroke
Wave height is limited by a number of factors (for the UNH system, these were copied from the "Regular Waves" LabVIEW program):
We can compute the maximum wave height limit by choosing the lowest limit computed using these criteria, given a desired wave period:
In [34]:
def calc_height_limit(period, max_stroke_amp=0.16, max_H_L=0.1, max_H_h=0.65, depth=2.44):
# Compute maximum wave height with maximum piston stroke amplitude
wh1 = stroke_amp_to_height(max_stroke_amp, period)
# Now using max H/L
L = 2*pi/wavenumber(2*pi/period, decimals=3)
wh2 = max_H_L*L
# Now max H/h
wh3 = max_H_h*depth
return np.min([wh1, wh2, wh3])
In [35]:
import matplotlib.pyplot as plt
%matplotlib inline
# Now lets plot these over our range of periods
periods = np.linspace(0.5, 5, num=20)
max_heights = []
for period in periods:
max_heights.append(calc_height_limit(period))
plt.figure()
plt.plot(periods, max_heights)
plt.xlabel("Period (s)")
plt.ylabel("Max wave height (m)")
plt.show()
Instantaneous paddle torque exerted by the piston is a result of the paddle's angular acceleration, and the torque imparted by the water, which we will assume is the pressure (neglecting hydrostatic) integrated over the wetted area multiplied by the center of pressure $(h - z_c)$, i.e.,
$$ F_\mathrm{piston} h_\mathrm{piston} = I_\mathrm{paddle} \ddot{\theta} + (h - z_c) \int p \, \mathrm{d} A. $$For the pressure distribution we will use the formula given in Dean and Dalrymple for linear propagating waves (Eq. 4.22)
$$ p_p = -\rho g z + \rho g \frac{H}{2}\frac{\cosh k(h+z)}{\cosh kh} \cos (kx - \sigma t). $$Assuming the paddle amplitude is small and hydrostatic pressure imparts no net force, our effective pressure
$$ p = \rho g \frac{H}{2}\frac{\cosh k(h+z)}{\cosh kh} \cos (-\sigma t). $$The force imparted by the water is then the average pressure multiplied by the paddle area, which for engineering purposes we will assume to be the tank cross section, giving us
$$ F_\mathrm{water} = p_\mathrm{ave} wh $$We can assume the paddle moves in a sinusoidal fashion, therefore
$$ F_\mathrm{piston} = \frac{1}{h_\mathrm{piston}} \left( I_\mathrm{paddle} \sigma^2 \frac{S_\mathrm{piston}}{2} \sin \sigma t - (h - z_c)F_\mathrm{water} \right). $$Now, we will write this into a function to compute the paddle dynamics for a give wave height and period:
In [36]:
def calc_dynamics(height, period, verbose=False):
sigma = 2*pi/period
k = wavenumber(sigma, depth=2.44, decimals=3)
c = sigma/k # Wave celerity
c_g = c*1/2*(1 + 2*k*h/np.sinh(2*k*h)) # Group velocity
piston_stroke_amp = height_to_stroke_amp(height, period)
S_p = 2*piston_stroke_amp
wave_energy = 1/8*rho*g*height**2*w
wave_power = wave_energy*c_g
# Create arrays to generate a time series for one period
t = np.linspace(0, period)
z = np.linspace(0, -h)
z_c = np.zeros(len(t))
p_ave = np.zeros(len(t))
p = np.zeros((len(z), len(t)))
for n in range(len(z)):
"""Calculate pressure at each height"""
p[n,:] = rho*g*height/2*np.cosh(k*(h + z[n]))/np.cosh(k*h)*np.cos(-sigma*t)
for n in range(len(t)):
A = np.trapz(p[:,n], z)
p_ave[n] = A/h
z_c[n] = 1/A * np.trapz(z*p[:,n], z)
force_water = p_ave*w*h
I_paddle = 1/3*m_paddle*h_piston**2
force_inertial = I_paddle*sigma**2*S_p/2*np.sin(sigma*t)
force_piston = (force_inertial - (h - z_c[0])*force_water)/h_piston
x_piston = S_p/2*np.sin(sigma*t)
vel_piston = sigma*S_p/2*np.cos(sigma*t)
acc_piston = -sigma**2*S_p/2*np.sin(sigma*t)
power_piston = force_piston*vel_piston
data = {"mean piston power" : power_piston.mean(),
"mean wave power" : wave_power,
"max piston power" : power_piston.max(),
"max piston force" : force_piston.max(),
"max piston displacement" : x_piston.max(),
"max piston acceleration" : acc_piston.max(),
"max piston velocity" : vel_piston.max(),
"max water force" : force_water.max(),
"max interial force" : force_inertial.max()}
if verbose:
for k,v in data.iteritems():
print(k, "=", v)
return data
# Test the function
d = calc_dynamics(0.1, 1.0, verbose=True)
In [37]:
# Now let's compute the dynamics over the entire envelope of operation at max wave height
periods = np.linspace(0.5, 5.0, num=30)
data = {"max piston displacement" : [],
"max piston force" : [],
"max piston velocity" : [],
"max piston acceleration" : [],
"mean piston power" : [],
"mean wave power" : []}
for period in periods:
height = calc_height_limit(period)
d = calc_dynamics(height, period)
for key, val in d.iteritems():
if key in data:
data[key].append(val)
units = {"force" : "N",
"velocity" : "m/s",
"acceleration" : r"m/s$^2$",
"power" : "W",
"displacement" : "m"}
plt.figure(figsize=(8,10))
n = 1
for key in data.keys():
plt.subplot("32{}".format(n))
plt.plot(periods, data[key])
plt.title(key.capitalize() + " (" + units[key.split()[-1]] + ")")
plt.xlabel("Period (s)")
n += 1
plt.tight_layout()
plt.show()
A simulation was performed using OpenFOAM's interDyMFoam solver in 2-D, where a max-height 2.2 second period wave was generated. Moment on the flapper was calculated, and then inertial forces were added using the above assumption of moment of inertia. Below is a video of the simulation. Note that there is no beach, so by the end we have set up a standing wave. OpenFOAM case files are available at https://github.com/petebachant/waveFlapper-OpenFOAM.
In [38]:
from IPython.display import YouTubeVideo
YouTubeVideo("lfOAxYMy60c", width=600, height=500, vq="large")
Out[38]:
Next we can compare the flapper moment with that of the theoretical calculations. The figure below shows the moment including the flap inertia. Note that for reference, the max moment from the theoretical calculation for a 2.2 second period wave is approximately
In [39]:
5500*h_piston
Out[39]:
In [40]:
from IPython.display import SVG
SVG("images/cfd.svg")
Out[40]:
The CFD results show greater varation of the force over each cycle, but agree pretty well with the theoretical calculation of max moment, which gives us some confidence that we can use these figures to size an actuator.
Dean, R.G. and Dalrymple, R.A., Water Wave Mechanics for Engineers and Scientists. World Scientific