Function to convert from Te to Tp for spectrum creation

[mbruggs]

Describe the feature:

There are scenarios in which we have an energy period (Te) and would like to create, for example, a JONSWAP spectrum. The functions for creating the spectrum use Tp rather Te, meaning a means of obtaining Tp for the target spectrum is needed.

DNV 205 (link) has approximates for the JONSWAP spectral moments in section 3.5.5.6 in terms of Hs and Tp

$$ M_{-1} = \frac{1}{16} H_s^2 \omega_p^{-1} \frac{4.2 + \gamma}{5 + \gamma} $$

$$ M_0 = \frac{1}{16} H_s^2 $$

We also have

$$ Te = 2\pi\frac{M{-1}}{M_0} $$

Combining the above, we then have

$$ T_e = T_p \frac{4.2 + \gamma}{5 + \gamma} $$

which we can use to get Tp from a Te.

Is this something you might be interesting in having in MHKit? I'm happy to contribute the code but wanted to gauge interest before opening a PR.

[ssolson]

I think this sounds useful. @akeeste & @cmichelenstrofer any thoughts?

[cmichelenstrofer]

I find these approximations extremely useful. But I usually use the ITTC's approximation:

See full document: ITTC2002.pdf - ITTC2002.pdf (Table A4, page 547)

I wonder how different these two are and which one should we choose. @ryancoe any thoughts?

[ryancoe]

That is very useful. I've always used that ITTC document as well...

[akeeste]

I agree that a conversion would be beneficial to add into MHKiT, though I'm not familiar with differences in the formulations. Is there an IEC definition that covers this?

[mbruggs]

Thanks for the input everyone. The plot below compares the results from calculating Te using the spectral moments and the two approximations in the thread. It looks like the approximations are very similar, so I'll implement the ITTC as that what people seem to be using anyway.

import numpy as np
import matplotlib.pyplot as plt
from mhkit.wave.resource import jonswap_spectrum, energy_period

Tp = 10.0
Hs = 5.0

f = np.linspace(1 / (10 * Tp), 3/ Tp, 100)
gamma = np.linspace(1,7, 10)

Te = []
for g in gamma:
    S = jonswap_spectrum(f, Tp, Hs, g)
    Te.append(energy_period(S).values[0].squeeze())

Te_DNV = Tp * (4.2 + gamma) / (5.0 + gamma)
Te_ITTC = Tp * (0.8255 + 0.03852*gamma - 0.005537*gamma**2 + 0.0003154*gamma**3)

plt.figure()
plt.plot(gamma, Te)
plt.plot(gamma, Te_DNV)
plt.plot(gamma, Te_ITTC)
plt.xlabel(r'$\gamma$')
plt.ylabel('Te [s]')
plt.legend([r'$M_{-1} / M_0$', 'DNV205', 'ITTC'])
plt.title(f"Tp = {Tp}")
plt.grid()

[ryancoe]

Well done Mark. Agree with your decision.

[ssolson]

Added in #210