TESIS DE MAESTRÍA EN HIDRÁULICA COMPUTACIONAL
UNIVERSIDAD NACIONAL DE MAYOR DE SAN MARCOS

  1. Solución Numérica de la ecuación de Manning

    • Valor incial
      •

      Para calcular el valor incial del tirante usamos la hipótesis de canal ancho. La ecuación de Manning viene dada por:

      \[Q=\frac{A}{n}R^{2/3}S_{o}^{1/2}\qquad[1.1]\]

      donde \(Q\): caudal, \(A\): área hidráulica, \(S_{o}\): pendiente de fondo, \(n\): rugosidad. Si \(b>>y\) (hipótesis de canal ancho) entonces \(R→y\), por lo tanto

      \[Q=\frac{A}{n}y_{o}^{2/3}S_{o}^{1/2}\qquad[1.2]\]
      \[\frac{Qn}{S_{o}^{1/2}}=y_{o}^{5/3}b\qquad[1.3]\]
      \[Z=y_{o}^{5/3}b\qquad[1.4]\]

      donde \(Z\) es el factor de sección, despejando \(y\) tenemos

      \[y_{o}=\left(\frac{Z}{b}\right)^{3/5}\qquad[1.5]\]

      el valor \(y_{o}\) se usará como aproximación en el método de Newton-Raphson para encontrar el valor \(y\) en la Ec. [1.1]


    • Método de Newton-Raphson
      •

      Reestructuramos la Ec. [1.1] a la forma \(f(y)=0\) y hallamos su derivada

      \[ \begin{aligned} &Q=\frac{A}{n}\left(\frac{A}{P}\right)^{2/3}S_{o}^{1/2} \\ &\frac{Qn}{S_{o}^{1/2}}=\left(\frac{A}{P}\right)^{2/3} \\ &Z=\left(\frac{A}{P}\right)^{2/3} \end{aligned}\]
      \[f(y)=\frac{A^{5/3}}{P^{2/3}}-Z \qquad[1.6]\]
      \[f'(y)=\frac{d}{dy}\left[A^{5/3}P^{-2/3}\right]\]
      \[f'(y)=A^{5/3}\left(\frac{-2}{3}\right)P^{-5/3}\frac{dP}{dy}+P^{-2/3}\left(\frac{5}{3}\right)A^{2/3}\frac{dA}{dy}\]

      para una sección trapezoidal se sabe que, \(A=y(b+zy)\) y \(P=b+2y\sqrt{z^{2}+1}\), y su dervidas son \(dA/dy=b+2zy\) y \(dP/dy=2\sqrt{z^{2}+1}\), respectivamente. Reordenando y reeemplazando expresiones obtenemos

      \[f'(y)=\frac{5}{3}(b+2zy)R^{2/3}-\frac{4}{3}(z^{2}+1)^{1/2}R^{5/3} \qquad[1.7]\]

      donde \(R\): radio hidráulico, \(z\): talud. El cálculo de \(y\) se realiza a partir de un tirante inicial \(y_{o}\) y luego se itera en la siguiente expresión

      \[ y_{n+1}=y_{n}-\frac{f(y)}{f'(y)} \qquad[1.8]\]

      donde \(y_{n+1}\) = tirante en la iteración \(n+1\), \(y_{n}\): tirante en la iteración \(n\), cuando \(n = 0\) entonces \(y = y_{o}\).


  2. Ecuaciones de Saint Venant


    3. Ecuaciones de Conservación

    Ecuacion de la conservación de la masa

    \[ \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}=0 \qquad[2.1] \]

    Ecuacion de la conservación del momento lineal

    \[ \frac{1}{g} \frac{\partial V}{\partial t} + \frac{V}{g} \frac{\partial V}{\partial x} + \frac{\partial y}{\partial x}+S_{f}-S_{o}=0 \]
    \[ gA\left\{ \frac{1}{g} \frac{\partial V}{\partial t} + \frac{V}{g} \frac{\partial V}{\partial x} + \frac{\partial y}{\partial x}+S_{f}-S_{o}=0 \right\} \]
    \[ A\frac{\partial V}{\partial t} + VA\frac{\partial V}{\partial x} + gA\frac{\partial y}{\partial x}+gA\left(S_{f}-S_{o}\right)=0 \]
    \[ A\frac{\partial V}{\partial t} + VA\frac{\partial V}{\partial x} + g(yb+zy^{2})\frac{\partial y}{\partial x}+gA\left(S_{f}-S_{o}\right)=0 \]
    3.

    Introducimos las variables \(V\) y \(y\) dentro de las derivadas parciales, de tal forma que, al aplicar la regla del producto, se recuperen los mismos términos externos. Esto aplica únicamente para el segundo y tercer término.

    \[ \frac{\partial}{\partial t} \left(VA\right) + \frac{\partial}{\partial x} \left(\frac{V^{2}A}{2}\right) + \frac{\partial}{\partial x} \left( \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right) + gA\left(S_{f}-S_{o}\right)=0 \]
    \[ \frac{\partial Q}{\partial t} + \frac{\partial}{\partial x} \left(\frac{Q^{2}}{2A}\right) + \frac{\partial}{\partial x} \left( \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right) + gA\left(S_{f}-S_{o}\right)=0 \]
    \[ \frac{\partial Q}{\partial t} + \frac{\partial}{\partial x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right) + gA\left(S_{f}-S_{o}\right)=0 \qquad[2.2] \]

    Reescribimos las ecuaciones de movimiento en forma vectorial, de la siguiente manera

    \[ \frac{\partial }{\partial t} \begin{bmatrix} A \\ Q\end{bmatrix} + \frac{\partial }{\partial x} \begin{bmatrix} Q \\ \frac{Q^{2}}{2A} +\frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \end{bmatrix} = \begin{bmatrix} 0 \\ gA\left(S_{o}-S_{f}\right)\end{bmatrix} \]
    \[ \frac{\partial \hat{U}}{\partial t}+\frac{\partial \hat{F}}{\partial x}=\hat{S} \qquad[2.3]\]

    La Ec. [2.3] nos servirá más adelante para aplicar los distintos esquemas de discretización.


    Método de diferencias finitas (FDM)


    Esquema Difuso o LAX


    El esquema de Difuso es un método explícito de primer orden de precisión en el tiempo y segundo orden en el espacio, que se emplea principalmente como método básico para la discretización de las ecuaciones de conservación. Los operadores diferenciales para el esquema de LAX son

    \[ \begin{aligned} &\frac{\partial \hat{F}}{\partial x}=\frac{\hat{F}^{k}_{i+1}-\hat{F}^{k}_{i-1}}{2\Delta x} \\ &\frac{\partial \hat{U}}{\partial t}=\frac{\hat{U}^{k+1}_{i}-\hat{U}^{k}_{i}}{\Delta t} \\ &\hat{U}^{k}_{i}=\frac{\hat{U}^{k}_{i+1}+\hat{U}^{k}_{i-1}}{2} \\ &\hat{S}^{k}_{i}=\frac{\hat{S}^{k}_{i+1}+\hat{S}^{k}_{i-1}}{2} \end{aligned} \qquad[2.4] \]

    reemplazando los operadores diferenciales en la Ec. [2.3] se tiene

    \[ \frac{\hat{U}^{k+1}_{i}-\hat{U}^{k}_{i}}{\Delta t} + \frac{\hat{F}^{k}_{i+1}-\hat{F}^{k}_{i-1}}{2\Delta x} = \frac{\hat{S}^{k}_{i+1}+\hat{S}^{k}_{i-1}}{2} \]
    \[ \hat{U}^{k+1}_{i}=\frac{\hat{U}^{k}_{i+1}+\hat{U}^{k}_{i-1}}{2}-\frac{\Delta t}{2\Delta x} \left( \hat{F}^{k}_{i+1}-\hat{F}^{k}_{i-1} \right) + \frac{\Delta t}{2} \left(\hat{S}^{k}_{i+1}+\hat{S}^{k}_{i-1} \right)\qquad[2.5] \]

    Finalmente las ecuaciones de movimiento quedan discretizadas así:


    Ecuación de la conservación de la masa

    \[ A^{k+1}_{i}=\frac{A^{k}_{i+1}+A^{k}_{i-1}}{2}-\frac{\Delta t}{2\Delta x} \left( Q^{k}_{i+1}-Q^{k}_{i-1} \right) \qquad[2.6] \]

    Ecuación de la conservación del momentum

    \[ Q^{k+1}_{i}=\frac{Q^{k}_{i+1}+Q^{k}_{i-1}}{2}-\frac{\Delta t}{2\Delta x} \begin{Bmatrix} \left( \frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{i+1} \\ -\left( \frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{i-1} \end{Bmatrix} \]

    \[ + \frac{\Delta t}{2} \left[ \left( gA(S_{o}-S_{f}) \right)^{k}_{i+1} + \left( gA(S_{o}-S_{f}) \right)^{k}_{i-1} \right] \qquad[2.7] \]

    donde la pendiente de fricción \(S_{f}\) se evalua con la siguiente expresión

    \[ S_{f}=\frac{Q|Q|n^{2}P^{4/3}}{A^{10/3}} \qquad[2.8] \]

    para \(A=y(b+zy)\) y \(P=b+2y\sqrt{z^{2}+1}\). Y el centroide \(y_{c}\) con

    \[ y_{c}=\frac{4zy^{2}+3by}{6(b+zy)} \qquad[2.9] \]

    usamos la siguiente expresión para calcular el tirante a partir del área :

    \[ \begin{aligned} &A=y(b+zy) \\ &zy^{2}+by-A=0 \end{aligned} \]
    \[ y=\frac{-b+\sqrt{b^{2}+4zA}}{2z} \qquad[2.10] \]

    Esquema Mac Cormack


    El esquema de MacCormack es un método explícito de diferencias finitas, de segundo orden en el espacio y en el tiempo, ampliamente utilizado para la resolución numérica de las ecuaciones de conservación. Se fundamenta en un enfoque predictor–corrector o de doble barrido, en el cual la discretización de los términos espaciales se adapta al sentido de propagación de la onda: para ondas que viajan en sentido positivo (+), la solución se predice mediante derivadas hacia atrás y se corrige con derivadas hacia adelante, mientras que para ondas que se propagan en sentido negativo (–) el uso de los operadores diferenciales se invierte. Este procedimiento permite mejorar la precisión y capturar adecuadamente la dinámica de propagación de ondas y discontinuidades.


    Predictor

    Para el nivel de tiempo \(*\) usamos la derivada hacia atrás para todos los nodos

    \[ U^{*}_{i}= U^{k}_{i} - \frac{\Delta t}{\Delta x} \left( F^{k}_{i}-F^{k}_{i-1} \right) + \frac{\Delta t}{2} \left(S^{k}_{i+1}+S^{k}_{i-1} \right)\qquad[2.11] \]

    Corrector

    Para el nivel de tiempo \(**\) usamos la derivada hacia adelante y lo valores del nivel de tiempo \(*\), nuevamente para todos lo nodos

    \[ U^{**}_{i}= U^{k}_{i} - \frac{\Delta t}{\Delta x} \left( F^{*}_{i+1}-F^{*}_{i} \right) + \frac{\Delta t}{2} \left(S^{*}_{i+1}+S^{*}_{i-1} \right)\qquad[2.12] \]

    Luego hallamos las variables desconocidas en el nivel de tiempo \(k+1\) mediante un promedio del predictor y corrector

    \[ U^{k+1}_{i}= \frac{U^{*}_{i}+U^{**}_{i}}{2}\qquad[2.13] \]

  3. Condiciones de contorno

    • Aguas arriba
      •

      Usamos un hidrograma de entrada con la forma de una curva gaussiana y con los siguientes parámetros.

      \[ Q=Q_{min}+(Q_{max}-Q_{min})e^{-\frac{1}{2}{\left(\frac{t-αT_{h}}{βT_{h}}\right)}^2} \qquad[3.1]\]

      donde \(Q\): caudal en función del tiempo, \(Q_{min}\): caudal mínimo, \(Q_{max}\): caudal máximo, \(T_{h}\): tiempo de duración del hidrograma, y \(α,β\) son parámetros que controlan la tendencia central y la dispersión del hidrograma.


    • Aguas abajo
      •

      Restamos de la Ec. [2.1] a la Ec. [2.2] multiplicada por \(\sqrt{1/gb}\) para hacer compatible sus dimensiones

      \[ \frac{\partial A}{\partial t}+\frac{\partial Q}{\partial x}-\sqrt{\frac{1}{gb}} \begin{Bmatrix} \frac{\partial Q}{\partial t} \\ +\frac{\partial }{\partial x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right) \\ -gA\left(S_{o}-S_{f}\right) \end{Bmatrix} = 0 \qquad[3.2] \]

      Usamos los operadores diferenciales hacia adelante en espacio y tiempo para el último y penúltimo punto de eje X

      \[ \begin{aligned} &\frac{\partial f}{\partial t} = \frac{f^{k+1}_{N}-f^{k}_{N}}{\Delta t} \\ &\frac{\partial f}{\partial x} = \frac{f^{k}_{N}-f^{k}_{N-1}}{\Delta x} \end{aligned} \qquad[3.3] \]

      Reemplazamos los operadores en la Ec. [3.2]

      \[ \frac{A^{k+1}_{N}-A^{k}_{N}}{\Delta t} + \frac{Q^{k}_{N}-Q^{k}_{N-1}}{\Delta x} -\sqrt{\frac{1}{gb}} \begin{Bmatrix} \frac{Q^{k+1}_{N}-Q^{k}_{N}}{\Delta t} \\ + \frac{1}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N} \\ - \frac{1}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N-1} \\ - \left[ gA\left(S_{o}-S_{f}\right) \right]_{N}^{k} \end{Bmatrix} = 0 \qquad[3.4] \]

      La ecuación anterior presenta dos incógnitas \( \left( Q_{N}^{k+1}, A_{N}^{k+1} \right) \) al aplicar la condición de contorno de profundidad normal aguas abajo, se incorpora la ecuación de Manning para vincular ambas variables y reducir el sistema a una sola incógnita.

      \[ Q_{N}^{k+1} = \frac{S_{o}^{1/2}}{n}\left( \frac{A^{5/3}}{P^{2/3}} \right)_{N}^{k+1} \qquad[3.5] \]

      Ahora, expresamos \( P_{N}^{k+1} \) en función de \( A_{N}^{k+1} \) de la siguiente manera.

      Despejamos \(y\) de la ecuación del perímetro para una sección trapezoidal

      \[ P=b+2y\sqrt{z^{2}+1} \]
      \[ y=\frac{(P-b)}{2\sqrt{z^{2}+1}} \]
      \[ y=\phi(P-b) \]

      Reemplazamos \(y\) en la ecuación del área

      \[ \begin{aligned} &A = zy^{2}+by \\ &A = z\phi^{2}(P-b)^{2} + b\phi(P-b) \\ &A = z\phi^{2}(P^{2}-2Pb+b^{2}) + b\phi(P-b) \\ &A = z\phi^{2}P^{2} - 2z\phi^{2}Pb + z\phi^{2}b^{2} + b\phi P - b^{2}\phi \\ &(z\phi^{2})P^{2} + (b\phi-2z\phi^{2}b)P + (z\phi^{2}b^{2} - b^{2}\phi - A) = 0 \\ &\alpha P^{2} + \beta P + \gamma = 0 \end{aligned} \]

      usamos la fórmula cuadrática para encontrar el valor de \( P \)

      \[ P = \frac{-\beta \pm \sqrt{\beta^{2} - 4\alpha\gamma}}{2\alpha} \]

      Simplificando y tomando el signo \( (+) \), obtenemos

      \[ P = b \left( 1-\frac{\sqrt{z^{2}+1}}{z} \right) + \frac{\sqrt{z^{2}+1}}{z} \sqrt{b^{2}+4zA} \qquad[3.6] \]

      reemplazamos \( P \) en la Ec. [3.5]

      \[ Q_{N}^{k+1} = \frac{S_{o}^{1/2}}{n} \left[ \frac{A^{5/3}}{\left( b \left( 1-\frac{\sqrt{z^{2}+1}}{z} \right) + \frac{\sqrt{z^{2}+1}}{z} \sqrt{b^{2}+4zA} \right)^{2/3}} \right]_{N}^{k+1} \qquad[3.7] \]

      y por último reemplazamos \( Q_{N}^{k+1} \) en la Ec. [3.4]

      \[ \frac{A^{k+1}_{N}-A^{k}_{N}}{\Delta t} + \frac{Q^{k}_{N}-Q^{k}_{N-1}}{\Delta x} -\sqrt{\frac{1}{gb}} \begin{Bmatrix} \frac{Q^{k+1}_{N}}{\Delta t} -\frac{Q^{k}_{N}}{\Delta t} \\ + \frac{1}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N} \\ - \frac{1}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N-1} \\ - \left[ gA\left(S_{o}-S_{f} \right) \right]_{N}^{k} \end{Bmatrix} = 0 \]
      \[ \frac{A^{k+1}_{N}-A^{k}_{N}}{\Delta t} + \frac{Q^{k}_{N}-Q^{k}_{N-1}}{\Delta x} -\sqrt{\frac{1}{gb}} \begin{Bmatrix} \frac{1}{\Delta t} \frac{S_{o}^{1/2}}{n} \left[ \frac{A^{5/3}}{\left( b \left( 1-\frac{\sqrt{z^{2}+1}}{z} \right) + \frac{\sqrt{z^{2}+1}}{z} \sqrt{b^{2}+4zA} \right)^{2/3}} \right]_{N}^{k+1} -\frac{Q^{k}_{N}}{\Delta t} \\ + \frac{1}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N} \\ - \frac{1}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N-1} \\ - \left[ gA\left(S_{o}-S_{f} \right) \right]_{N}^{k} \end{Bmatrix} = 0 \qquad[3.8] \]

      La ecuación [3.8] es implícita para \(A_{N}^{k+1}\), procedemos a resolverla por el método de Newton.

      \[ A^{k+1}_{N}=A^{k}_{N} - \frac{\Delta t}{\Delta x} \left( Q^{k}_{N}-Q^{k}_{N-1} \right) +\sqrt{\frac{1}{gb}} \begin{Bmatrix} \frac{S_{o}^{1/2}}{n} \left[ \frac{A^{5/3}}{\left( b \left( 1-\frac{\sqrt{z^{2}+1}}{z} \right) + \frac{\sqrt{z^{2}+1}}{z} \sqrt{b^{2}+4zA} \right)^{2/3}} \right]_{N}^{k+1} - Q^{k}_{N} \\ + \frac{\Delta t}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N} \\ - \frac{\Delta t}{\Delta x} \left(\frac{Q^{2}}{2A} + \frac{gby_{c}^{2}}{2}+\frac{gzy_{c}^{3}}{3} \right)^{k}_{N-1} \\ - g\Delta t \left[ A\left(S_{o}-S_{f} \right) \right]_{N}^{k} \end{Bmatrix} = 0 \qquad[3.9] \]