Goal of this notebook:
We have a solution for the over-pressure due to the surface wave, which is valid for the full domain under the hypothesis that the depth of the propagation domain is very large compared to the wavelength of the wave:
$$ p(x, z, t) = \Re \left(\exp(kz) \exp(i(kx - \omega t) \right) $$If we choose a wavenumber, then we can deduce the frequency from the dispersion relation: $$ D(\omega, k) = \omega^2 - gk = 0 \Leftrightarrow \omega = \sqrt{ gk } $$
Given these equation, we can model our first case study: a sinusoidal pressure wave proapgating from left to right.
Let's first plot the pressure as a function of space.
import numpy as np
import holoviews as hv
hv.extension('matplotlib')
qm_opts = hv.opts.QuadMesh(fig_size=300, aspect=2.5, cmap='seismic')
x = np.linspace(-6.25, 6.25, num=400).reshape(1, -1)
z = np.linspace(0, -5, num=100).reshape(-1, 1)
X, Z = np.meshgrid(x, z)
k = 2 * np.pi / 4.
rho = 1000
def omega_of_k(k):
return np.sqrt(9.81 * k)
p = np.real(np.exp(k * z) * np.exp(1j * (k * x - omega_of_k(k) * 0)))
hv.QuadMesh((X, Z, p)).opts(qm_opts)
We can easily vary the value of the wavenumber in this plot to see different behaviours.
def p_of_k(k):
p = np.real(np.exp(k * z) * np.exp(1j * (k * x - omega_of_k(k) * 0)))
return hv.QuadMesh((X, Z, p)).opts(qm_opts)
dmap_p = hv.DynamicMap(p_of_k, kdims=['k']).redim.range(k=(1, 10))
(dmap_p[1.] + dmap_p[6]).cols(1).opts(fig_inches=7)
As seen above, low wavenumbers (large wavelengths) go deep, while large wavenumbers (small wavelengths) stay at the surface.
Of course, this pressure wave propagates from left to right due to the change of time. Let's do an animation of the pressure field as a function of its wavenumber.
def animate_pressure(k, Nframes=10):
omega = omega_of_k(k)
T = 2 * np.pi / omega
ts = np.arange(Nframes) / Nframes * T
ps = [np.real(np.exp(k * z) * np.exp(1j * (k * x - omega_of_k(k) * t))) for t in ts]
return {t: hv.QuadMesh((X, Z, p)).opts(qm_opts) for p,t in zip(ps, ts)}
%%output holomap='scrubber' fig='auto'
hmap_p_lowk = hv.HoloMap(animate_pressure(1))
hmap_p_lowk
Now the question is: what are the two displacement components due to the pressure field? We already know that the surface displacement is proportional to the value of the pressure on the surface.
def surface_displacement(k):
return hv.Curve((x, np.real(np.exp(1j * (k * x - omega_of_k(k) * 0)))))
dmap_y = hv.DynamicMap(surface_displacement, kdims=['k']).redim.range(k=dmap_p.range('k'))
(dmap_p * dmap_y)[k]
((dmap_p * dmap_y)[1] + (dmap_p * dmap_y)[6]).cols(1).opts(fig_inches=7)
Now this explains the displacement at the surface. But what about inside the fluid? It turns out that we can use the momentum balance equations with the known pressure to derive the displacement components and simple "Fourier maths":
def compute_solution(k, t=0, nx=28, nz=14):
x = np.linspace(-6.25, 6.25, num=nx).reshape(1, -1)
z = np.linspace(0, -5, num=nz).reshape(-1, 1)
X, Z = np.meshgrid(x, z)
p = np.exp(k * Z) * np.exp(1j * (k * X - omega_of_k(k) * t))
U = 1/(rho * omega_of_k(k) ** 2) * k * p
V = 1/(rho * omega_of_k(k) ** 2) * 1j * k * p
return X, Z, np.real(p), np.real(U), np.real(V)
X, Z, p, U, V = compute_solution(k, t=0)
We can turn this computation into a vector field quite easily.
vector_opts = hv.opts.VectorField(color='Magnitude',
magnitude=hv.dim('Magnitude').norm(),
rescale_lengths=False,
fig_size=300, aspect=2.5)
def make_vector_field(X, Z, U, V, vector_opts=None):
mag = np.sqrt(U**2 + V**2)
angle = (np.pi/2.) - np.arctan2(U/mag, V/mag)
field = hv.VectorField((X, Z, angle, mag))
if vector_opts is not None:
field = field.opts(vector_opts)
return field
make_vector_field(X, Z, U, V, vector_opts)
And we can even animate this:
def animate_vector_field(k, Nframes=10):
omega = omega_of_k(k)
T = 2 * np.pi / omega
ts = np.arange(Nframes) / Nframes * T
frames = {}
for t in ts:
X, Z, p, U, V = compute_solution(k, t)
frames[t] = make_vector_field(X, Z, U, V, vector_opts)
return frames
How does this look at low wavenumbers?
%%output holomap='scrubber'
hv.HoloMap(animate_vector_field(.2, Nframes=20))
Something that is not completely apparent in the above animation is that the trajectories of the particles are little circles. Let's try to superpose the points obtained over several frames on a still image.
def compute_trajectory_points(k, Nframes=10, viz_amp=3000, **kwargs):
omega = omega_of_k(k)
T = 2 * np.pi / omega
ts = np.arange(Nframes) / Nframes * T
positions = []
for t in ts:
X, Z, p, U, V = compute_solution(k, t, **kwargs)
positions.append(np.c_[(X + U * viz_amp).flat, (Z + V * viz_amp).flat])
return np.concatenate(positions)
point_opts = hv.opts.Points(color='k', fig_size=300, aspect=2.5, s=np.array([2]).reshape(-1, 1))
(hv.Points(compute_trajectory_points(0.1, nx=14, nz=7, Nframes=50)).opts(point_opts) + \
hv.Points(compute_trajectory_points(1., nx=14, nz=7, Nframes=50)).opts(point_opts)).cols(1).opts(fig_inches=7)
What we see here (for a low wavenumber wave) are the particle displacements plotted over one period of time, and they're all little circles. This does not change for high wavenumbers, however the circle amplitude gets smaller and smaller due to depth.
Let's now put everything together and do an animation of the wave, its particles and the circles:
def animate_everything(k, Nframes=20):
trajectory = hv.Points(compute_trajectory_points(k, Nframes=50)).opts(point_opts)
omega = omega_of_k(k)
T = 2 * np.pi / omega
ts = np.arange(Nframes) / Nframes * T
frames = {}
for t in ts:
X, Z, p, U, V = compute_solution(k, t)
frames[t] = hv.QuadMesh((X, Z, p)).opts(qm_opts) * trajectory * make_vector_field(X, Z, U, V, vector_opts)
return frames
%%output holomap='scrubber'
hv.HoloMap(animate_everything(0.4))
%%output holomap='scrubber'
hv.HoloMap(animate_everything(1.))
Let's now conclude with a recreation of pictures to be found in Milton Van Dyke's Album of fluid motion on the chapter on surface waves. I copied three of these together to get the following image:
from IPython.display import Image
Image(filename='files/VanDykeSurfaceWaves.png')
We will use the previous code to visualize the effect of two waves of opposite directions propagating at the same time. We will trace the resulting pressure field, the vectors as well as the trajectories all on top of each other.
def compute_solution_two_waves(reflected_part, k, t=0, nx=28, nz=14):
x = np.linspace(-6.25, 6.25, num=nx).reshape(1, -1)
z = np.linspace(0, -5, num=nz).reshape(-1, 1)
X, Z = np.meshgrid(x, z)
p_to_left = np.exp(k * Z) * np.exp(1j * (k * X - omega_of_k(k) * t))
p_to_right = reflected_part * np.exp(k * Z) * np.exp(1j * (k * X + omega_of_k(k) * t))
p = p_to_left + p_to_right
U = 1/(rho * omega_of_k(k) ** 2) * k * p
V = 1/(rho * omega_of_k(k) ** 2) * 1j * k * p
return X, Z, np.real(p), np.real(U), np.real(V)
def compute_trajectory_points_two_waves(reflection_part, k, Nframes=10, viz_amp=3000, **kwargs):
omega = omega_of_k(k)
T = 2 * np.pi / omega
ts = np.arange(Nframes) / Nframes * T
positions = []
for t in ts:
X, Z, p, U, V = compute_solution_two_waves(reflection_part, k, t, **kwargs)
positions.append(np.c_[(X + U * viz_amp).flat, (Z + V * viz_amp).flat])
return np.concatenate(positions)
We can now do a plot of all this animated information together.
def animate_two_waves(k, reflected_part, Nframes=10, viz_amp=4000):
omega = omega_of_k(k)
trajectories = hv.Points(compute_trajectory_points_two_waves(reflected_part,
k,
Nframes=100,
viz_amp=viz_amp)).opts(point_opts).opts(alpha=.3)
T = 2 * np.pi / omega
ts = np.arange(Nframes) / Nframes * T
frames = {}
circle_indices = None
for t in ts:
X, Z, p, U, V = compute_solution_two_waves(reflected_part, k, t)
frames[t] = hv.QuadMesh((X, Z, p)).opts(cmap='plasma', alpha=.2) * \
make_vector_field(X, Z, U, V) * \
hv.Points(((X + U * viz_amp).flat,
(Z + V * viz_amp).flat)).opts(color='black',
s=np.array([20]).reshape(-1, 1)) * \
trajectories
return frames
%%output holomap='scrubber'
hv.HoloMap(animate_two_waves(k=.62, reflected_part=0, Nframes=15), kdims='t').opts(show_frame=False, show_title=False)
%%output holomap='scrubber'
hv.HoloMap(animate_two_waves(k=.62, reflected_part=.53, Nframes=15), kdims='t').opts(show_frame=False, show_title=False)
%%output holomap='scrubber'
hv.HoloMap(animate_two_waves(k=.62, reflected_part=.999, Nframes=15), kdims='t').opts(show_frame=False, show_title=False)
The comparison to the still pictures I mentioned above are interesting:
They're is always more to learn. Here, my principal source of learning was the MOOC Fundamentals of waves and vibrations, for free on Coursera. Feynman's chapter on waves is also fantastic (http://www.feynmanlectures.caltech.edu/I_51.html)
http://courses.washington.edu/mengr543/handouts/Album-Fluid-Motion-Van-Dyke.pdf