Calculate the dispersion of flexural gravity waves (FGWs)

fgws.py

'''
Code to calculate the dispersion relation of Flexural Gravity Waves
'''
import numpy as np
import matplotlib.pyplot as plt

g=9.81      # Gravity, m/s**2
hw=100      # Water depth, m
hi = 300    # Ice thickness, m
rhow=1000   # Water density, kg/m**3
rhoi = 910  # Ice density, kg/m**3
nu = 0.3    # Poisson ratio
E = 8e9     # Youngs modulus (of ice), Pa

D = E * hi**3 / (1-nu) / 12         # Flexural rigidity
L = np.logspace(1,3,1000) * hi   # Range of wavelengths as short as ~10*hi
k = 2*np.pi / L                     # Wavenumber vector

# The flexural gravity dispersion relation (see Lipovsky, 2018 and references therein)
gamma = 1 / (k*np.tanh(k*hw)) # Part of the relation due to gravity waves
omega = np.sqrt( (D *k**4 + rhow * g)/(rhoi*hi + rhow*gamma) ) # Angular frequency (rad/s)
f = omega/(2*np.pi) # Ordinary frequency (cycles/s)
csw = np.sqrt(g*hw)

# Make plots
fig,ax=plt.subplots(1,2,figsize=(10,5))
fig.patch.set_facecolor('w')
plt.subplot(1,2,1)
plt.plot(f,L,label='fgw dispersion relation')
plt.plot((f[0],f[-1]),(hi,hi),'--',label='ice thickness')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Wavelength (m)')
plt.legend()
plt.xscale('log')
# plt.xlim([1e-3,1])
# plt.ylim([100,10000])
plt.yscale('log')
plt.grid('on')

plt.subplot(1,2,2)
plt.plot(f,omega/k,label='fgw dispersion relation')
plt.plot((f[0],f[-1]),(csw,csw),'--',label='Shallow water wave limit')
plt.grid('on')
# plt.yscale('log')
plt.xscale('log')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Wavespeed (m/s)')
plt.legend()

plt.tight_layout()
plt.show()