Repositorio Institucional de Documentos

Modelo de redes de flujo con transición dinámica de régimen lámina libre/presión Trabajo Fin de Máster Máster Universitario en Mecánica Aplicada Programa Oficial de Postgrado en Ingeniería Mecánica y Materiales Curso 2011/12 Septiembre 2012 Autor: Javier Fernández Pato Directora: Dra. Pilar García-Navarro Universidad de Zaragoza Escuela de Ingeniería y Arquitectura Modelo de redes de flujo con transición dinámica de régimen lámina libre/presión Resumen La simulación numérica de flujos de agua en sistemas de drenaje urbano es uno de los ámbitos en donde se pone de manifiesto la necesidad de combinar flujos en lámina libre a presión atmosférica con situaciones en las que el conducto se encuentra presurizado, tanto en régimen estacionario como transitorio. Las ecuaciones que gobiernan los dos tipos de flujo son diferentes, por lo tanto, es necesario tener en cuenta el cambio lámina libre/presión a la hora de programar un modelo numérico completo que sea capaz de resolver transitorios independientemente de la región de trabajo. En este trabajo se desarrolla un modelo de simulación numérica capaz de resolver redes de tuberías cuyo régimen mayoritario de funcionamiento sea el de lámina libre, pero que se puedan ver presurizadas ante situaciones puntuales. Para ello, se propone adaptar la formulación matemática a través del método de la rendija de Preissmann, mediante en cuál se consigue una estimación razonable de la presión del agua en el conducto. El método numérico empleado se basa en el esquema de Roe de primer orden, enmarcado dentro de la familia de los métodos de volúmenes finitos. Se trata de un método adaptado a las situaciones transitorias bruscas, capaz de trabajar en régimen subcrítico, supercrítico y mixto. Para su validación, se han resuelto varios casos con solución analítica o datos empíricos correspondientes a experimentos de laboratorio. Mediante la aplicación a casos más complejos, como confluencias o redes de tuberías, se ha evaluado la sensibilidad del método a los cambios de régimen en este tipo de situaciones más realistas. Índice general 1. Introducción y objetivos....................................................................................................11 2. Modelo matemático...........................................................................................................13 2.1. Ecuaciones de conservación...................................................................................................13 2.1.1. Ecuación de conservación de la masa...........................................................................13 2.1.2. Ecuación de conservación del momento lineal...........................................................14 2.2. Flujo en lámina libre (shallow water)...................................................................................14 2.3. Sistemas transitorios presurizados (water hammer)..........................................................21 3. Modelo de la rendija de Preissmann...............................................................................27 4. Discretización mediante volúmenes finitos...................................................................29 4.1. Esquema de Roe explícito de primer orden.........................................................................29 4.2. Discretización de los términos fuente...................................................................................31 4.3. Condiciones de contorno........................................................................................................32 4.3.1. Condiciones de contorno físicas....................................................................................32 4.3.2. Condiciones de contorno y régimen de flujo...............................................................32 4.3.3. Confluencias.....................................................................................................................33 4.4. Condición de estabilidad........................................................................................................34 5. Validación del modelo.......................................................................................................37 5.1. Fondo con obstáculo y flujo estacionario.............................................................................37 5.2. Estado estacionario en un canal.............................................................................................39 5.3. Rotura de presa........................................................................................................................40 5.4. Test de Wiggert.........................................................................................................................43 5.5. Propagación de discontinuidades en flujo mixto................................................................45 5.6. Transitorio en flujo completamente presurizado................................................................46 6. Aplicación a redes..............................................................................................................49 6.1. Estado estacionario en una unión de conductos.................................................................49 6.2. Flujo transitorio en una unión de conductos.......................................................................51 6.3. Flujo transitorio en una red de tuberías...............................................................................54 7. Conclusiones y trabajo futuro..........................................................................................61 Apéndice A. Diagrama de flujo............................................................................................65 Apéndice B. Aplicación del método de Roe a las ecuaciones de lámina libre..............67 Índice de figuras Figura 1. Nomenclatura y sistema de coordenadas........................................................................15 Figura 2. Volumen de control con las fuerzas de presión, fricción y peso..................................15 Figura 3. Volumen de control.............................................................................................................21 Figura 4. Esquema de fuerzas............................................................................................................22 Figura 5. Esquema de la rendija de Preissmann..............................................................................27 Figura 6. Situación de lámina libre (izquierda) y presurización (derecha).................................27 Figura 7. Esquema de la propagación de la señal a través de las paredes de una celda...........30 Figura 8. Condiciones de contorno a la entrada (izq) y a la salida (dcha) para el caso subcrítico............................................................................................................................................33 Figura 9. Condiciones de contorno a la entrada (izq) y a la salida (dcha) para el caso supercrítico........................................................................................................................................33 Figura 10. Esquema de la confluencia...............................................................................................33 Figura 11. Ejemplo de un caso estable (izq) con CFL<1 y un caso inestable (dcha) con CFL>1. .............................................................................................................................................................35 Figura 12. Caso test #1.1 - Flujo transcrítico con onda de choque. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal...............................38 Figura 13. Caso test #1.2 - Flujo subcrítico. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal..................................................................38 Figura 14. Caso test #1.3 - Flujo transcrítico sin onda de choque. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal...............................38 Figura 15. Caso test #2.1. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal............................................................................................39 Figura 16. Caso test #2.2. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal............................................................................................40 Figura 17. Caso test #3.1......................................................................................................................41 Figura 18. Caso test #3.2. Nivel de agua (arriba) y caudal (abajo) en los tiempos t=0.3 s (rojo), t=1 s (verde) y t=2 s (azul). Estado inicial (linea discontinua)....................................................41 Figura 19. Caso test #3.3. Nivel de agua (arriba) y caudal (abajo) en los tiempos t=0.3 s (rojo), t=1 s (verde) y t=2 s (azul). Estado inicial (linea discontinua)....................................................42 Figura 20. Caso test #3.3 - Nivel de agua (arriba) y caudal (abajo) para el estado estacionario. .............................................................................................................................................................42 Figura 21. Montaje experimental de Wiggert..................................................................................43 Figura 22. Condiciones de contorno aguas arriba (izq) y aguas abajo (dcha).............................43 Figura 23. Resultados numéricos obtenidos para las cuatros sondas (CFL=0.9)........................44 Figura 24. Resultados numéricos obtenidos para la sonda 2 con CFL=0.6 (izq) y CFL=0.75 (dcha)..................................................................................................................................................44 Figura 25. Resultados numéricos obtenidos para la sonda 2 con anchuras de 3 cm (izq) y 1,5 cm (dcha)............................................................................................................................................44 Figura 26. Resultados numéricos obtenidos para la sonda 2 con Δx=0.5 m (arriba, izq), Δx=0.125 m (arriba, dcha), Δx=0.05 m (abajo, izq) y Δx=0.025 m (abajo, dcha).......................45 Figura 27. Propagación de la discontinuidad. Nivel de agua en los tiempos t=100 s (rojo), t=200 s (verde) y t=300 s (azul). Estado inicial (linea discontinua). Techo y fondo del conducto (gris)...................................................................................................................................46 Figura 28. Simulación del golpe de ariete. Nivel de agua en los tiempos t=3 s (rojo), t=6 s (verde) y t=9 s (azul). Estado inicial (linea discontinua). Techo del conducto (gris)...............47 Figura 29. Efectos en la solución al modificar la anchura de la rendija. Nivel de agua en los tiempos t=3 s (rojo), t=6 s (verde) y t=9 s (azul). Estado inicial (linea discontinua). Techo del conducto (gris)...................................................................................................................................47 Figura 30. Esquema de la unión.........................................................................................................49 Figura 31. Caso #6.1. Flujo subcrítico. Estado inicial (linea discontinua). Fondo del conducto (gris)....................................................................................................................................................50 Figura 32. Caso #6.2. Flujo subcrítico en la confluencia. Estado inicial (linea discontinua). Fondo del conducto (gris)................................................................................................................50 Figura 33. Caso #6.3. Flujo supercrítico en la confluencia. Estado inicial (linea discontinua). Fondo del conducto (gris)................................................................................................................51 Figura 34. Señal triangular para el caudal de entrada....................................................................52 Figura 35. Caso transitorio sin presurización (QMÁX=2.8 m3/s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=100 celdas..................................................................................................................................................52 Figura 36. Caso transitorio con presurización (QMÁX=3.2 m3/s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=100 celdas..................................................................................................................................................53 Figura 37. Caso transitorio sin presurización (QMÁX=2.8 m3/s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=200 celdas..................................................................................................................................................53 Figura 38. Caso transitorio con presurización (QMÁX=3.2 m3/s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=200 celdas..................................................................................................................................................54 Figura 39. Vista en planta y en perfil de la red de siete tuberías..................................................55 Figura 40. Estado estacionario para la red de siete tuberías. Estado inicial (linea discontinua). Fondo de los conductos (gris).........................................................................................................56 Figura 41. Calado en función del tiempo en el centro (rojo) y al final (azul) de cada tramo. Caudal máximo = 2.0 m3/s. Δx=10 m. N=10 celdas......................................................................57 Figura 42. Caudal en función del tiempo en el centro (rojo) y al final (azul) de cada tramo. Caudal máximo = 2.0 m3/s. Δx=10 m. N=10 celdas......................................................................58 Figura 43. Calado en función del tiempo en el centro (rojo) y al final (azul) de cada tramo. Caudal máximo = 3.0 m3/s. Δx=10 m. N=10 celdas......................................................................59 Figura 44. Caudal en función del tiempo en el centro (rojo) y al final (azul) de cada tramo. Caudal máximo = 3.0 m3/s. Δx=10 m. N=10 celdas......................................................................60 Capítulo 2. Modelo matemático ∫ t 1 t 2 [ (ρA u) x 2 (ρ A u) x 1 ] dt+∫ x 1 x 2 [ (ρ A) t 2 (ρ A) t 1 ] dx=0 (9) Por otro lado, la ecuación de conservación del momento lineal (5) aplicada a este volumen de control es: ∑ F= d dt ∫ V ρ v dV +∫ A ρ v (  v· n ) dA (10) La variación de momento lineal en el volumen de control entre t 1 y t 2 ha de ser igual a la suma de las fuerzas exteriores más el flujo neto de cantidad de movimiento que entra en dicho volumen durante ese intervalo de tiempo. Tomando la componente paralela a la dirección del flujo ( x ), el momento lineal por unidad de longitud será ρAu y su flujo a través de una sección transversal ρAu 2 . Entonces, el flujo neto de cantidad de movimiento entre t 1 y t 2 es: ∫ t 1 t 2 [ (ρ Au 2 ) x 2 (ρ Au 2 ) x 1 ] dt (11) Por otro lado, el incremento de momento lineal en el volumen de control será la suma de los incrementos infinitesimales en cada elemento diferencial de volumen: ∫ x 1 x 2 [ (ρ Au) t 2 (ρ A u) t 1 ] dx (12) Por último, es necesario calcular las fuerzas sobre el fluido, asumiendo que solamente actúan el peso, la fuerza de presión y la fricción con las paredes del canal. La componente x del peso se puede expresar en términos del ángulo θ . Para ello, recurrimos a la hipótesis inicial que asume que el ángulo de la pendiente del fondo es pequeño: S 0 =tg θ=∂ z ∂x≃senθ (13) Entonces, la componente del peso paralela a la dirección de la corriente es: ∫ x 1 x 2 ρg A senθdx≃∫ x 1 x 2 ρg A S 0 dx (14) 16 Capítulo 2. Modelo matemático Por lo tanto, la acción de dicha fuerza en el intervalo temporal dado será: ∫ t 1 t 2 ∫ x 1 x 2 ρg A S 0 dx dt (15) La fuerza de fricción entre el fluido y las paredes del cauce por unidad de longitud se puede expresar como: ρ g A S f (16) donde S f es la pendiente de la linea energética. Entonces, la acción de la fuerza de rozamiento será: ∫ t 1 t 2 ∫ x 1 x 2 ρg A S f dx dt (17) Para la fuerza de presión en la dirección del flujo sobre el volumen de control, distinguiremos entre las fuerzas sobre las paredes líquidas ( x 1 y x 2 ) y las fuerzas sobre las paredes sólidas. Estas últimas solamente estarán presentes en canales no prismáticos. La fuerza de presión ejercida sobre cualquier sección transversal sigue una distribución hidrostática: g I 1 = ∫ 0 h ρg(hη)dη (18) donde h es el calado (ver figura 1). La acción de la fuerza de presión neta entre t 1 y t 2 es: ∫ t 1 t 2 g [ ( I 1 ) x 2  ( I 1 ) x 1 ] dt (19) Por último, evaluaremos la componente de la reacción en las paredes debida a variaciones en la anchura del canal sobre un volumen diferencial: g I 2 dx= [ ∫ 0 h ρg(hη) ∂σ(x ,η) ∂xdη ] dx (20) Integrando a todo el volumen de control: 17 Capítulo 2. Modelo matemático ∫ x 1 x 2 g I 2 dx (21) La acción de esta componente de la fuerza es: ∫ t 1 t 2 ∫ x 1 x 2 g I 2 dx dt (22) Entonces, la ecuación de conservación para el momento lineal queda de la siguiente forma: ∫ x 1 x 2 [ (ρ Au) t 2 (ρ A u) t 1 ] dx+∫ t 1 t 2 [ (ρ A u 2 ) x 2 (ρ A u 2 ) x 1 ] dt+⋯ ⋯+∫ t 1 t 2 g [ ( I 1 ) x 2  ( I 1 ) x 1 ] dt∫ t 1 t 2 ∫ x 1 x 2 g I 2 dx dt+⋯ ⋯+∫ t 1 t 2 ∫ x 1 x 2 ρg A S f dx dt∫ t 1 t 2 ∫ x 1 x 2 ρg A S 0 dx dt=0 (23) Para obtener la forma diferencial de las ecuaciones (9) y (23), supondremos que las variables son funciones continuas y diferenciables. Entonces, mediante desarrollos en serie de Taylor: (ρ A) t 2 =(ρ A) t 1 +∂ρ A ∂ t ∆t+O(∆t 2 ) (ρ Au) x 2 =(ρ A u) x 1 +∂ρ A u ∂ t ∆x+O(∆ x 2 ) (24) y despreciando los términos de segundo orden, eliminando la densidad y suponiendo que ∆x→0 y ∆t→0 : lim t 1 →t 2 ∫ x 1 x 2 ( A t 2 A t 1 ) dx=∫ t 1 t 2 ∫ x 1 x 2 ∂ A ∂tdx dt (25) Procediendo de forma análoga en la ecuación de conservación de momento lineal obtenemos el sistema completo en forma diferencial conservativa. ∂ A ∂ t +∂ Q ∂ x =0 (26) ∂Q ∂t+∂ ∂x ( Q 2 A+g I 1 ) =g I 2 +g A ( S 0 S f ) (27) 18 Capítulo 2. Modelo matemático donde Q=Au es el caudal. El sistema puede formularse en forma vectorial: ∂ U ∂ t +∂ F ∂ x =R (28) donde los vectores U, F y R representan las variables conservadas, los flujos y los términos fuente, respectivamente: U=( A , Q ) T (29) F= ( Q , Q 2 A+g I 1 ) T (30) R= ( 0, g I2+g A ( S0Sf ) ) T (31) Además, se puede demostrar, mediante la regla de Leibniz, que: ∂ I 1 ∂ x =I 2 +A∂h ∂ x (32) Para el término de fricción se utilizará el modelo de Manning, descrito por: S f =Q∣Q∣n 2 A 2 R 4/3 (33) donde n es el número de Manning y R es el radio hidráulico, definido en términos del perímetro mojado P : R= A P (34) La matriz jacobiana del sistema es: J=∂F ∂U= ( 0 1 c 2 u 2 2u ) (35) donde c es la velocidad de propagación de las ondas en el modelo de aguas poco profundas: 19 Capítulo 2. Modelo matemático c= √ g∂I 1 ∂ A (36) Los autovalores y autovectores del sistema son los siguientes: λ 1,2 = u ± c , e 1,2 =( 1, u ± c ) T (37) Como la matriz jacobiana (35) es cuadrada de orden m=2, posee dos autovalores reales y es diagonalizable (sus autovectores son linealmente independientes), se dice que el sistema pertenece a la familia de ecuaciones hiperbólicas. Para el caso particular de un canal prismático rectangular: c= √ gA b = √ g h , I 1 =A 2 2 b , ∂ I 1 ∂ A =A b , I 2 = 0 (38) En este tipo de problemas, resulta común la caracterización del tipo de flujo a través del número de Froude (análogo al número de Mach en flujo compresible). Se trata de una magnitud adimensional definida por: Fr= u c (39) El número de Froude representa el balance entre la inercia del fluido y las fuerzas gravitatorias. En función de su valor, se pueden dar tres situaciones distintas: •Flujo subcrítico (Fr<1): la velocidad del fluido es menor que la velocidad de las perturbaciones, por lo que el flujo se encuentra controlado principalmente por la fuerza gravitatoria. •Flujo supercrítico (Fr>1): la velocidad del fluido es mayor que la velocidad de las perturbaciones. Se trata de un flujo rápido, donde el fluido no recibe ninguna información proveniente de la región de aguas abajo. Como se verá más adelante, esta transición puede derivar en la aparición de fenómenos como el salto hidráulico. •Flujo crítico (Fr=1): representa la frontera entre las dos situaciones anteriores. Es importante resaltar la importancia de este balance entre ambas velocidades, ya que los autovalores de la matriz jacobiana del sistema dependen directamente de ello. En concreto, el signo de dichos autovalores representa el sentido en el que se puede propagar la 20 Capítulo 2. Modelo matemático información en el flujo. 2.3. Sistemas transitorios presurizados (water hammer) En este apartado se presentan las ecuaciones diferenciales que gobiernan el comportamiento de un sistema presurizado en régimen transitorio (Chaudry et al. 1994). Al igual que las ecuaciones de lámina libre, es habitual escribirlas como una simplificación de las ecuaciones generales de conservación, bajo una serie de hipótesis: •Los flujos radiales de masa y cantidad de movimiento se suponen despreciables frente a los axiales, por lo tanto, se considerará que el flujo es unidireccional. •Se asumirá un flujo compresible y una tubería elástica. •Se tomará el promedio en la sección transversal de las variables conservadas, presión p (o altura piezométrica H ) y velocidad v. •Se aplicará un modelo de fricción para flujo estacionario (Darcy-Weissbach) para considerar las pérdidas energéticas con las paredes del conducto. Bajo estas condiciones, seguiremos un proceso análogo al de la sección anterior para escribir las ecuaciones de conservación de masa y momento lineal para un sistema presurizado. Consideremos un tramo de tubería elástica a través del cuál circula un fluido compresible, con un volumen de control asociado representado en la figura 3. Si suponemos que el fluido siempre está en contacto con las paredes de la tubería, se cumple que v=w en dicha región. De la ecuación (4): d dt ∫ V(t) ρdV +∫ A(t) ρ v R ·n dA=0⇒ d dt ∫ 0 ∆x ρA dx+∑ρ v R ·n A=0(40) 21 Figura 3. Volumen de control. Capítulo 2. Modelo matemático El último sumando solamente es distinto de cero en la entrada y en la salida, como se puede ver en la figura 3, por lo tanto: ∂ ∂t ( ρA∆x ) +∑ salida ρv R A∑ entrada ρv R A=0⇒ ⇒ ∂ ∂ t ( ρA∆x ) +ρ A v+ ∂ ∂ x ( ρA v ) ∆xρ A v=0 (41) Simplificando, obtenemos la ecuación de continuidad en forma diferencial: ∂ ∂ t ( ρA ) + ∂ ∂ x ( ρAv ) =0(42) Si reescribimos la ecuación anterior en términos del flujo másico: ˙ m= dm dt (43) ∂ ∂ t ( ρA ) = ∂ ∂ x ( ρA v ) =∂˙ m ∂ x (44) A la vista de la expresión anterior, queda claro que la compresibilidad del fluido (a través de la densidad) y la elasticidad del conducto (a través de la sección) pueden afectar al gradiente de flujo másico en la tubería. Para la deducción de la ecuación de conservación de la cantidad de movimiento, es necesario estudiar el balance total de fuerzas sobre el sistema, es decir, las fuerzas de presión 22 Figura 4. Esquema de fuerzas. Capítulo 2. Modelo matemático p y el esfuerzo τ 0 sobre las paredes de la superficie de control y el peso del fluido contenido en el volumen de control. En la figura 4 se muestra el esquema de fuerzas, así como los flujos de entrada y salida de cantidad de movimiento. Partiendo de la ecuación general (5) y analizando la componente en la dirección del flujo (x): ∑F x = d dt ∫ V(t) ρv x dV +∫ A(t) ρv x (  v R ·n ) dA⇒ ⇒p A+1 2 [ p+p+∂p ∂x∆x ] ∂A ∂x∆xp A∂ ∂x(p A)∆ x⋯ ⋯W senθτ 0 P∆x= d dt ∫ V(t) ρv x A dx+∑ρv x (  v R ·n ) A (45) donde P es el perímetro de la tubería y W el peso del fluido contenido en el volumen de control. Teniendo en cuenta que W=ρgA∆x y senθ=∆z/∆x=∂z/∂x : 1 2 ·∂ p ∂ x ·∂ A ∂ x (∆ x) 2 +p∂ A ∂ x ∆x∂ ∂ x (p A)∆ xρ g A ∆zτ 0 P∆x=⋯ ⋯= ∂ ∂ t (ρ A v)∆ x+ρ A v2+ ∂ ∂ x ( ρA v2 ) ∆xρ A v2 (46) Por último, dividiendo por ∆x y tomando el límite ∆x→0 : A ( ∂p ∂x+ρg∂z ∂x ) τ 0 P=∂ ∂t ( ρA v ) +∂ ∂x ( ρA v 2 ) (47) Es posible reescribir la ecuación (47) de la siguiente forma: ρ g A ( 1 ρg·∂p ∂x+∂z ∂x ) τ 0 P=v [ ∂(ρ A) ∂t+∂(ρ A v) ∂x ] +ρ A ( ∂v ∂t+v∂v ∂x ) (48) El tercer término de la expresión (48) es idénticamente cero, por la ecuación de continuidad (42). Introduciendo la altura piezométrica: H= p ρg+z(49) 23 Capítulo 2. Modelo matemático ρ g A ( ∂H ∂x ) τ 0 P=ρ A ( ∂v ∂t+v∂v ∂x ) (50) Si se particulariza la ecuación (50) para el caso de un conducto circular de diámetro D , obtenemos la siguiente expresión: P =π D ,A=π D 2 4 →1 g ( ∂v ∂t+v∂v ∂x ) =∂H ∂x 4 τ 0 ρg D (51) Es habitual encontrar las ecuaciones de conservación deducidas anteriormente en términos de las variaciones de presión. Para poder expresarlas de esta forma, definimos el módulo de compresibilidad del fluido: K= dp dρ/ρ (52) Además, en una tubería circular de espesor e en la que se cumpla que p/2<<E/(D/e) , el módulo elástico se puede definir de la siguiente forma: E= D e dp dA / A (53) Con las definiciones (52), (53) y un pequeño desarrollo matemático, es posible reescribir la ecuación de continuidad en función de las variaciones de presión: ∂(ρ A ) ∂ t =∂ρ ∂ t A+ρ ∂A ∂ t =∂p ∂ t ρ K A+ρ ∂p ∂ t ·D A E e (54) ∂(ρ A v ) ∂ t =∂ρ ∂ t Av+ρ ∂A ∂ t v+ρ A∂v ∂ x =ρ K ·∂p ∂ x A v+ρ D A E e ·∂p ∂ x v+ρ A∂v ∂ x (55) Sustituyendo (54) y (55) en la ecuación de continuidad (42): 24 Capítulo 2. Modelo matemático ∂p ∂ t ·ρ A K +ρ ∂p ∂ t D A E e +ρ K ∂p ∂ x A v+ρ D A E e ·∂p ∂ x v+ρ A∂v ∂ x =0⇒ ⇒∂p ∂tρ ( 1+D K e E ) +v∂p ∂xρ ( 1+D K e E ) +ρ K∂v ∂x=0⇒ ⇒ρ ∂p ∂t+ρv∂p ∂x+ρ K ( 1+D K e E ) ∂v ∂x=0 (56) Simplificando: ∂ p ∂t+v∂ p ∂x+ρc WH 2 ∂ v ∂x=0(57) donde c WH es la velocidad de propagación de las ondas de presión: c WH = √ K/ρ 1+D K e E (58) En el caso particular de una tubería inelástica: E→∞⇒CWH → √ K/ρ (59) La ecuación (59) corresponde a la velocidad máxima de las ondas de presión en la tubería. Siguiendo un procedimiento análogo para la ecuación de cantidad de movimiento, llegamos a la siguiente expresión: ∂v ∂t+v∂v ∂x+1 ρ·∂p ∂x+g senθ+ 4 τ 0 ρD=0 (60) También es frecuente expresar las ecuaciones de conservación en términos de la altura piezométrica, definida en (49): H= p ρg+z⇒p=ρ g(Hz)⇒dp=dρ· g(Hz)+ρ g(dH dz)(61) Introduciendo el módulo de compresibilidad del fluido (52): 25 Capítulo 4. Discretización mediante volúmenes finitos U i n+1 =U i n ∆ t ∆x [ ( ∑ k+ (  λ k  α k  β k )  e k ) i1/2 + ( ∑ k- (  λ k  α k  β k )  e k ) i+1/2 ] (86) Al igual que con los flujos, la influencia de los términos fuentes sobre la celda i se divide en la contribuciones a través de las paredes izquierda (i-1/2) y derecha (i+1/2). 4.3. Condiciones de contorno 4.3.1. Condiciones de contorno físicas. Para poder resolver completamente un problema de forma numérica es necesario establecer una serie de condiciones de contorno en los extremos del dominio de discretización. Para ello, si fuera necesario y en función del régimen de flujo, se impondrán condiciones de contorno físicas en las variables de área y/o caudal a la entrada y salida del sistema. Dichos valores en la frontera son aquellos que vienen impuestos por las condiciones reales del problema y que influirán en gran medida en la solución del mismo. En este trabajo se han empleado fundamentalmente tres tipos: nivel de agua constante, caudal uniforme e hidrograma de entrada en función del tiempo para representar, por ejemplo, una onda de avenida. Extendiendo el cálculo de las contribuciones disponibles (78) a cada celda hasta los propios contornos y usando el hecho de que en las paredes exteriores no existe intercambio de información (ver figuras 8 y 9), es posible calcular todos los puntos del mallado de forma numérica. Posteriormente, se deben imponer las condiciones de contorno físicas que correspondan. 4.3.2. Condiciones de contorno y régimen de flujo. El número de condiciones de contorno físicas que se deben imponer en los contornos depende del tipo de régimen (subcrítico o supercrítico), ya que éste se encuentra estrechamente relacionado con el sentido en el que se propaga la información. En concreto, se pueden dar las cuatro posibilidades que se detallan a continuación: •Flujo subcrítico a la entrada: se debe imponer solamente una condición de contorno física. La otra variable vendrá determinada por el esquema numérico (figura 8). •Flujo supercrítico a la entrada: se han de fijar las dos variables, ya que es imposible recibir información de la celda adyacente (figura 9). 32 Capítulo 4. Discretización mediante volúmenes finitos •Flujo subcrítico a la salida: al igual que en el caso de flujo subcrítico a la entrada, se impone una condición de contorno física y otra numérica (figura 8). •Flujo supercrítico a la salida: no es necesario imponer condiciones de contorno físicas. Toda la información proviene de la penúltima celda (figura 9). Figura 8. Condiciones de contorno a la entrada (izq) y a la salida (dcha) para el caso subcrítico. Figura 9. Condiciones de contorno a la entrada (izq) y a la salida (dcha) para el caso supercrítico. 4.3.3. Confluencias En los casos en los que se considere la confluencia de tres tuberías necesitamos establecer tres condiciones de contorno adicionales para poder resolver por completo el problema. En la figura 10 se muestra un esquema de las celdas que participan en la confluencia: Por un lado, se exigirá que la altura de agua a la entrada de las tuberías 2 y 3 sea igual a la salida de la tubería 1: 33 Figura 10. Esquema de la confluencia. Capítulo 4. Discretización mediante volúmenes finitos h 2 (1)=h 3 (1)=h 1 (I MAX ) En cuanto al caudal, debemos distinguir si el flujo es subcrítico o supercrítico: Caso subcrítico: Q 1 (I MAX )=Q 2 (1)+Q 3 (1) Caso supercrítico: Q 2 (1)=Q 3 (1)= 1 2Q 1 (I MAX ) Para un caso más general con N tuberías: h 1 = h 2 =⋯= h N , ∑ i=1 N Q i =0 (87) Si en la confluencia se considera la existencia de un pozo con una sección en planta A w la condición de contorno para el caudal se ve modificada de la siguiente manera: h 1 = h 2 =⋯= h N = H w , ∑ i=1 N Q i =A w dH w dt (88) donde H w representa la altura de agua en el pozo. 4.4. Condición de estabilidad En general, un método numérico se considera estable si las perturbaciones de la solución se mantienen acotadas. De lo contrario, el error en la solución crece de forma exponencial y la calidad de los resultados numéricos se ve seriamente comprometida. Una de las principales causas de las inestabilidades numéricas es el hecho de que la región de influencia numérica sea menor que la región de influencia física (ver figura 11). En los esquemas explícitos, como el empleado en este trabajo, la región de influencia numérica viene determinada por el tamaño de celda ∆x , ya que el valor de una variable depende de los valores de las celdas contiguas. Por otro lado, la región de influencia física viene dada por la distancia a la cual se ha podido propagar la información a velocidad c , es decir, ( |u|±c)∆t . Por lo tanto, una primera forma de escribir la condición de estabilidad es: ∆ x ⩾ ( ∣ u ∣± c ) ∆ t (89) De la expresión anterior, se puede deducir que el tamaño máximo para el paso temporal es: 34 Capítulo 4. Discretización mediante volúmenes finitos ∆t máx =∆ x máx ( ∣u∣+c ) (90) Definiendo el número de Courant-Friedrichs-Lewy (CFL) como: CFL=∆ t ∆t máx (91) se obtiene que la condición de estabilidad para un esquema explícito consiste en mantener este coeficiente por debajo de la unidad. 35 Figura 11. Ejemplo de un caso estable (izq) con CFL<1 y un caso inestable (dcha) con CFL>1. Capítulo 5. Validación del modelo 5. Validación del modelo En este capítulo se aplicará el modelo descrito anteriormente a casos de prueba que, o bien tienen solución exacta o han sido medidos en el laboratorio, con el objetivo de sopesar la validez del mismo. Para todos los casos de este capítulo se han empleado canales o conductos prismáticos de sección rectangular. 5.1. Fondo con obstáculo y flujo estacionario A continuación se reproducirán los tres casos test propuestos en Murillo et al. 2012, en los cuales el fondo del canal presenta una elevación dada por la siguiente función: z ( 8 ⩽ x ⩽ 12 )= 0.2  0.05 ( x  10 ) 2(92) La longitud y anchura del canal son 25 m y 1 m, respectivamente, y las condiciones iniciales son: h ( x ,0 )= 0.5  z ( x ) , u ( x ,0 )= 0 (93) En función de las condiciones de contorno que se impongan a la entrada y a la salida, se obtendrá flujo subcrítico o una transición sub-supercrítico, con o sin onda de choque. Los tres casos test se presentan en la siguiente tabla: Test Q aguas arriba h aguas abajo #1.1 (flujo transcrítico con onda de choque) 0.18 m 3 /s 0.33 m #1.2 (flujo transcrítico) 4.42 m 3 /s 2.0 m #1.3 (flujo transcrítico sin onda de choque) 1.53 m 3 /s 0.66 m (sub) Tabla 1. Fondo con obstáculo y flujo estacionario. En ninguno de los tres casos se ha considerado la fricción con el fondo. Los resultados numéricos obtenidos para el estado estacionario se observan en las figuras 12, 13 y 14: 37 Capítulo 5. Validación del modelo 38 Figura 12. Caso test #1.1 - Flujo transcrítico con onda de choque. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal. Figura 13. Caso test #1.2 - Flujo subcrítico. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal. Figura 14. Caso test #1.3 - Flujo transcrítico sin onda de choque. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal. 0 1 2 3 4 5 0 5 10 15 20 25 Q (m3/s) x (m) 0 0.5 1 1.5 2 2.5 0 5 10 15 20 25 h+z (m) x (m) 0 0.5 1 1.5 2 0 5 10 15 20 25 Q (m3/s) x (m) 0 0.2 0.4 0.6 0.8 1 1.2 0 5 10 15 20 25 h+z (m) x (m) 0 0.2 0.4 0.6 0.8 1 0 5 10 15 20 25 Q (m3/s) x (m) 0 0.1 0.2 0.3 0.4 0.5 0 5 10 15 20 25 h+z (m) x (m) Capítulo 5. Validación del modelo Se puede demostrar matemáticamente (y comprobar de forma numérica) que la transición sub-supercrítico tiene lugar en la parte más alta del obstáculo. 5.2. Estado estacionario en un canal Los siguientes casos se realizan en García-Navarro et al. 1993 y consisten en reproducir los estados estacionarios en un canal prismático bajo diferentes condiciones de entrada/salida, así como rozamiento. Como condiciones iniciales se tomará un caudal de 3 m 3 /s y un calado de 2 m. En la siguiente tabla se reflejan los parámetros empleados para cada uno de los casos: Test Fondo Nº de Manning CC aguas arriba CC aguas abajo #2.1 z(x) = 4.0 – 0.01 0.03 Q=3.0 m 3 /s h=2.0 m #2.2 z(x) = 4.0 – 0.01 0.009 Q=3.0 m 3 /s h=3.0 m Tabla 2. Estado estacionario en un canal 39 Figura 15. Caso test #2.1. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal. 0 1 2 3 4 5 0 50 100 150 200 250 300 350 400 Q (m3/s) x (m) 0 1 2 3 4 5 6 0 50 100 150 200 250 300 350 400 h+z (m) x (m) Capítulo 5. Validación del modelo En el primer caso (figura 15) se observa el perfil de un flujo completamente subcrítico, mientras que en el segundo caso (figura 16) los regímenes subcrítico y supercrítico se conectan a través de un salto hidráulico. 5.3. Rotura de presa El tercer caso test corresponde a la evolución temporal de una rotura de presa, es decir, partiremos de un desnivel inicial de agua en reposo situado en la mitad del canal y estudiaremos la evolución de las ondas de choque y rarefacción. Se trata de un problema clásico con solución exacta (Stoker 1957) para la validación de esquemas en régimen transitorio sin término de fricción. A continuación, se muestran los resultados obtenidos: Test Ratio Nº de Manning #3.1 1m : 0.5m 0 #3.2 10m : 1m 0 #3.3 10m : 1m 0.03 Tabla 3. Parámetros para las diferentes roturas de presa. 40 Figura 16. Caso test #2.2. Arriba: nivel de agua (azul), estado inicial (linea discontinua) y altura del fondo (gris). Abajo: caudal. 0 1 2 3 4 5 0 50 100 150 200 250 300 350 400 Q (m3/s) x (m) 0 1 2 3 4 5 6 0 50 100 150 200 250 300 350 400 h+z (m) x (m) Capítulo 5. Validación del modelo En la figura 17, correspondiente al primer caso sin rozamiento se muestra también la comparación con la solución exacta propuesta por Stoker (1957). En las figuras 18 y 19 se puede apreciar la evolución temporal de las ondas de choque y rarefacción: 41 Figura 18. Caso test #3.2. Nivel de agua (arriba) y caudal (abajo) en los tiempos t=0.3 s (rojo), t=1 s (verde) y t=2 s (azul). Estado inicial (linea discontinua). Figura 17. Caso test #3.1 0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1 0,4 0,5 0,6 0,7 0,8 0,9 1 t=0,05 s (numérica) t=0,05 s (analítica) t=0 s x (m) h (m) 0 10 20 30 40 50 60 0 20 40 60 80 100 Q (m3/s) x (m) 0 2 4 6 8 10 0 20 40 60 80 100 h+z (m) x (m) Capítulo 5. Validación del modelo Como se puede ver en la figura 29, la solución numérica es este caso está muy alejada de la realidad (tanto la presión como la velocidad de la discontinuidad), debido a que la velocidad de las ondas superficiales es mucho menor (87,8 m/s) que la velocidad de las ondas de presión en la tubería (1000 m/s). 48 Capítulo 6. Aplicación a redes 6. Aplicación a redes En este apartado aplicaremos el modelo a distintas conexiones de tuberías que, bajo ciertas condiciones de contorno aguas arriba, pueden verse sometidas a presurización en algunos puntos. Se analizarán tanto casos estacionarios como transitorios. 6.1. Estado estacionario en una unión de conductos Comenzaremos con un caso estacionario de unión en Y (figura 30), propuesto en García-Navarro et al. 1993, en el que una tubería principal se bifurca en dos ramas secundarias iguales. Las tres tuberías tienen geometría rectangular y comparten longitud (L=400 m), anchura (b=1 m) y coeficiente de Manning (n=0.009). Las pendientes y las condiciones iniciales y de contorno se especifican en la tabla 4: Caso S 0 (1) S 0 (2) S 0 (3) Condiciones iniciales CC entrada CC salida #6.1 0.001 0.001 0.001 Q 1 (i)=3.0 m 3 /s; Q 2 (i)=Q 3 (i)=1.5 m 3 /s; h 1 (i)=h 2 (i)=h 3 (i)=2.0 m Q=3.0 m 3 /s h=3.0 m #6.2 0.01 0.001 0.001 Q=3.0 m 3 /s h=3.0 m #6.3 0.01 0.01 0.01 Q=3.0 m 3 /s h=3.0 m Tabla 4. Estado estacionario en una unión de tuberías. En las figuras 31, 32 y 33 se presentan los resultados para el nivel de agua en los tres casos propuestos: 49 Figura 30. Esquema de la unión. Capítulo 6. Aplicación a redes 50 Figura 31. Caso #6.1. Flujo subcrítico. Estado inicial (linea discontinua). Fondo del conducto (gris). 0 0.5 1 1.5 2 2.5 3 3.5 0 50 100 150 200 250 300 350 400 h+z (m) (T3) x (m) 0 0.5 1 1.5 2 2.5 3 3.5 0 50 100 150 200 250 300 350 400 h+z (m) (T2) x (m) 0 0.5 1 1.5 2 2.5 3 0 50 100 150 200 250 300 350 400 h+z (m) (T1) x (m) Figura 32. Caso #6.2. Flujo subcrítico en la confluencia. Estado inicial (linea discontinua). Fondo del conducto (gris). 0 0.5 1 1.5 2 2.5 3 3.5 0 50 100 150 200 250 300 350 400 h+z (m) (T3) x (m) 0 0.5 1 1.5 2 2.5 3 3.5 0 50 100 150 200 250 300 350 400 h+z (m) (T2) x (m) 0 1 2 3 4 5 6 0 50 100 150 200 250 300 350 400 h+z (m) (T1) x (m) Capítulo 6. Aplicación a redes En estos dos últimos casos, aparece un salto hidráulico en el punto donde se produce el cambio de régimen super-subcrítico. 6.2. Flujo transitorio en una unión de conductos A continuación se desarrollará un caso propuesto en Wixcey 1990, en el cual se considera una bifurcación similar a la del apartado anterior. La longitud de las tuberías en este caso es de 5 km, la anchura es de 1 m, y el número de Manning es 0.01 para todos los tramos. La altura de los conductos también se considerará constante e igual a 1 m. Las pendientes de las tuberías principal y secundarias son 0.002 y 0.001, respectivamente. Primero se ha calculado un estado estacionario partiendo de las siguientes condiciones iniciales y de contorno: Q 1 (i)=0.1m 3 /s,Q 2 (i)=Q 3 (i)=0.05m 3 /s, h 1 ( i )= h 2 ( i )= h 3 ( i )= 0.2 m ,Q 1 (1)=0.1m 3 /s Tomando como condición inicial dicho estado estacionario, se impondrá como condición de contorno a la entrada una función de onda triangular para el caudal con las características mostradas en la figura 34: 51 Figura 33. Caso #6.3. Flujo supercrítico en la confluencia. Estado inicial (linea discontinua). Fondo del conducto (gris). 0 1 2 3 4 5 6 0 50 100 150 200 250 300 350 400 h+z (m) (T3) x (m) 0 1 2 3 4 5 6 0 50 100 150 200 250 300 350 400 h+z (m) (T2) x (m) 0 1 2 3 4 5 6 0 50 100 150 200 250 300 350 400 h+z (m) (T1) x (m) Capítulo 6. Aplicación a redes Como valor mínimo para el caudal tomaremos el mismo que en la simulación del estado estacionario, Q MIN =0.1 m 3 /s. En cuanto al valor de pico, se realizarán simulaciones con dos valores distintos. En el primer caso, tomaremos un valor de 2.8 m 3 /s, de forma que en ningún momento la altura de agua sea superior al techo de la tubería y, por lo tanto, el sistema no entre en condiciones de presurización. Posteriormente, se repetirá la simulación con un valor de 3.2 m 3 /s, donde entrará en juego el método de la rendija de Preissmann para estimar el cálculo de la presión. En las figuras 35 y 36 se presentan los resultados obtenidos con CFL=0.9 y una malla de 100 celdas para los dos casos descritos. En las figuras 37 y 38 se muestran los resultados de la simulación con una malla de 200 celdas. Vemos que los resultados son similares, aunque la forma de la onda se captura con mayor precisión, debido al refinamiento de la malla. 52 Figura 34. Señal triangular para el caudal de entrada. Figura 35. Caso transitorio sin presurización (Q MÁX =2.8 m 3 /s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=100 celdas. 0 0.2 0.4 0.6 0 1000 2000 3000 4000 5000 6000 h(m) (T3) t (s) 0 0.2 0.4 0.6 0 1000 2000 3000 4000 5000 6000 h(m) (T2) t (s) 0 0.2 0.4 0.6 0.8 1 0 1000 2000 3000 4000 5000 6000 h(m) (T1) t (s) Capítulo 6. Aplicación a redes 53 Figura 36. Caso transitorio con presurización (Q MÁX =3.2 m 3 /s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=100 celdas. 0 0.2 0.4 0.6 0.8 0 1000 2000 3000 4000 5000 6000 h(m) (T3) t (s) 0 0.2 0.4 0.6 0.8 0 1000 2000 3000 4000 5000 6000 h(m) (T2) t (s) 0 0.2 0.4 0.6 0.8 1 1.2 0 1000 2000 3000 4000 5000 6000 h(m) (T1) t (s) Figura 37. Caso transitorio sin presurización (Q MÁX =2.8 m 3 /s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=200 celdas. 0 0.2 0.4 0.6 0 1000 2000 3000 4000 5000 6000 h(m) (T3) t (s) 0 0.2 0.4 0.6 0 1000 2000 3000 4000 5000 6000 h(m) (T2) t (s) 0 0.2 0.4 0.6 0.8 1 0 1000 2000 3000 4000 5000 6000 h(m) (T1) t (s) Capítulo 6. Aplicación a redes 6.3. Flujo transitorio en una red de tuberías En este apartado se repetirá el experimento anterior en una red de siete tuberías, dispuestas como se indica en la figura 39. Todos los tramos tienen una longitud de 100 m y un coeficiente de Manning igual a 0,01. Las pendientes son uniformes para cada tramo y sus valores son los siguientes: S 01 = S 07 = 0.002 , S 02 = S 03 = S 05 = S 06 = 0.001 , S 04 = 0 Al igual que en el caso anterior, se realizará un primer cálculo hasta conseguir un estado estacionario, partiendo de las siguientes condiciones iniciales y de contorno: Q 1 ( i )= Q 7 ( i )= 0.1 m 3 / s,Q 2 ( i )= Q 3 ( i )= Q 5 ( i )= Q 6 ( i )= 0.05m 3 / s,Q 4 ( i )= 0m 3 / s h 1 ( i )= h 2 ( i )= h 3 ( i )= h 4 ( i )= h 5 ( i )= h 6 ( i )= h 7 ( i )= 0.2 m , Q 1 ( 1 )= 0.1m 3 / s Las condiciones para las confluencias J 1 y J 2 son las mismas que las empleadas en el apartado 6.1, mientras que en las uniones W 1 y W 2 se ha supuesto la existencia de un pozo con una sección en planta A w =5 m 2 por lo que la condición de contorno para el caudal se ve modificada como se especificó en la sección 4.3.3. 54 Figura 38. Caso transitorio con presurización (Q MÁX =3.2 m 3 /s). Calado en función del tiempo para los puntos x=500 m (rojo), x=1000 m (verde) y x=5000 m (azul). CFL=0,9. N=200 celdas. 0 0.2 0.4 0.6 0.8 0 1000 2000 3000 4000 5000 6000 h(m) (T3) t (s) 0 0.2 0.4 0.6 0.8 0 1000 2000 3000 4000 5000 6000 h(m) (T2) t (s) 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0 1000 2000 3000 4000 5000 6000 h(m) (T1) t (s) Capítulo 6. Aplicación a redes Partiendo del estado estacionario obtenido, se ha modificado la condición de contorno aguas arriba, imponiendo la función triangular descrita en la figura 34. El caudal máximo se ha establecido en 2.0 m 3 /s para una primera simulación en la que no se llega a presurizar la red en ningún punto. Si se aumenta el caudal de pico hasta 3.0 m 3 /s se conseguirá una presurización parcial en algunos puntos del sistema. En las figuras 40, 41, 42, 43 y 44 se presentan los resultados obtenidos para el calado y el caudal de agua para los estados estacionario, transitorio sin presurizar y transitorio presurizado, con CFL=0,9. En las figuras correspondientes al caudal de agua en la red, se puede ver que el caudal en el centro de la tubería es idénticamente cero, mientras que en los extremos opuestos es igual y de signo contrario, por lo que se pone de manifiesto la simetría del problema. Debido a dicha simetría, las tuberías 3 y 6 no son representadas, ya que los resultados son idénticos a los de los conductos 2 y 5. 55 Figura 39. Vista en planta y en perfil de la red de siete tuberías. Figura 40. Estado estacionario para la red de siete tuberías. Estado inicial (linea discontinua). Fondo de los conductos (gris). 0 0.1 0.2 0.3 0.4 0 20 40 60 80 100 h7+z7 (m) x (m) 0 0.1 0.2 0.3 0 20 40 60 80 100 h5+z5 (m) x (m) 0 0.1 0.2 0 20 40 60 80 100 h4+z4 (m) x (m) 0 0.1 0.2 0.3 0 20 40 60 80 100 h2+z2 (m) x (m) 0 0.1 0.2 0.3 0.4 0 20 40 60 80 100 h1+z1 (m) x (m) Figura 41. Calado en función del tiempo en el centro (rojo) y al final (azul) de cada tramo. Caudal máximo = 2.0 m 3 /s. Δx=10 m. N=10 celdas. 0 0.2 0.4 0.6 0.8 0 500 1000 1500 2000 h7(m) t (s) 0 0.2 0.4 0.6 0.8 0 500 1000 1500 2000 h5(m) t (s) 0 0.2 0.4 0.6 0.8 0 500 1000 1500 2000 h4(m) t (s) 0 0.2 0.4 0.6 0.8 0 500 1000 1500 2000 h2(m) t (s) 0 0.4 0.8 1.2 0 500 1000 1500 2000 h1(m) t (s) Bibliografía 4368, 2010 [15] Murillo, J., García-Navarro, P., Augmented versions of the HLL and HLLC Riemann Solvers including source terms in one and two dimensions for shallow flow applications, Journal of Computational Physics, 2012 [16] Rocha Felices, Arturo, Hidráulica de tuberías y canales, Universidad Nacional de Ingeniería, 2007 [17] Stoker, J.J., Water waves, Wiley Interscience, 1957 [18] Toro, E.F., Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer-Verlag Berlin Heidelberg, 1999 [19] Trajkovic, B., Ivetic, M., Calomino, F., D'Ippolito, A., Investigation of transition from free surface to pressurized flow in a circular pipe, Water Science and Technology, 39(9), 105-112, 1999 [20] Villanueva Lacabrera, Ignacio, Simulación numérica de flujos estacionarios y transitorios en ríos y canales, 1999 [21] Wiggert, D., Transient flow in free-surface, pressurized systems, Journal of the Hydraulics Division, Proceedings of the American Society of Civil Engineers 98 (1)(1972) 11-26, 1972 [22] Wixcey, J.R., An investigation of algorithms for open channel flow calculations, Numerical Analysis Internal Report, 21, Departament of Mathematics, University of Reading, 1990 64