Altura e direção de onda usando a lei de snell e batimetria realística

Esse bloco apenas lê os dados da planilha do Excel, convert a batimetria para metros e inverte o sinal (o modelo entra com "h" positivo).



In [ ]:

    
from pandas import read_excel
xls = -1e-2 * read_excel('./batimetria.xlsx', sheetname='Plan1', skiprows=1, index_col=0, parse_cols=(0, 1, 4, 7))
xls.index.name = 'Distance [m]'
xls.columns = ['T1', 'T2', 'T3']
xls



In [ ]:

    
ax = xls.plot()
ax.invert_yaxis()
ax.set_ylabel('Depth [m]')

O bloco abaixo é o modelo. Temos que entrar com,

T: Período de onda típico;
h: Batimetria que vocês mediram;
Ho: Altura de onda inicial;
thetao: Angulo de incidência das ondas.



In [ ]:

    
import numpy as np
import numpy.ma as ma

# Constantes.
g = 9.81
twopi = 2 * np.pi
    
# Relação de dispersão (Resolvida pelo método de Newton-Raphson).
def dispersion_relationship(T, h, Ho, thetao):
    T, h, Ho, thetao = np.broadcast_arrays(T, h, Ho, thetao)
    omega = twopi / T
    Lo = (g * T ** 2) / twopi
    # Começa pela relação de dispersão de águas profundas.
    k = omega / np.sqrt(g)
    # Vamos aproximando por incrementos de "f".
    f = g * k * np.tanh(k * h) - omega ** 2
    while np.abs(f.max()) > 1e-7:
        dfdk = (g * k * h * (1 / (np.cosh(k * h))) ** 2 + g * np.tanh(k * h))
        k = k - f / dfdk
        f = g * k * np.tanh(k * h) - omega ** 2

    # Com o número de onda resolvido podemos calcular:
    L = twopi / k
    C = omega / k
    Co = Lo / T
    G = 2 * k * h / np.sinh(2 * k * h)
    
    # Aqui é basicamente a lei de Snell e o shoaling.
    theta = np.rad2deg(np.arcsin(C / Co * np.sin(np.deg2rad(thetao))))
    Kr = np.sqrt(np.cos(np.deg2rad(thetao)) / np.cos(np.deg2rad(theta)))
    Ks = np.sqrt(1 / (1 + G) / np.tanh(k * h))
    H = Ho * Ks * Kr
    # Retorna Altura de onda e ângulo.
    return map(ma.asanyarray, (theta, H))



In [ ]:

    
T = 5.  # Segundos
Ho = 0.1  # Metros.
thetao = 0  # Graus.

# Start at the deepest bathymetry and goes from there up to the coast.
h = np.fliplr(np.flipud(xls.values)[3:,:])
rows, cols = h.shape

theta, H = [], []
for row in range(rows):
    thetao, Ho = dispersion_relationship(T=T, h=h[row, :], Ho=Ho, thetao=thetao)
    H.append(Ho)
    theta.append(thetao)



In [ ]:

    
# Flip then back.
H = np.array(H)
theta = np.array(theta)



In [ ]:

    
import matplotlib.pyplot as plt

x = np.arange(3)
y = xls.index.values[:-3]

fig, (ax0, ax1, ax2) = plt.subplots(ncols=3, sharex=True, sharey=True, figsize=(12, 4))
cs0 = ax0.pcolorfast(x, y, h)
fig.colorbar(cs0, ax=ax0, extend='both')
ax0.set_title('Batimetria [m]')
ax0.set_xlabel('Perfil [sem unidade]')

cs1 = ax1.pcolorfast(x, y, H, vmax=1., vmin=0.1)
fig.colorbar(cs1, ax=ax1, extend='both')
ax1.set_ylabel(u'Distância [m]')
ax1.set_title('Altura de onda [m]')

cs2 = ax2.pcolorfast(x, y, theta, vmax=20, vmin=0)
fig.colorbar(cs2, ax=ax2, extend='both')
ax2.set_title(u'Ângulo [graus]')



In [ ]:

    
from IPython.display import Image
Image('./images/perfis.png', retina=True)



In [ ]:

    
ma.masked_invalid(H).max()



In [ ]: