The goal of this notebook is to compute the potential flow and forces on a body in waves, but let's start with a recap.
Recall from the first notebook that the scaled potential $\phi_k$ due to motion in degree of freedom $k$ is defined as a superposition of the contributions from each panel $i$
$$\phi_k = \sum_i q_{k,i} F_i^\phi \quad\text{where}\quad F_i^\phi=\int_{s_i}^{s_{i+1}} G(s) ds $$where $q_{k,i}$ is the panel strength and $G$ is the panel Green's function. Once the strength's are set to satisfy the body boundary condition, $\phi_k$ determines the whole flow and the added mass matrix.
The first notebook developed a numerical method to compute this potential around any submerged 2D geometry using a potential flow source as the Green's function, as well as code to plot the flow and compute the added mass. Those functions are collected in this source file which we can execute using the requests package.
import requests
url = 'https://raw.githubusercontent.com/weymouth/MarineHydro/master/src/PanelMethod.py'
exec(requests.get(url).content)
In this notebook we will consider the oscillation of a body near the free surface, with displacement
$$\xi_k = |\xi_k|\cos\omega t$$where $|\xi_k|,\omega$ are the amplitude and frequency of motion.
This motions near the free surface will generate waves, changing the flow and forces significantly compared to the case considered in the first notebook. Despite the huge change to the physical set-up and results, we only need to make a single modification to the numerical method:
We need to modify the source Green's function $G$ to include the potential flow waves which it generates. That's it!
Here is new the Green's function, taken from this modern paper:
$$ G(x,y,x_s,y_s) = \log(r_s)-2\pi \imath e^Z-f(Z) $$There are three terms:
The first term is the oscillating source where $r_s^2 = (x-x_s)^2+(y-y_s)^2$ is the square-distance from the point $(x_s,y_s)$ on the segment to a point $(x,y)$ in space, as before. As in the first notebook, the Green's function only includes the spacial dependence, the source strength and the oscillation time dependance are multiplied later.
The second term is the traveling wave written using complex variables. The standard wave potential is
where $a,\Phi$ are the wave amplitude and the phase relative to the motion, and $K=\omega^2/g$ is the deep water wave number. In complex variables we have $e^{\imath\theta} = \cos\theta+\imath\sin\theta$, where $\imath=\sqrt{-1}$. Therefore we can write $$ \phi = \Re\left[\underbrace{\zeta e^{K(y+\imath x)}}_\text{space} \underbrace{e^{-\imath\omega t}}_\text{time}\right] $$ where $\zeta$ is the complex amplitude which includes the phase. The traveling wave term in $G$ matches the spacial part of this potential using $Z=K(y+y_s+\imath(x-x_s))$ (shifting the origin of the wave) and setting $\zeta=-2\pi\imath$.
Let's look at the velocity field induced by a single source panel near the free surface.
from scipy.special import exp1
def wave_source(x, y, xs, ys, K):
"Source plus generated free surface waves"
r2 = (x-xs)**2 + (y-ys)**2 # source square-distance
Z = K*(y+ys + 1j*abs(x-xs)) # wave number scaled complex vector
eZ = np.exp(Z) # traveling wave potential
fZ = 2*np.real(eZ*exp1(Z)) + np.log(np.abs(Z)/K) # standing wave potential
return 0.5*np.log(r2) - 2j*np.pi*eZ - fZ
XY = np.meshgrid(np.linspace(-3,3,66), np.linspace(-2,0,22))
K = 4
plot_flow([-.1,.1], [-.77,-0.65], [1], XY,wave_source, (K,), size=(16,6), ymax=0.6)
We can clearly see the source and the waves above, but why are there two sets of vector fields and what do they mean?
The Green's function is now complex, which makes the potential $\phi_k$ and the velocity complex as well. These complex numbers are storing the phase information, just like the complex wave amplitude discussed above.
This is easier to understand using a movie!
Note: It takes a while for Python to generate the videos, so I've commented those lines out. Uncomment and play around with the flows in the following cells as much as you like.
from matplotlib.animation import FuncAnimation, PillowWriter
from IPython.display import Image
def wave_video(x, y, q, XY, G=wave_source, args=(4,), size=(16,6)):
# Get complex velocity
def uv(i): return q[i]*velocity(*XY, x[i], y[i], x[i+1], y[i+1], G, args)
UV = sum(uv(i) for i in range(len(x)-1))
# Plot flow and segments
fig, ax = plt.subplots(1, 1, figsize=size)
Q = ax.quiver(*XY, *UV)#, pivot='mid')
ax.plot(x, y, c='b')
ax.set_ylim(None, 0.5)
ax.set_aspect('equal', adjustable='box')
plt.close()
# run through a wave period and take the real part
def update_quiver(num, Q):
Q.set_UVC(*np.real(UV*np.exp(-2j*np.pi*num/101)))
return Q,
# create the animation
return FuncAnimation(fig, update_quiver, fargs=(Q,))
# anim = wave_video([-.1,.1],[-.77,-0.65],[1],XY)
# anim.save('temp.gif', writer=PillowWriter(fps=20))
# Image(open('temp.gif', 'rb').read())
Looks good. The source, standing wave above it, and the radiating waves on either side are all visible.
Note: Now that you've seen the video, go back and look at the plot above. The black and green arrows are snapshots of the video at the start and a quarter of the way through a cycle. While the movie is certainly easier to digest, it's important to realize the scaled potential $\phi_k$ has all the information needed to make the movie.
Note: The plots above use the original
velocity
function, they just apply the newwave_source
singularity instead of thesource
singularity. Similarly, thepotential
,construct_A
andadded_mass
functions can all by applied to free surface flows without modification!
This is why writing re-usable functions is so powerful. Let's use our new Green's function to produce the flow around a semi-circular ship section in sway.
x, y = ellipse(64, theta1=0)
circleXY = mask_grid(*XY, mask=lambda x,y: x**2 + y**2 < 1)
K = 2
A, nx, ny = construct_A(x, y, wave_source, (K,))
plot_flow(x, y, np.linalg.solve(A,nx), circleXY, wave_source, (K,), size=(14,6), ymax=0.5)
# anim = wave_video(x,y,np.linalg.solve(A,nx),circleXY,wave_source,(K,))
# anim.save('temp.gif', writer=PillowWriter(fps=20))
# Image(open('temp.gif', 'rb').read())
That's amazing! You can clearly see the waves radiating away from the ship hull in either direction. When the ship moves to port, it forms a new crest on the port side, and vice-versa!
Now let's check the added mass:
ma = added_mass(x, y, wave_source, (K,))
print(np.array_str(ma, precision=3, suppress_small=True))
[[ 0.294+0.608j -0. -0.j ] [-0. +0.j 1.135+0.197j]]
Look carefully at this result because two important things have changed since we moved into free surface flows
You can test this by changing the value of $K$ above and you will see a huge variation. Below, I have plotted the dimensionless added-mass and damping coefficients for a range of dimensionless frequencies $\omega \sqrt{R/g}=\sqrt{KR}$.
Note: We've switched to a section in the $y,z$ plane since this is standard for strip theory.
def force_plots(y, z, omega, rhoV=0.5*np.pi):
@np.vectorize
def force_coeffs(omega):
Ca = added_mass(y, z, wave_source, (omega**2,))/rhoV
return Ca[0,0].real,omega*Ca[0,0].imag,Ca[1,1].real,omega*Ca[1,1].imag
A22,B22,A33,B33 = force_coeffs(omega)
plt.plot(omega,A22,'C0',label=r'$\frac{A_{22}}{\rho V}$')
plt.plot(omega,B22,'C0.',label=r'$\frac{B_{22}}{\rho V}\sqrt{\frac{R}{g}}$')
plt.plot(omega,A33,'C1',label=r'$\frac{A_{33}}{\rho V}$')
plt.plot(omega,B33,'C1.',label=r'$\frac{B_{33}}{\rho V}\sqrt{\frac{R}{g}}$')
plt.xlabel(r'$\omega\sqrt{R/g}$')
plt.ylim(0, 2);plt.legend()
force_plots(*ellipse(64,theta1=0), np.linspace(0.1,2))
The added-mass predictions vary from 20% to 200% of the displacement! But is this correct?
Note: You weren't going to trust this numerical method before it was validated... were you?
Here are some experimental and analytic results for the added mass and damping on a semi-circle from Vugts 1968.
Heave | Sway |
---|---|
Our results match the analytic results in the paper really well (other than at two special points that we will discuss in a moment). This validates our free surface numerical predictions!
The potential flow methods both match the experiments nicely other than the zero frequency heave added mass. The linear equations themselves are invalid in the limit of $\omega\rightarrow 0$, causing the discrepancy.
The two blips in our added mass prediction results are due to resonance inside the geometry!
Note: These are the frequencies which could potentially cause issues with sloshing cargo.
The flow inside the hull is plotted below - showing an internal standing wave. This wave totally dominates the flow and causes a singularity in the added mass calculation. However, we're typically only interested in the forces caused by the flow external to the hull, meaning we can ignore these artifacts.
K = 3.25 # sway sloshing wavenumber from plot above
y, z = ellipse(64, theta1=0)
A, ny, nz = construct_A(y, z, wave_source, (K,))
YZ = np.meshgrid(np.linspace(-1,1,22), np.linspace(-1,0,11))
inside = mask_grid(*YZ, mask=lambda y,z: y**2+z**2>1)
plot_flow(y, z, np.linalg.solve(A, ny), inside, wave_source, (K,), size=(8,6), ymax=0.5)
The functions above only compute the flow and forces around a single shape - they take in one pair of coordinate arrays x,y
. How can we compute the flow induced by more than one shape, like a two boats in proximity or a multi-hull? Here are some ideas:
No. If we solve for the flows separately, we will miss the influence of the bodies on each other. When we add them up, the body BC won't be satisfied and so the flow and forces will be wrong.
Yes... in some cases. This is the classic approach (since it play nice with analytic models) but there are some caveats. Importantly, this only works when you have a plane of symmetry! So you can't use this to model two independent hulls - even if the geometries are identical, their motions won't be. It would also require using a new Green's function.
Yes... if we are careful. This approach includes the influence of each body on the other, but it will also create fake panels linking the bodies together. We need to zero the influence of those "non-panels". This is the approach I've used below.
def circle_cat(N,gap):
ya, za = ellipse(N//2, theta1=0)
yb, zb = ellipse(N//2, theta1=0)
return np.append(ya+gap, yb-gap), np.append(za, zb), [N//2]
def concat_sway(y,z,blank,G,args):
yc, zc, ny, _, h = properties(y[:-1], z[:-1], y[1:], z[1:])
A, ny, _ = construct_A(y, z, G, args)
# *zero* the connecting panels before solving
for i in blank:
A[i,:] = np.zeros(len(ny)); A[i,i] = 1; ny[i] = 0
qy = np.linalg.solve(A, ny) # sway source
B = [h*potential(yc, zc, y[i], z[i], y[i+1], z[i+1], G, args) for i in range(len(yc))]
myy = -ny@B@qy # sway added-mass
return qy, myy
N, K, gap = 32, 2, 3
y, z, blank = circle_cat(N, gap)
catYZ = np.meshgrid(np.linspace(-6,6,108), np.linspace(-2,0,18))
catYZ = mask_grid(*catYZ, mask=lambda y,z: (abs(y)-gap)**2 + z**2<1)
qy, myy = concat_sway(y, z, blank,wave_source, (K,))
print(f'Sway added mass = {myy:.3g}')
plot_flow(y, z, qy, catYZ, wave_source, (K,), size=(14,6), ymax=0.5)
Sway added mass = 1.12+0.579j
Look at the wave in the center of the plot, between the two circles. Since the arrows are all black instead of green, this is a standing wave caught between the hulls. As with the internal wave mentioned above, this external standing wave will also produce large forces and we therefore need to be careful to consider it's impact on multi-hull performance.
We will run a geometric parameter sweep on the distance between the centers $Y$=2gap
to investigate the impact of the geometry and on the forces.
@np.vectorize
def sway(omega, gap, rhoV=np.pi):
y, z, blank = circle_cat(32, gap)
_ , ma = concat_sway(y, z, blank, wave_source, (omega**2,))
return np.real(ma/rhoV)
def myplot(x,y,label,**kwargs): # remove segment connecting +/- singularities
i = np.argmax(np.diff(y) < -4) + 1
plt.plot(x[:i], y[:i], label=label, **kwargs)
plt.plot(x[i:], y[i:], **kwargs)
omega = np.linspace(0.1, 1.3, 21)
plt.figure(figsize=(4, 5))
for Y in range(3, 7):
myplot(omega, sway(omega, Y/2), label=Y, color=plt.cm.Greens((Y-2)/4))
plt.legend(title=r'$Y/R$')
plt.axhline(y=0, c='k')
plt.xlabel(r'$\omega\sqrt{R/g}$')
plt.ylabel(r'$\frac{A_{22}}{\rho V}$', rotation=0, size=14)
plt.ylim(-1.5, 3.5)
plt.xlim(0, 1.8);
This plot (formatted to match the figure in Faltisen) shows that the resonant standing wave between the hulls causes a singularity in the linear added-mass. As before, the singularity indicates a breakdown of linear theory but there will still be very large (positive or negative) added mass near this frequency.
Note: The resonant frequency decreases as the gap is increased, corresponding to the longer wavelength of the standing wave. Building a cat with a spacing similar to the dominant water wavelength in a region would lead to serious sea keeping and stability issues.
That's it for the second notebook. Let's review the main takeaway points: