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_ss = radix(qq,EE,gg)
    
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.