# Roots of the depth – specific energy function¶

In free surface flows, the specific energy is a fundamental quantity governing the physics of relevant phenomena, e.g, Chow (2009). Recently, Computational Community (e.g Le Floch & Duc Thanh (2011), Castro et al. (2013), Bernetti et al. (2008)) has re-discovered its crucial role in solving some Shallow Water problems on discontinuous bottom. In the work of Valiani & Caleffi (2008) the analytical inversion of the depth - specific energy relationship is given.

The reference:

Valiani, A. & Caleffi, V., Depth-energy and depth-force relationships in open channel flows: Analytical findings (2008) Advances in Water Resources, 31 (3), pp. 447-454.

In [1]:
# preparation of the numerical environment
%matplotlib inline
from numpy import sqrt, cos, pi, arctan, linspace
from matplotlib import pyplot as plt


## The problem to solve¶

For a fixed value of the specific discharge $q$, the classical expressions of the specific energy $E$ as functions of the depth $Y$ is:

$$E = Y + \frac{{q^2 }}{{2g\,Y^2 }};$$

where $g$ the is gravity acceleration, and can be interpreted as a third-order equation in $Y$.

The problem to solve is to analytically found the two pysically meaningful depths $Y_{sb}$ and $Y_{sp}$ (subcritical and supercritical, respectivelly) for given values of $E=E_0$ and $q=q_0$.

## The analytical solution¶

The analytical solution described in Valiani & Caleffi (2008) is given by the following steps:

1. for the given discharge $q_0$, the critical values of the specific energy $E_c$ and depth $Y_c$ are: $$Y_c = \sqrt[3]{{\frac{{q_0^2 }}{g}}};\quad E_c = \frac{3}{2}Y_c;$$

2. the non-domensional depth, $\eta$, and the non-domensional energy, $\Gamma$, are obtained using their critical values: $$\eta = \frac{Y}{{Y_c }};\quad \Gamma = \frac{E}{{E_c }};$$

3. the non-domensional form of the depth - specific energy relationship becomes: $$\Gamma = \frac{2}{3}\eta + \frac{1}{3}\left( {\frac{1}{\eta }} \right)^2;$$

4. this equation has two meaningful solutions (see Valiani & Caleffi (2008) for the details): $$\eta _{sb} = \frac{{\Gamma _0 }}{2}\left[ {1 + 2\cos \left( {\frac{\pi -2 \alpha}{3}} \right)} \right]; \quad \text{subcritical depth}$$ $$\eta _{sp} = \frac{{\Gamma _0 }}{2}\;\left[ {1 + 2\cos \left( {\frac{\pi+2\alpha}{3}} \right)} \right]; \quad \text{supercritical depth}$$ where $\alpha = \arctan \left( {\sqrt {\Gamma _0^3 - 1} } \right)$;

5. corresponding to the dimensional form: $$Y_{sb} = \eta _{sb}Y_c; \qquad Y_{sp} = \eta _{sp}Y_c.$$

The following python code is proposed to highlight the simplicity of the method.

In [2]:
def radix(q_0,E_0,g):
"""
Compute analitically the roots (Y_sb, Y_sp) of the equation:
E_0 - Y - q_0^2/(2 g Y^2) = 0

Parameters
----------
Input:                             |  Output:
q_0 : float                    |     Y_sb : float
[m^2/2] specific discarge    |       [m] subcritical depth
E_0 : float                    |     Y_sp : float
[m] specific energy          |       [m] supercritical depth
g : float                      |
[m/s^2] gravity              |
"""

Y_c = (q_0**2/g)**(1.0/3.0)   # critical depth
E_c = 3.0/2.0*Y_c             # crtitical energy

# check the esistence of physical solutions!
if E_0 > E_c:
Gamma_0 = E_0/E_c
else:
print('no physical solutions!')
return 0.0, 0.0

alpha = arctan(sqrt(Gamma_0**3 - 1.0))

# non-dimensional subcritical depth
eta_sb = Gamma_0 / 2.0 * (1.0 + 2.0*cos(( pi - 2.0*alpha)/3.0))
# non-dimensional supercritical depth
eta_sp = Gamma_0 / 2.0 * (1.0 + 2.0*cos(( pi + 2.0*alpha)/3.0))

Y_sb = eta_sb * Y_c # subcritical depth
Y_sp = eta_sp * Y_c # supercritical depth

return Y_sb, Y_sp


## Test of the solution¶

To test the solution we find the two physical roots of the specific energy function for $q$ = 2.00 [m$^2$/2] and $E$ = 2.50 [m]:

In [5]:
gg = 9.81    # [m/s^2] gravity
qq = 2.00    # [m^2/2] specific discarge
EE = 2.50    # [m]     specific energy

Y = linspace(0.2,3,100)
E = Y + qq**2 / (2.0*gg*Y**2)

plt.figure(2,figsize=(10, 7))
plt.title('Speficic Energy vs. Depth')
plt.plot(E,Y,label='E = E(Y)')
plt.plot((EE,EE),Y_ss,'or',label='solutions')
plt.xlabel('E [m]')
plt.ylabel('Y [m]')
plt.axis((1.0,4.0,0.0,3.0))
plt.grid('on')
plt.legend()
plt.show()


Bibliography:

1. Bernetti, R., Titarev, V.A. & Toro, E.F. Exact solution of the Riemann problem for the shallow water equations with discontinuous bottom geometry (2008) Journal of Computational Physics (227), pp. 3212–3243.
2. Castro Díaz, M.J., López-García, J.A. & Parés C. High order exactly well-balanced numerical methods for shallow water systems (2013) Journal of Computational Physics, 246, pp. 242–264.
3. Chow, V.T., Open-Channel Hydraulics, 2009, Blackburn Press.
4. LeFloch P.G. & Duc Thanh M., A Godunov-type method for the shallow water equations with discontinuous topography in the resonant regime (2011) Journal of Computational Physics, 230, pp. 7631–7660.
5. Valiani, A. & Caleffi, V., Depth-energy and depth-force relationships in open channel flows: Analytical findings (2008) Advances in Water Resources, 31, pp. 447-454.