Contents

Potential equation

Una de las posibles soluciones del potencial \(\Phi\) es el de una onda propagándose en el sentido postivo del eje \(x\), que se expresa como:

\[ \Phi(x,z,t)=-\frac{Ag}{\omega}\frac{\cosh k(h+z)}{\cosh kh}\sin(kx-\omega t) \]

Pero las soluciones correspondientes a ondas suelen expresarse también haciendo uso de la variable compleja (facilitando esta el álgebra considerablemente), por lo que la expresión se reduce a:

\[ \Phi(x,z,t) = \textrm{Re} \left\{ -\frac{ig}{\omega}A\frac{\cosh k(h+z)}{\cosh kh}\exp{-i(kx-\omega t)}\right\} \]

dónde Re implica la parte real de la expresión entre corchetes.

Importamos las librerías necesarias

# import maths
import numpy as np
from sympy import *

# import matplotlib / plotting
import matplotlib
matplotlib.rcParams['animation.embed_limit'] = 2**32
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
plt.style.use('dark_background')
from IPython.display import HTML # diplay anim

import warnings
warnings.filterwarnings("ignore")

cargamos las variables para SymPy

# save variables for SymPy
A, g, w, k, h, x, z, t = symbols('A g w k h x z t')

y definimos las ecuaciones

# potential equation
pot = -(A*g/w) * (cosh(k*(h+z))/cosh(k*h)) * sin(k*x-w*t)

pot
\[\displaystyle - \frac{A g \sin{\left(k x - t w \right)} \cosh{\left(k \left(h + z\right) \right)}}{w \cosh{\left(h k \right)}}\]
# calculate eta, u, w, ax y az
eta = 1/g * diff(pot,t).evalf(subs={'z':0})
u = - diff(pot, x)
w = - diff(pot, z)
ax = - diff(u, t)
az = - diff(w, t)

eta
\[\displaystyle A \cos{\left(k x - t w \right)}\]
# evaluate u, w and eta
def evaluate_potential(xs,zs,tt,wave_parameters):
    us = []
    ws = []
    etas = []
    u_lamb = lambdify((x,z,t),u.subs(wave_parameters),'math')
    w_lamb = lambdify((x,z,t),w.subs(wave_parameters),'math')
    eta_lamb = lambdify((x,z,t),eta.subs(wave_parameters),'math')
    for xss in xs:
        for zss in zs:
            us.append(u_lamb(xss,zss,tt))
            # us.append(float(u.evalf(subs={
            #     **{'x':xss,'z':zss,'t':tt},**wave_parameters
            # })))
            ws.append(w_lamb(xss,zss,tt))
            # ws.append(float(w.evalf(subs={
            #     **{'x':xss,'z':zss,'t':tt},**wave_parameters
            # })))
        etas.append(eta_lamb(xss,zss,tt))
        # etas.append(float(eta.evalf(subs={
        #     **{'x':xss,'z':zss,'t':tt},**wave_parameters
        # })))
        
    return np.array(us), np.array(ws), np.array(etas)

Evaluate and plot results

# fixed parameters for the wave
A_value = 1.5 # meters
g_value = 9.8 # m / s^2
T_value = 6 # seconds
h_value = 10 # mean depth (meters)
# derived parameters for the wave
w_value = 2*np.pi/T_value
k_value = w_value**2/g_value # wave number

wave_parameters = {
    'A': A_value, 
    'g': g_value, 
    'w': w_value, 
    'k': k_value, 
    'h': h_value
}

# plot parameters
fig = plt.figure(figsize=(12,6))
x_range = (0,100,24)
x_values = np.linspace(*x_range)
y_range = (-10,4,10)
y_values = np.linspace(*y_range)
ax = plt.axes(xlim=(x_range[0]-5,x_range[1]+5), 
              ylim=(y_range[0]-0.5,y_range[1]+0.5))
line, = ax.plot([], [], lw=4)
u0, w0, eta0 = evaluate_potential(x_values,y_values,0,wave_parameters)
line.set_data(x_values, eta0)
idxs = []
for i in range(len(eta0)-1):
    idxs += list(
        np.where(y_values<eta0[i])[0]+(i*len(y_values))
    ) # get idxs below wave
quiv = ax.quiver(
    np.meshgrid(y_values,x_values)[1].reshape(-1)[idxs],
    np.meshgrid(y_values,x_values)[0].reshape(-1)[idxs],
    u0[idxs], w0[idxs], u0[idxs], width=0.002, cmap='bwr_r'
)

def init():
    line.set_data([],[])
    
    return line, quiv

def animate(tt):
    global quiv
    quiv.remove()
    us, ws, etas = evaluate_potential(x_values,y_values,tt/5,wave_parameters)
    line.set_data(x_values, etas)
    idxs = []
    for i in range(len(etas)-1):
        idxs += list(
            np.where(y_values<etas[i])[0]+(i*len(y_values))
        ) # get idxs below wave
    quiv = ax.quiver(
        np.meshgrid(y_values,x_values)[1].reshape(-1)[idxs],
        np.meshgrid(y_values,x_values)[0].reshape(-1)[idxs],
        us[idxs], ws[idxs], us[idxs], width=0.002, cmap='bwr_r'
    )

    return line, quiv
plt.close()

# we create the animation, based in the previous cell
anim = FuncAnimation(fig, animate, init_func=init, frames=20, interval=200, blit=True)

# visualize animation
HTML(anim.to_jshtml())

# anim.save('potencial.gif')