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^{2/3}S_{o}^{1/2}\qquad[1.2]\]
      \[\frac{Qn}{S_{o}^{1/2}}=y^{5/3}b\qquad[1.3]\]
      \[Z=y^{5/3}b\qquad[1.4]\]

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

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

      el valor de \(y\) se usará como condicón inicial para encontrar el valor del tirante en la Ec. [1.1]


    • Método de Newton-Raphson
      •

      Para poder usar este método necesitamos reformular la Ec. [1.1] y hallar su derivada de la siguiente forma

      \[f(y)=0=\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}\) (ver Ec. [1.5]) y luego se reemplaza 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. Discretización de las ecuaciones de movimiento


    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{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^{2}}{A}+gAy_{c}\right)=gA\left(S_{o}-S_{f}\right) \qquad[2.2]\]

    Las ecuaciones de movimiento se puede reescribir en forma compacta o vectorial como sigue

    \[ \frac{\partial }{\partial t} \begin{bmatrix} A \\ Q\end{bmatrix} + \frac{\partial }{\partial x} \begin{bmatrix} Q \\ \frac{Q^{2}}{A}+gAy_{c} \end{bmatrix} = \begin{bmatrix} 0 \\ gA\left(S_{o}-S_{f}\right)\end{bmatrix}\]
    \[ \frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}=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 Lax es un método explícito de diferencias finitas, de primer orden en el tiempo y segundo orden en el espacio, que se emplea principalmente como herramienta introductoria para el análisis numérico de ecuaciones de conservación hiperbólicas. Su formulación incorpora un promedio espacial de las variables, lo que introduce disipación numérica artificial, otorgándole estabilidad y robustez frente a discontinuidades. Debido a su carácter difusivo y baja precisión, no suele utilizarse en aplicaciones prácticas avanzadas, pero resulta fundamental para comprender el comportamiento, las limitaciones y la evolución hacia esquemas más sofisticados.Los operadores diferenciales para el esquema de LAX son

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

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

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

    donde la pendiente de fricción \(S_{f}\) se evalua en la posición \(i+1\) e \(i-1\) con la siguiente expresión

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

    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.7]\]

    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.8] \]

    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.9] \]

    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.10] \]

  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}} \left[ \frac{\partial Q}{\partial t}+\frac{\partial }{\partial x}\left(\frac{Q^{2}}{2A}+gAy_{c}\right)-gA\left(S_{o}-S_{f}\right)\right]=0 \qquad[3.2] \]

      reemplazamos las derivadas parciales por los operadores diferenciales en la Ec. [3.2]