Repositorio Institucional de Documentos

Algoritmos de regulaci´ on en flujos transitorios bidimensionales Mario Morales Hern´andez M´aster en Mec´anica Aplicada Programa Oficial de Posgrado en Ingenier´ıa Mec´anica y de Materiales Septiembre 2010 Dirigido por: Dra. Pilar Garc´ıa Navarro Curso 2009-2010 Centro Polit´ ecnico Superior Universidad de Zaragoza 2 Algoritmos de regulaci´on en flujos transitorios bidimensionales Resumen El uso de m´etodos num´ericos para predecir la forma de la l´amina de agua o las variaciones de caudal en caso de flujos tanto estacionarios como no estacionarios es, hoy en d´ıa, una pr´actica necesaria dentro de la tecnolog´ıa hidr´aulica moderna, ya que ofrece la posibilidad de evaluar, de una forma no muy costosa, la respuesta de los sistemas hidr´aulicos frente a una gran variedad de situaciones pr´acticas. En particular, la simulaci´on num´erica de flujos transitorios permite acceder, entre otras cosas, a la predicci´on de tiempos de tr´ansito de las ondas, alcanzados por los valores m´aximos en distintos puntos del dominio de c´alculo. Esta informaci´on es esencial para el dise˜no de estrategias de prevenci´on y puede ser combinada con algoritmos de regulaci´on y optimizaci´on. Los problemas que acarrean las grandes precipitaciones y las consiguientes inundaciones es un tema de inter´es creciente en los ´ultimos a˜nos. Grandes p´erdidas econ´omicas y sociales son las principales consecuencias que traen consigo estas inundaciones. Para poder actuar de una manera efectiva contra esto, una opci´on es la utilizaci´on de ´areas de inundaci´on controlada para almacenar la mayor parte de la crecida y devolver el agua al r´ıo una vez haya finalizado. Esto es, por otra parte, un reto, debido a que la aplicaci´on de estos conceptos a cuencas urbanizadas que fueron desarrolladas hace mucho tiempo con zonas bajas ocupadas por infraestructura industrial, comercial y residencial se antoja muy complicada. El trabajo consiste en familiarizarse con modelos de simulaci´on transitoria para flujos bidimensionales existentes en el grupo de investigaci´on, recopilar informaci´on de elementos y algoritmos de regulaci´on, establecer un caso test de aplicaci´on y programar la combinaci´on de ambos elementos para ensayar su potencial utilidad. 3 4 ´ Indice general 1. Introducci´on 11 2. Modelo 2D de flujo de l´amina libre con promedio en la vertical 15 2.1. Ecuaciones generales ............................. 15 2.2. Ecuaciones promediadas en la vertical .................... 18 2.2.1. T´ermino de fricci´on y modelos de turbulencia ........... 19 2.2.2. Versi´on com´un de las ecuaciones de aguas poco profundas . . . . 20 2.3. Esquema num´erico .............................. 22 2.4. Modelos de compuertas ............................ 24 3. Caso test 27 3.1. Cauce ..................................... 27 3.2. Elevaci´on topogr´afica y malla ........................ 28 3.2.1. Elevaci´on topogr´afica ......................... 28 3.2.2. Malla ................................. 29 3.3. Condiciones iniciales y condiciones de contorno ............... 29 3.3.1. Condiciones iniciales ......................... 30 3.3.2. Condiciones de contorno ....................... 30 4. Influencia de determinados factores en la reducci´on del pico de caudal 33 4.1. An´alisis dimensional ............................. 33 4.2. Consideraciones previas ............................ 35 4.3. Caso 1 : Π∆tG0= 0.9 ΠV= 0.00817 ..................... 35 4.4. Caso 2: Πtac = 0.4 ΠV=0.00817 ....................... 39 4.5. Caso 3: Π∆tG0= 0.9 Πtac =0.4 ........................ 42 5 5. Algoritmos de regulaci´on 45 5.1. Controlador b´asico On/Off ......................... 45 5.2. Controlador PID: Proporcional, Integral y Diferencial ........... 46 5.2.1. El algoritmo .............................. 47 5.2.2. Representaci´on discreta del PID ................... 48 5.3. Implementaci´on en el c´odigo de simulaci´on y funcionamiento ....... 49 5.4. Resultados ................................... 51 6. Conclusiones y trabajo futuro 57 6.1. Conclusiones .................................. 57 6.2. Trabajo futuro ................................ 59 Bibliograf´ıa 60 A. Esquema num´erico 65 B. Teorema Πde Buckingham 69 B.1. Demostraci´on ................................. 69 B.2. Aplicaci´on del teorema Pi .......................... 71 C. Condiciones de contorno 75 C.1. Condiciones de contorno en la ecuaci´on lineal escalar 2D ......... 75 C.1.1. Celda ficticia ............................. 77 C.2. Condiciones de contorno para el sistema de ecuaciones 2D ........ 78 C.2.1. Aplicaci´on al modelo bidimensional de las ecuaciones de aguas poco profundas ............................... 78 C.3. Condiciones de contorno internas ...................... 80 D. Ajuste de controladores PID 83 D.1. M´etodo de la curva de reacci´on ....................... 84 D.2. Aplicaci´on a nuestro caso test ........................ 86 6 ´ Indice de figuras 2.1. Perfil del cauce. ................................ 16 2.2. Ejemplo para la aplicaci´on de la ecuaci´on de Bernoulli .......... 24 2.3. Niveles de agua para el caudal en la ecuaci´on (2.38) ............ 25 2.4. Niveles de agua para el caudal en la ecuaci´on (2.39) ........... 25 2.5. Niveles de agua para el caudal en la ecuaci´on (2.40) ............ 26 2.6. Niveles de agua para el caudal en la ecuaci´on (2.41) ........... 26 3.1. Secci´on transversal del r´ıo .......................... 28 3.2. Caso test .................................... 29 3.3. Malla ...................................... 30 3.4. Hidrograma gaussiano y SCS ......................... 32 4.1. Influencia del instante de apertura ...................... 36 4.2. Reducci´on del tiempo que est´a abierta la compuerta ............ 37 4.3. Influencia real del instante de apertura ................... 38 4.4. Influencia del tiempo en el que est´a abierta la compuerta ......... 40 4.5. Influencia tiempo en el que est´a abierta la compuerta ........... 41 4.6. Influencia real del tiempo que est´a abierta la compuerta ......... 42 4.7. Influencia del volumen en el amortiguamiento del pico de caudal ..... 43 4.8. Hidrograma entrante y laminaci´on (a) .................... 44 4.9. Hidrograma entrante y laminaci´on (b) ................... 44 5.1. Compuerta tipo Narmix ........................... 46 5.2. Ejemplo de una instalaci´on de una compuerta Narmix .......... 46 5.3. Celdas involucradas en la compuerta 3 ................... 50 5.4. Esquema de funcionamiento ......................... 51 7 5.5. Hidrogramas elegidos para las simulaciones ................. 52 A.1. Representaci´on constante de las varaibles en cada celda. ......... 66 A.2. Selecci´on de la informaci´on necesaria en el m´etodo descentrado. ..... 68 D.1. Control de una planta ............................ 83 D.2. Control de una planta ............................ 84 D.3. Curva de reacci´on t´ıpica ........................... 85 D.4. Curva de reacci´on en nuestro caso test ................... 86 8 ´ Indice de cuadros 3.1. Hidrograma SCS unitario ........................... 31 4.1. Casos a simular ................................ 34 D.1. Determinaci´on de las constantes propuestas por Ziegler-Nichols y Cohen- Coon ..................................... 85 D.2. Determinaci´on de las constantes propuestas por Ziegler-Nichols y Cohen- Coon ..................................... 86 9 tipos: cinem´aticas y din´amicas. Primero vamos a plantear las condiciones cinem´aticas en el fondo y en la superficie libre y despu´es haremos lo mismo con las condiciones de frontera de tipo din´amico. Las condiciones cinem´aticas est´an relacionadas con la velocidad y nos dicen que las part´ıculas de agua en su movimiento no pueden cruzar ninguna frontera. Para el fondo significa que la componente de la velocidad normal a la superficie s´olida debe ser cero (fondo s´olido, impermeable y fijo). vˆ nb=u∂zb ∂x +v∂zb ∂y −w= 0 (2.3) con ˆ nb= (∂zb/∂x, ∂zb/∂y, −1) el vector normal a la superficie l´ıquida en contacto con la s´olida hacia afuera en z=zb(x, y), donde zbes la cota del fondo medida desde un nivel de referencia horizontal (Fig. 2.1). En la superficie libre las cosas son un poco m´as complicadas ya que ´esta se puede mover. En este caso lo que debe ser nula es la velocidad normal relativa a esa superficie y representa que el fluido no se puede salir del propio fluido, no puede atravesar tampoco la superficie libre ∂H ∂t +u∂H ∂x +v∂H ∂y −w= 0 (2.4) en z=H(x, y, t), donde Hes el nivel de la superficie libre medido desde el nivel de referencia (Fig. 2.1). Zb h H Figura 2.1: Perfil del cauce. Queremos hacer notar aqu´ı que el cauce puede tener una pendiente elevada tanto en direcci´on xcomo en direcci´on yy por ello hay que elegir bien los ejes de referencia sobre los cuales se van a tomar medidas. En este caso, lo que se hace es suponer una referencia 16 horizontal arbitraria para simplificar notablemente la formulaci´on y se considera que el ´angulo que define la pendiente es casi despreciable. Las condiciones de contorno din´amicas nos dan informaci´on sobre las fuerzas que act´uan en los contornos. Si el flujo es viscoso y el fondo fijo, la fuerza que act´ua es la de la viscosidad y por lo tanto las part´ıculas que se encuentran en contacto con el fondo est´an pegadas a ´el, por lo cual puede imponerse la condici´on de no deslizamiento que conduce a: u=v= 0 (2.5) en z=zb(x, y). En la superficie libre se supone continuidad de esfuerzos; es decir, los esfuerzos en el fluido justo por debajo de la superficie libre son los mismos que los del aire justo encima. Estos esfuerzos pueden ser de dos tipos: normales y cortantes ´o tangentes a la superficie. En el caso de los esfuerzos normales a la superficie donde interviene el t´ermino de presi´on, despreciando los efectos de la tensi´on superficial, P=Pa(2.6) donde Paes la presi´on atmosf´erica. El nivel absoluto de presiones no es importante y se puede tomar como cero. Las diferencias s´olo podr´ıan ser importantes en el caso de que se quisiera estudiar el efecto de las variaciones de presi´on atmosf´erica en el movimiento del agua. Los esfuerzos tangenciales que act´uan en la superficie libre son debidos a esfuerzos viscosos y la condici´on de contorno nos dice que deben ser iguales al esfuerzo tangencial aplicado al otro lado de la frontera de la superficie libre y que puede ser provocado por el viento. De este modo, el modelo contempla que en la superficie del mar, por ejemplo, puede actuar un esfuerzo cortante debido al viento. Este esfuerzo cortante externo τs= (τsx, τsy) tangente a la superficie del agua es τsx = (T ·nH)x=−τxx ∂H ∂x −τxy ∂H ∂y +τxz (2.7) en z=Hy de forma similar para la direcci´on y. El vector de esfuerzos producido por el viento se supone conocido y se trata como una fuerza externa. La magnitud y la direcci´on de la fuerza del viento en la superficie del mar vienen determinadas por el flujo en la atm´osfera. Normalmente se supone que el m´odulo de la velocidad del viento es conocida Wy se acepta la f´ormula semi-emp´ırica dada por Gill [24], τs=ρcWW2(2.8) 17 y la direcci´on se supone que es la direcci´on de la velocidad que lleva el viento. El coeficiente cWno es constante, depende de la velocidad del viento; por ejemplo, es del orden de 0,001 si la velocidad del viento se mide a unos 10 mde altura. Resolver este sistema de ecuaciones es muy costoso computacionalmente ya que hace falta resolver la ecuaci´on tridimensional de Poisson [48] para obtener la distribuci´on de presiones en cada paso temporal. La condici´on de contorno de la superficie libre implica una no linealidad extra en el problema y para muchas aplicaciones se prefiere resolver el sistema de ecuaciones de aguas poco profundas que se obtiene siguiendo una serie de aproximaciones. Despu´es de realizar un estudio de las escalas caracter´ısticas del problema [9] llegamos a las ecuaciones de aguas poco profundas tridimensionales que reescribimos a continuaci´on: ∂u ∂x +∂v ∂y +∂w ∂z = 0 (2.9) ∂u ∂t +∂u2 ∂x +∂uv ∂y +∂uw ∂z =−g∂h ∂x +∂τxx ∂x +∂τxy ∂y +∂τxz ∂z (2.10) ∂v ∂t +∂uv ∂x +∂v2 ∂y +∂vw ∂z =−g∂h ∂y +∂τyx ∂x +∂τyy ∂y +∂τyz ∂z (2.11) ∂p ∂z +ρg = 0 (2.12) Se hace notar que con las ecuaciones de movimiento (2.10), (2.11) podemos conocer los valores de uyv,wla obtendr´ıamos a partir de la ecuaci´on de conservaci´on de la masa (2.9) y, por ´ultimo, la variable htendr´ıa que ser determinada por las condiciones de contorno de la superficie libre (2.4), siendo a´un as´ı un procedimiento todav´ıa complicado. Por esto se recurre al promedio en la vertical cuya hip´otesis fundamental es que Las ondas que se producen en la superficie var´ıan suavemente, lo cual es equivalente a decir que la distribuci´on de presiones en la vertical es hidrost´atica o que la aceleraci´on en la vertical es peque˜na. 2.2. Ecuaciones promediadas en la vertical Para pasar a la forma bidimensional de las ecuaciones, dejando la profundidad como variable dependiente, hay que dar un paso m´as. Con el objetivo de eliminar de las ecuaciones la informaci´on del movimiento en la direcci´on vertical zse promedian las ecuaciones en 18 esta direcci´on haciendo uso de las definiciones de los promedios de las variables. ¯u=1 hZH zb udz (2.13) ¯v=1 hZH zb vdz (2.14) El proceso de promediar en la vertical las ecuaciones convierte el problema tridimensional en uno bidimensional de grosor variable hdonde los contornos ya no est´an en la superficie libre y el fondo sino en el per´ımetro. El promedio en la vertical de las ecuaciones del flujo de superficie libre bajo las hip´otesis del modelo de aguas poco profundas (ver [9]) conduce a una versi´on muy com´un del sistema de ecuaciones en 2D que repetimos aqu´ı: ∂h ∂t +∂(hu) ∂x +∂(hv) ∂y = 0 ∂(hu) ∂t +∂(hu2) ∂x +∂(huv) ∂y =−gh∂H ∂x +cfu√u2+v2+hνT∇2u ∂(hv) ∂t +∂(huv) ∂x +∂(hv2) ∂y =−gh∂H ∂y +cfv√u2+v2+hνT∇2v 2.2.1. T´ermino de fricci´on y modelos de turbulencia El coeficiente cfque aparece en el t´ermino de fricci´on se expresa habitualmente en t´erminos del coeficiente de rugosidad de Manning no de Ch´ezy [14], cfu√u2+v2=n2u√u2+v2 h4 3 (2.15) cfv√u2+v2=n2v√u2+v2 h4 3 (2.16) El coeficiente de rugosidad nen la pr´actica se determina a partir de medidas experimentales o se estima a partir de valores que ya han sido almacenados en tablas [14]. La ecuaci´on de Manning aqu´ı descrita es de naturaleza emp´ırica y por tanto es el resultado de un proceso de ajuste a una curva de datos experimentales. La primera dificultad que surge a la hora de usar este coeficiente de rugosidad es la precisi´on con la que ha sido estimado. El coeficiente ndepende en principio del n´umero de Reynolds del flujo, de la rugosidad de los contornos y de la forma geom´etrica de la cuenca. La rugosidad de la 19 superficie del contorno representa un valor cr´ıtico a la hora de estimar n, con valores peque˜nos si el material es fino y valores altos en el caso contrario. El valor de ntambi´en debe de dar cuenta de la vegetaci´on retardando el flujo y proporcionando valores altos de n, dependiendo tambi´en de la altura de agua. El modelo de fricci´on dado por (2.15) y (2.16) se basa en la teor´ıa de capa l´ımite estacionaria sobre pared rugosa. Con el objeto de calcular las variables hidrodin´amicas, es necesario fijar el valor del coeficiente de Manning, como ya hemos dicho, y el valor de la viscosidad cinem´atica de remolino. La viscosidad turbulenta νTdepende de las caracter´ısticas del flujo y puede variar de un punto a otro del dominio. Por tanto, es necesario plantear un modelo de turbulencia que nos permita evaluar el valor de νTen cada punto del dominio. En [9] se presentan algunos de los modelos de turbulencia existentes. 2.2.2. Versi´on com´un de las ecuaciones de aguas poco profundas El t´ermino que proviene de promediar en la vertical el gradiente de presi´on ha dado lugar a los t´erminos g∂H/∂x,g∂H/∂y, que a su vez se pueden descomponer, teniendo en cuenta que H=h+zben g∂H ∂x =g∂h ∂x +g∂zb ∂x (2.17) g∂H ∂y =g∂h ∂y +g∂zb ∂y (2.18) Los t´erminos ∂h/∂x,∂h/∂y se agrupan junto a las otras derivadas del mismo tipo (t´erminos convectivos). Las variaciones del fondo se expresan en forma de pendiente S0x=−∂zb ∂x (2.19) S0y=−∂zb ∂y (2.20) Los t´erminos de fricci´on del agua con el fondo del cauce se representan por Sf, pendiente de la l´ınea de energ´ıa en cada direcci´on Sfx =cfu√u2+v2 gh (2.21) Sfy =cfv√u2+v2 gh (2.22) 20 dando lugar al siguiente sistema de ecuaciones que es la forma m´as conocida de representaci´on del modelo de aguas poco profundas ∂h ∂t +∂hu ∂x +∂hv ∂y = 0 (2.23) ∂hu ∂t +∂hu2 ∂x +gh∂h ∂x +∂huv ∂y =gh(S0x−Sfx) (2.24) ∂hv ∂t +∂huv ∂x +∂hv2 ∂y +gh∂h ∂y =gh(S0y−Sfy) (2.25) Este sistema de ecuaciones en su forma conservativa [1,17] es decir, escritas las ecuaciones de la forma m´as cercana posible a un sistema de leyes de conservaci´on de masa y cantidad de movimiento, es ∂U ∂t +∇E=S⇒∂U ∂t +∇·(F,G) = S(2.26) con U=   h hu hv   ,F=   hu hu2+gh2 2 huv   ,G=   hv huv hv2+gh2 2   , S=   0 gh (S0x−Sfx) gh (S0y−Sfy)   (2.27) Urepresenta el vector de variables conservadas (hprofundidad del agua (Fig. 2.1), hu yhv caudales unitarios a lo largo de las direcciones coordenadas x,yrespectivamente), FyGson los flujos de las variables conservadas a trav´es de los lados de un volumen de control, y contienen el flujo convectivo y los gradientes de presi´on hidrost´atica. La parte derecha de la igualdad en el sistema de ecuaciones, S, contiene las fuentes y sumideros de la cantidad de movimiento a lo largo de las dos direcciones coordenadas, provenientes de las variaciones del fondo del cauce y de las p´erdidas por fricci´on que deben estar relacionadas con el campo de velocidades. De esta manera, (2.26) representa un sistema hiperb´olico de ecuaciones diferenciales en derivadas parciales acopladas y no lineales. Si escribimos el sistema de ecuaciones en formulaci´on no conservativa ∂U ∂t + (A,B)·∇U=S(2.28) 21 las matrices Jacobianas de los vectores de flujo son A=∂F ∂U=   0 1 0 c2−u22u0 −uv v u   ,B=∂G ∂U=   0 0 1 −uv v u c2−v20 2v   (2.29) y la matriz Jacobiana del flujo normal a una direcci´on dada por ˆ nse puede escribir como A=Aˆnx+Bˆny=   0 ˆnxˆny −u(u·ˆ n) + c2ˆnxu·ˆ n+uˆnxuˆny −v(u·ˆ n) + c2ˆnyvˆnxu·ˆ n+uˆny   (2.30) Los valores propios del Jacobiano Jnson a1=u·ˆ n+c a2=u·ˆ n a3=u·ˆ n−c(2.31) y sus vectores propios e1=   1 u+cˆnx v+cˆny   ,e2=   0 −cˆny cˆnx   ,e3=   1 u−cˆnx v−cˆny   (2.32) 2.3. Esquema num´erico El dominio donde se mueve el flujo, se subdivide, en un conjunto de celdas para su resoluci´on num´erica. En el modelo presentado hay libertad a la hora de elegir el tipo de celdas: hex´agonos, cuadril´ateros, tri´angulos, etc... y adem´as pueden formar parte de una malla estructurada o de una malla no estructurada. La elecci´on de la malla es un factor importante en la simulaci´on num´erica. Respecto a la t´ecnica de resoluci´on de las ecuaciones, se ha usado un m´etodo de vol´umenes finitos porque combina lo mejor de los m´etodos de elementos finitos y su flexibilidad geom´etrica, con lo mejor de los m´etodos en diferencias finitas, su flexibilidad en la definici´on del flujo discreto (valores discretos de las variables dependientes y sus flujos asociados). El primer paso es escribir (2.28) en la forma ∂U ∂t +−→ ∇E=S(2.33) 22 donde E= (F,G)T. Aplicando el teorema de Gauss sobre la celda de c´alculo Ωifija en el tiempo, (2.33) se escribe como: ∂ ∂t ZΩi UdΩi+I∂Ωi Endl =I∂Ωi Tndl (2.34) donde Tn es un vector que expresa el t´ermino fuente a trav´es de la superficie ∂Ω, ldenota la variable de integraci´on de la superficie alrededor del volumen Ωiynes el vector exterior normal unitario. Una vez formulado el problema en vol´umenes finitos se ha de elaborar una estrategia adecuada para calcular el flujo num´erico a trav´es de la superficie. La forma hiperb´olica del sistema de ecuaciones hace que este problema sea resuelto adecuadamente utilizando un esquema num´erico perteneciente a la familia de los m´etodos de Godunov [29,46]. Este tipo de m´etodo calcula el flujo num´erico que actualiza el valor de cada celda de c´alculo ipromediando el valor de las diferentes soluciones aproximadas que aparecen al definir un problema de Riemann en la superficie entre el volumen y cada uno de los vol´umenes vecinos j. Dentro las posibles opciones que permiten generar una soluci´on aproximada, en este trabajo se utiliza la aproximaci´on propuesta por Roe [39], que a diferencia de otras considera todas las velocidades de propagaci´on de informaci´on contenidas en el Jacobiano de la matriz. Cuando aparecen t´erminos fuente formularlos a trav´es de una matriz en la pared, como en (A.1), permite desarrollar soluciones aproximadas m´as complejas adecuadamente [35]. De esta manera, el flujo normal En y su jacobiano cobran protagonismo. Se eval´ua la matriz jacobiana del flujo normal y se diagonaliza, permitiendo que el esquema num´erico se base en vectores y valores propios. De esta forma, el esquema num´erico utilizado para resolver la ecuaciones de l´amina libre se detalla en el Anexo A. En este trabajo se incluye la modelizaci´on de compuertas como agente regulador del flujo de agua. El flujo a trav´es una compuerta no puede ser definido como flujo de superficie libre y la hip´otesis de presi´on hidrost´atica ya no es v´alida. El sistema de ecuaciones de conservaci´on de masa y momento no es adecuado, se requiere la participaci´on de leyes de conservaci´on de energ´ıa. Este cambio en el sistema de ecuaciones se evita modelando las compuertas como una discontinuidad entre las superficies de las celdas, donde el flujo num´erico de Godunov entre dos celdas no es calculado. Con este fin, se definen condiciones de contorno internas donde hay que imponer un n´umero de variables adecuadas. La discretizaci´on de este tipo de condici´on de contorno que modela el flujo a presi´on en una compuerta y garantiza una correcta conservaci´on de la masa se detalla en el Anexo C. 23 2.4. Modelos de compuertas Para la realizaci´on de este trabajo se implementan algoritmos de regulaci´on de compuertas. Las compuertas se definen a trav´es de una condici´on de contorno interna en la que se impone un flujo m´asico que las atraviesa. Este flujo m´asico en la compuerta es modelado partiendo de un principio b´asico: el caudal unitario que atraviesa la compuerta viene gobernado por la diferencia de niveles superficiales [26], d=h+z, existentes a ambos lados de la compuerta, donde hes el calado y zla elevaci´on superficial. Para comprobar esto con un ejemplo sencillo supongamos que nos encontramos con un escenario como el que se ilustra en la Figura 2.2 en el cual el nivel de referencia es igual aguas arriba y aguas abajo de la compuerta. Figura 2.2: Ejemplo para la aplicaci´on de la ecuaci´on de Bernoulli Aplicamos la ecuaci´on de Bernoulli entre el punto 1 y en el punto 2 y llegamos a p ρg +z+v2 2g1 =p ρg +z+v2 2g2 (2.35) Ponemos p=pat +ρgh y simplificando se llega a que h1+v2 1 2g=h2+v2 2 2g(2.36) Imponemos por continuiddad h1v1=h2v2con h1≫h2, luego v1≪v2. Por lo tanto simplificando y despejando v2de la expresi´on (2.36) se tiene v2=p2g(h1−h2) (2.37) Para determinar el caudal unitario s´olo tendremos que multiplicar por la apertura de la compuerta. 24 En lo que sigue nos referiremos como d2al nivel aguas arriba de la compuerta y d1al nivel aguas abajo, as´ı como G0a la apertura de la compuerta. Esto permite incluir la presencia de discontinuidades en la elevaci´on del terreno exactamente donde est´a situada la compuerta. Obviamente el caso en el que G0= 0, la compuerta se convierte autom´aticamente en un pared s´olida. A partir de esta aclaraci´on b´asica, y sin p´erdida de generalidad, asumimos d2> d1(el caso d1≥d2es an´alogo). Con esta suposici´on, se generan cuatro posibles escenarios que dependen de los posibles niveles superficiales (d2yd1) y elevaci´on del terreno (z2yz1) aguas arriba y aguas abajo de la compuerta. Caso 1: z1>z2,d2−z1>G0,d1>z1+G0 Este caso se ilustra en la Figura 2.3 y el caudal que atraviesa la compuerta viene dado por la expresi´on q=G0K1(d2−d1)1 2(2.38) donde K1es una constante emp´ırica [26]. Caso 2: z1>z2,d2−z1>G0,d1≤z1+G0 Se ilustra en la Figura 2.4 y el caudal unitario se expresa mediante q=G0K2(d2−z1)1 2(2.39) donde K2es una constante emp´ırica [26]. Figura 2.3: Niveles de agua para el caudal en la ecuaci´on (2.38) Figura 2.4: Niveles de agua para el caudal en la ecuaci´on (2.39) 25 100 200 300 400 500 600 700 800 900 0 10000 20000 30000 40000 50000 60000 70000 80000 90000 Discharge (m3/s) Time (s) Hydrographs Scs hydrograph Gaussian hydrograph Figura 3.4: Hidrograma gaussiano y SCS Curva de aforo Para la condici´on de contorno a la salida se ha utilizado una ley de flujo normal: S0=1 1000 Sf=Q|Q|n2 A(h)2Rh(h)4/3(3.5) donde A(h) = ´area mojada, Rh(h) = Radio hidr´aulico que a su vez se calcula mediante Rh=A Pm donde Pmes el per´ımetro mojado. Imponemos S0=Sf, y realizando una tabla con valores discretos de (h,Q) obtenemos nuestra condici´on de contorno a la salida del cauce que introduciremos dentro de uno de los archivos de configuraci´on del simulador. Compuertas Las compuertas se definen a trav´es de una condici´on de contorno interna en la que se impone un flujo m´asico que las atraviesa. Este flujo m´asico en la compuerta es modelado utilizando K1= 3,33 y K2= 2,248 en (2.38-2.41). 32 Cap´ıtulo 4 Influencia de determinados factores en la reducci´on del pico de caudal 4.1. An´alisis dimensional Antes de abordar cualquier situaci´on es necesario saber qu´e par´ametros y en qu´e medida influyen en nuestro problema. Nuestro objetivo es intentar reducir el pico de caudal de un hidrograma mediante la utilizaci´on de ´areas de inundaci´on controlada. Una manera de abordar este primer contacto con el problema es la utilizaci´on del an´alisis dimensional. Nos referimos al an´alisis dimensional como aquellos procedimientos que, basados en el an´alisis de las variables y par´ametros que gobiernan un fen´omeno, y m´as espec´ıficamente en las magnitudes f´ısicas que dichas variables involucran, permiten encontrar relaciones entre par´ametros adimensionales. El problema f´ısico queda entonces descrito, con el mismo grado de fidelidad, por este nuevo conjunto reducido de par´ametros adimensionales. Enfatizamos la palabra reducido, dado que ´esta es una de las ventajas del an´alisis dimensional. Al ser menor el n´umero de variables o par´ametros, es posible organizar y expresar m´as eficientemente los resultados de la experimentaci´on. Una herramienta muy valiosa en el an´alisis dimensional es el teorema Πde Buckingham. Gracias a este teorema, es posible reducir el n´umero de par´ametros o variables de los cuales depende un fen´omeno f´ısico, mediante la generaci´on de grupos adimensionales que involucran dichas variables. Resulta particularmente valioso cuando no se conoce la ecuaci´on que gobierna un fen´omeno y se busca encontrar dicha relaci´on a trav´es de la experimentaci´on de laboratorio. En el anexo Cse enuncia y demuestra el mencionado teorema Π de Buckingham y se 33 aplica a nuestro caso de estudio. Supondremos que el descenso en el caudal pico depende de cinco variables, a saber, ∆Qp=f(∆tG0, G0, V, Qpe, ta) donde Qps= pico de caudal a la salida en m3/s Qpe= pico de caudal a la entrada m3/s tac = tiempo de apertura de la compuerta en s tcc = tiempo de cierre de la compuerta en s G0= apertura de la compuerta en m V = volumen del hidrograma ∆Qp=|Qps−Qpe|∆tG0=tcc −tac y habiendo adimensionalizado, probamos que ∆Qp Qpe =F(∆tG0Qpe V,tacQpe V,G0 vol1 3 ) Por lo tanto, s´olo nos queda estudiar la influencia de cada uno de estos par´ametros Π∆tG0=∆tG0Qpe V Πtac =tacQpe V ΠV=G0 vol1 3 en el descenso del pico de caudal. Para ello vamos a estudiar c´omo influyen estos par´ametros en el amortiguamiento del pico de una manera muy sencilla: mantendremos dos de ellos constantes y haremos variar el otro para ver la relaci´on existente entre el amortiguamiento y este ´ultimo par´ametro (ver Tabla 4.1). Π∆tG0Πtac ΠV Caso 1 cte var cte Caso 2 var cte cte Caso 3 cte cte var Cuadro 4.1: Casos a simular 34 4.2. Consideraciones previas Este an´alisis se va a realizar teniendo en cuenta la influencia de un solo dep´osito o ´area de inundaci´on. En concreto, ser´a siempre el tercer dep´osito (el que est´a aguas abajo). La compuerta permanecer´a, o bien cerrada, o bien abierta 2 m. Las simulaciones realizadas son ´unicamente de 40000 segundos. Los picos de caudal elegidos para realizar este estudio van desde 500 m3/s a 1000 m3/s aumentando de 50 en 50 m3/s considerando que partimos de un estado estacionario con caudal = 100 m3/s . Los hidrogramas utilizados para este tipo de simulaciones son hidrogramas de tipo SCS. El volumen del hidrograma unitario SCS es k = 1.357. Cuando damos valores a los pico caudal y a los tiempos de pico, el volumen del nuevo hidrograma se transforma de la manera siguiente: V=kTpeQpe(4.1) 4.3. Caso 1 : Π∆tG0= 0.9 ΠV= 0.00817 En este primer caso hemos intentado ver la relaci´on que hay entre el amortiguamiento en el pico de caudal y el instante de tiempo en el que se abre la compuerta, permaneciendo constante el tiempo durante el que est´a abierta esta misma compuerta. Imponer Π∆tG0y ΠVconstantes implicaba tomar por un lado un volumen constante del hidrograma y al variar el pico de caudal a la entrada, ir variando el tiempo que est´a abierta la compuerta para conseguir un valor constante en Π∆tG0. El volumen que se ha tomado es 14653980 m3. En este primer caso nos encontramos en la situaci´on de que Π∆tG0es constante, es decir, para cada pico de caudal tenemos un valor constante de ∆tG0; en otras palabras, la compuerta permanece abierta durante un per´ıodo fijo de tiempo. Los tiempos de apertura que se han elegido para las simulaciones viene de aplicar la siguiente f´ormula: Para cada pico de caudal, y para cada 1 ≤i≤10 hacemos 35 (ta)i=(Tpe+ 5000)i 10 Esta elecci´on de los tiempos de apertura se debe a intentar ver la influencia de elegir un tiempo de apertura anterior y posterior al tiempo de pico. Con todo esto, y haciendo un peque˜no recuento, para cada pico de caudal, hemos simulado 10 casos. Como tenemos 11 picos distintos de caudal hemos simulado 110 casos. En la Figura 4.1 se muestran los resultados. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ∆Q/Qp tap Qp/vol Qp=500 Qp=550 Qp=600 Qp=650 Qp=700 Qp=750 Qp=800 Qp=850 Qp=900 Qp=950 Qp=1000 Figura 4.1: Influencia del instante de apertura Conlusiones 1. En el eje de ordenadas se representa el descenso en el pico de caudal dividido por el caudal pico a la entrada. Este dato nos da una idea en tanto por uno de la cantidad de agua que se almacena en el dep´osito en comparaci´on con la que discurre por el cauce. Es decir, lo que buscamos nosotros es que ese valor sea lo m´as cercano a 1 posible. Cuanto m´as se aproxime este valor a 0 querr´a decir que la apertura que estamos considerando es menos eficiente. 36 2. En el eje de abscisas se representa el par´ametro adimensional Πtac 3. La primera conclusi´on es algo que ya pod´ıamos sospechar: si mantenemos constante el tiempo en el que est´a abierta la compuerta y abrimos la compuerta o bien muy pronto o bien muy tarde, el dep´osito no “reduce el hidrograma entrante” de una manera ´optima. ´ Esto lo vemos en los dos primeros tiempos de apertura y en los dos ´ultimos (se abre demasiado pronto y demasiado tarde respectivamente). 4. Tambi´en podemos observar que la funci´on meseta que se intuye en cada uno de los caudales pico es muy amplia, es decir, hay muchas aperturas en las cuales se produce un efecto ´optimo en la reducci´on del hidrograma. Esto posiblemente es debido a que el par´ametro en el que est´a incluido el tiempo que est´a abierta la apertura ( Π∆tG0=∆tG0Qpe V) es muy grande. Para probar este hecho se han realizado otras 440 simulaciones esta vez con un pico de caudal de 850 m3/s pero con Π∆tG0= 0.3, 0.5, 0.7 y 0.9 (es decir, hemos reducido el tiempo que est´a abierta la compuerta). Los resultados se muestran en la Figura 4.2. Se aprecia que conforme el valor Π∆tG0es m´as peque˜no se reduce el abanico de aperturas ´optimas. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.5 1 1.5 2 ∆Q/Qp ta Qp/vol Π2 = 0.3 Π2 = 0.5 Π2 = 0.7 Π2 = 0.9 Figura 4.2: Reducci´on del tiempo que est´a abierta la compuerta 37 5. Por ´ultimo la curvatura que se experimenta en los valores a medida que vamos aumentando el tiempo de apertura se debe a que (tac)i=(Tpe+ 5000)i 10 Luego por (4.1) (tac)iQp V=i 10k Tpe+ 5000 Tpe Qp Qp−100 Como QpyTpeson inversamente proporcionales se tiene (m´as o menos) que el cociente (Tpe+ 5000)Qp Tpe(Qp−100) est´a pr´oximo a 1 , por lo que no influye mucho el par´ametro i (ya veremos que en el caso de Πtac y ΠVconstantes es diferente). 0 0.05 0.1 0.15 0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 ∆Q / Ql tap Ql/vol Qp = 500 Qp = 550 Qp = 600 Qp = 650 Qp = 700 Qp = 750 Qp = 800 Qp = 850 Qp = 900 Qp = 950 Qp = 1000 Figura 4.3: Influencia real del instante de apertura Como se observa en la Figura 4.1, no hemos conseguido “adimensionalizar bien” puesto que las gr´aficas son muy similares pero siguen dependiendo del pico de caudal (no colapsan). En un intento por hacer que esto ´ultimo ocurra, se va a representar la Figura 4.1 de nuevo, pero, en vez deen el eje de ordenadas Qp−Qs Qp, se ha representado Ql−Qs Qldonde 38 Qlrepresenta el caudal a la salida cuando la compuerta est´a cerrada. Es decir, cuando representamos este par´ametro Ql−Qs Qlestamos midiendo verdaderamente la influencia del dep´osito puesto que dejamos a un lado lo que el propio cauce lamina cuando no hay ning´un dep´osito. En la Figura 4.3 se representa la influencia del instante de apertura en el amortiguamiento del pico (s´olo midiendo la influencia del dep´osito, dejando aparte lo que lamina el cauce). Se observa que, a partir de cierto pico de caudal (750 m3/s) todos los puntos colapsan, con lo que se genera una curva ´unica que relaciona el instante de apertura con el amortiguamiento a la salida. Si hubi´eramos representado m´as puntos, como hemos hecho en la Figura 4.2, ver´ıamos que la pendiente de subida es m´as suave de lo que aqu´ı se muestra. Tanto esta pendiente como la zona en la que es ´optimo el amortiguamiento dependen del par´ametro ∆tG0. 4.4. Caso 2: Πtac = 0.4 ΠV=0.00817 An´alogamente al Caso 1, se ha elegido un valor constante en el volumen del hidrograma de 14653980 m3. Como se ha variado el pico de caudal de igual manera que en el caso anterior y ahora se pretende que sea Πtac el que quede fijo, se propone para este Caso 2 variar el tiempo de apertura con el pico de caudal para conseguir un valor constante ( en este caso se ha elegido un valor de Πtac = 0.4 ). Lo que vamos a representar en la Figura 4.4 es la relaci´on que hay entre el amortiguamiento del pico de caudal y el tiempo en el que permanece abierta la compuerta. Estos tiempos se han variado de la manera siguiente: (∆tG0)i=(39500 −tac)i 10 para 1 ≤i≤10. Al igual que en el caso anterior hemos simulado por lo tanto 10 ∆tG0 diferentes para cada pico de caudal. Como tenemos 11 picos distintos de caudal, pues 110 simulaciones cuyos resultados se muestran en la Figura 4.4. Conlusiones 1. Igual que ocurr´ıa en el caso 1, en el eje de ordenadas se representa el amortiguamiento en el pico de caudal definido como ∆Qp Qpe; en el eje de abscisas es Π∆tG0lo que se dibuja. De nuevo, parecen intuirse las funciones meseta que se mencionaban anteriormente aunque han sido cortadas debido a que el tiempo de simulaci´on era solamente de 40000 segundos. Evidentemente si el instante de apertura es fijo y 39 el tiempo que se deja la compuerta abierta es muy peque˜no (valores de m´as a la izquierda en la gr´afica) el amortiguamiento es muy bajo. Sin embargo cuando nos movemos entre valores de 0.8 y 1 (siempre considerando que Πtac = 0.4 constante), en todos los casos simulados con caudales pico distintos, el aprovechamiento del dep´osito ya es ´optimo. 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 ∆Q/Qp ∆tap Qp/vol Qp=500 Qp=550 Qp=600 Qp=650 Qp=700 Qp=750 Qp=800 Qp=850 Qp=900 Qp=950 Qp=1000 Figura 4.4: Influencia del tiempo en el que est´a abierta la compuerta 2. Para comprobar el aspecto que tienen estas funciones meseta hemos simulado 100 casos con un pico de caudal de 850 m3/s y un tiempo de apertura de la compuerta constante (Πtac = 0.7). Se pueden ver los resultados en la Figura 4.5. La manera de interpretar este gr´afico es sencilla: como el tiempo en el que abrimos la compuerta es constante y lo que vamos variando es el tiempo que dejamos abierta la compuerta, cuanto m´as tiempo dejemos abierta la compuerta m´as agua entrar´a al dep´osito, con lo que el amortiguamiento del pico de caudal ser´a mayor. Esto es, cuando estamos m´as a la izquierda en el eje de las x los tiempos que dejamos abierta la compuerta son peque˜nos y el amortiguamiento es menor. Conforme nos vamos desplazando a la derecha lo que estamos haciendo es dejar la compuerta abierta durante m´as tiempo, con lo cual entra m´as agua al dep´osito. 40 3. Adem´as se observa que, conforme m´as grande es el par´ametro Πtac (esto es, el instante en el que se abre es mayor) m´as pronunciada es la pendiente de la funci´on meseta hasta alcanzar el amortiguamiento ´optimo. 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.5 1 1.5 2 ∆Q/Qp ∆tap Qp/vol Πtac = 0.7 Πtac =0.4 Figura 4.5: Influencia tiempo en el que est´a abierta la compuerta 4. En cuanto al desplazamiento hacia la derecha que se observa en los datos cuando aumentamos el caudal pico, es debido a que ahora el par´ametro que var´ıa lo hace de la siguiente manera: (∆tG0)i=(39500 −tac)i 10 Sabiendo que tacQp V= 0,4 se llega a que, (∆tG0)iQp V=i 10 (39500 −tac)0,4 tac = 0,04i(39500 tac −1) Ahora 39500 tac −1 es un n´umero mayor que uno que al multiplicarlo por ise har´a m´as grande, de ah´ı el desplazamiento que sufre la funci´on conforme m´as elevado es el pico de caudal. 41 a los cambios actuales en la se˜nal de error. Las partes integral, proporcional y derivada pueden interpretarse como acciones de control basadas en el pasado, el presente y el futuro [36]. Para poder utilizar este m´etodo es necesario, aparte de discretizar la ecuaci´on (5.1), ajustar los par´ametros K,TiyTdpara cada sistema mediante m´etodos tunning de controladores PID. En el anexo Dse explican los dos m´etodos que se han utilizado para el ajuste emp´ırico de nuestro sistema, y se proporcionan valores para las constantes mencionadas anteriormente. 5.2.2. Representaci´on discreta del PID Para discretizar el controlador es necesario aproximar la integral y la derivada – ecuaci´on (5.1) - por formas manejables para la computaci´on en un ordenador. La derivada se aproxima mediante diferencias finitas y se llega a que: de dt ≈e(tk)−e(tk−1) tk−tk−1 =e(tk)−e(tk−1) Ts (5.2) donde tkdenota los instantes de muestreo, es decir, los momentos en que el ordenador lee la variable continua y Tses el valor de dicho per´ıodo, llamado per´ıodo de muestreo. Es cada per´ıodo de muestreo cuando el controlador ejecuta la acci´on. Es decir, la nueva posici´on para la compuerta. Para la integral: Zt 0 e(t)≈Ts· t X i=0 e(i) (5.3) siendo Tsconstante. Por lo tanto el algoritmo discretizado del PID queda: u(tk) = K·e(tk) + K·Ts Td tk X i=0 e(i) + KTd·(e(tk)−e(tk−1)) Ts (5.4) Para poder obtener el algoritmo que implementaremos en el c´odigo se escribe la ecuaci´on (5.4) en terminos de k−1: u(tk−1) = K·e(tk−1) + K·Ts Td tk−1 X i=0 e(i) + K·Td·(e(tk−1)−e(tk−2)) Ts (5.5) Por lo tanto, el incremento de la acci´on generada por el controlador se obtiene restando las ecuaciones (5.4) y (5.5): 48 ∆u(tk) = u(tk)−u(tk−1) = K[e(tk)−e(tk−1)] + K·Ts Ti e(tk) +K·Td Ts [e(tk)−2e(tk−1) + e(tk−2)] (5.6) Agrupando t´erminos se tiene: ∆u(tk) = u(tk)−u(tk−1) = K·1 + Ts Ti +Td Tse(tk)−K·1 + 2Td Tse(tk−1) +K·Td Ts e(tk−2) (5.7) Con el fin de intentar obtener m´as formas de estabilizar el algoritmo en caso de problemas de inestabilidad se a˜naden al algoritmo PID los factores de peso α1,α2yα3, para obtener finalmente la representaci´on discreta de la ecuaci´on del PID: u(tk) = u(tk−1) + α1·K·1 + Ts Ti +Td Ts(href (tk)−h(tk)) −α2·K·1 + 2Td Ts(href (tk−1)−h(tk−1)) +α3·K·Td Ts (href (tk−2)−h(tk−2)) (5.8) donde href es el valor objetivo para la variable regulada (el nivel de agua) α1,α2,α3son los factores para ponderar cada paso temporal (´ındice n, n-1 y n-2) hes el valor actual de la variable controlada 5.3. Implementaci´on en el c´odigo de simulaci´on y funcionamiento Una vez ajustado el controlador y obtenidas las constantes (ver anexo D), solo falta introducir todas estas modificaciones en el c´odigo de simulaci´on de flujos superficiales que tenemos en el grupo de investigaci´on. Esto se ha hecho a trav´es de dos nuevos archivos de configuraci´on para el simulador. Por un lado,un archivo que permite identificar cu´antas compuertas existen en nuestro sistema y d´onde est´an colocadas. Para esto ´ultimo se necesitan proporcionar la informaci´on del numero de parejas de celdas que componen la compuerta, y identificador de cada celda involucrada. Podemos ver un ejemplo en la Figura 5.3. 49 Figura 5.3: Celdas involucradas en la compuerta 3 Este caso, que hace referencia a la compuerta 3, las parejas de celdas que intervienen son: 1224-11211 y 798-1225. Por otro lado, un nuevo archivo que involucra todos los aspectos relatados anteriormente referidos a los algoritmos de regulaci´on. En primer lugar se tiene que discernir entre que tipo de regulaci´on se quiere (On/Off o PID), hay que establecer el tiempo de muestreo, es decir, el tiempo para el cual se da una nueva apertura de la compuerta, as´ı como las aperturas m´axima y m´ınima de la compuerta. Adem´as tambi´en se especificaran los valores elegidos para los par´ametros αique provienen de la discretizaci´on del algoritmo, y las coordenadas xeydel punto elegido para ser la referencia y nivel de referencia (o setpoint) escogido. Igualmente se generan autom´aticamente dos archivos de salida, uno que nos proporciona las aperturas que se registran de cada una de las compuertas y otro en el que se detallan los niveles de referencia preestablecidos y los registrados en cada instante de tiempo. Para clarificar como se comporta el simulador una vez implementados los correspondientes algoritmos de regulaci´on se incluye el siguiente diagrama de flujo, en el que se esquematiza el funcionamiento del bloque de c´alculo. 50 Definición de aproximación a las derivadas en las interceldas. Definición de valores promediados en interceldas. Cálculo de flujos de información desde intercelda hacia celdas vecinas. Cálculo de (h,u,v) n+1 en función de (h,u,v) ny de información de las interceldas vecinas, y de descarga de las compuertas. Actualización del resto de variables. t = t +! t t = 0 K = 1 Obtención de apertura ( K +1) a partir de ecuación del controlador t = t + ! t K = K+1 ENTRADA DE DATOS ¿ t = KT muestreo ? SÍ NO ¿ t > T cálculo ? NO SÍ SALIDA DE DATOS Figura 5.4: Esquema de funcionamiento 5.4. Resultados En este apartado se presentan los resultados que se han obtenido a ra´ız de la implementaci´on en el c´odigo de simulaci´on de estos dos tipos de algoritmos de control diferentes. Aunque se han hecho numerosas simulaciones con diferentes tipos de hidrogramas, se incluyen tres ejemplos, cada uno con sus correspondientes m´etodos de control, para poder observar las diferencias entre unos y otros. 51 En la Figura 5.5 se pueden ver los tres hidrogramas mencionados. La elecci´on de estos tres hidrogaramas se debe a las numerosas diferencias que presentan entre s´ı, como por ejemplo los picos de caudal registrado, la forma de cada uno... De esta manera se ampl´ıa el abanico de posibles actuaciones frente a hidrogramas diferentes. 0 500 1000 1500 2000 2500 3000 0 20000 40000 60000 80000 100000 120000 Caudal (m3/s) Tiempo(s) Hidrograma 1 Hidrograma 2 Hidrograma 3 Figura 5.5: Hidrogramas elegidos para las simulaciones Con cada uno de los tres hidrogramas que aqu´ı se presentan, se han simulado 4 casos: sin ´areas de inundaci´on controlada, es decir, las compuertas cerradas, sin regulaci´on, o sea, las compuertas completamente abiertas, con una regulaci´on On/Off y con un controlador PID. Como resultados se presentan las figuras siguientes: para cada hidrograma y cada mecanismo de regulaci´on se incluyen dos figuras, una en la que se representan los caudales a la entrada y a la salida, as´ı como las aperturas que se obtienen y otra en la que se puede ver los niveles previamente establecidos (setpoint) y los registrados durante el per´ıodo de avenida. A continuaci´on se presentan dichos resultados. 52 0 200 400 600 800 1000 1200 1400 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 0 200 400 600 800 1000 1200 1400 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 Apertura=0.0 On/Off 0 200 400 600 800 1000 1200 1400 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 0 200 400 600 800 1000 1200 1400 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 PID Apertura=8.0 0 0.5 1 1.5 2 2.5 3 3.5 4 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 0 0.5 1 1.5 2 2.5 3 3.5 4 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 Apertura=0.0 On/Off 0 0.5 1 1.5 2 2.5 3 3.5 4 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 0 0.5 1 1.5 2 2.5 3 3.5 4 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 PID Apertura=8.0 53 0 500 1000 1500 2000 2500 3000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 0 500 1000 1500 2000 2500 3000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 Apertura=0.0 On/Off 0 500 1000 1500 2000 2500 3000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 0 500 1000 1500 2000 2500 3000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 PID Apertura=8.0 0 2 4 6 8 10 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 0 2 4 6 8 10 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 Apertura=0.0 On/Off 0 2 4 6 8 10 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 0 2 4 6 8 10 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 PID Apertura=8.0 54 0 500 1000 1500 2000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 0 500 1000 1500 2000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 Apertura=0.0 On/Off 0 500 1000 1500 2000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 0 500 1000 1500 2000 0 20000 40000 60000 80000 100000 120000 0 1 2 3 4 5 6 7 8 9 10 Caudal (m3/s) Apertura compuerta (m) Tiempo(s) Caudal a la entrada Caudal a la salida Compuerta 1 Compuerta 2 Compuerta 3 PID Apertura=8.0 0 1 2 3 4 5 6 7 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 0 1 2 3 4 5 6 7 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 Apertura=0.0 On/Off 0 1 2 3 4 5 6 7 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 0 1 2 3 4 5 6 7 0 20000 40000 60000 80000 100000 120000 Setpoint 1 Nivel 1 Setpoint 2 Nivel 2 Setpoint 3 Nivel 3 PID Apertura=8.0 55 Las conclusiones se detallan en el cap´ıtulo siguiente, aunque solamente viendo las gr´aficas podemos afirmar que: La habilitaci´on de zonas de inundaci´on controlada para intentar reducir el pico de caudal en las avenidas es un m´etodo muy eficiente para la gesti´on de las avenidas, puesto que consigue laminar el hidrograma entrante de una manera considerable. En el caso de caudales bajos, se aprecia que no es imprescindible la instalaci´on de ning´un algoritmo de regulaci´on, ya que el hidrograma laminado alcanza pr´acticamente el mismo pico de caudal que si hubiese instalado alg´un mecanismo regulador. La laminaci´on que producen los mecanismos de regulaci´on On/Off y el controlador PID son muy parecidos, aunque con caudales altos se aprecia una mejor´ıa en el amortiguamiento del pico si actuamos con un controlador PID. Si tenemos un hidrograma entrante cuyo pico de caudal es muy grande y adem´as se prolonga mucho en el tiempo, como ocurre en el hidrograma no2, la diferencia entre aplicar un controlador y dejar que el agua entre por s´ı sola a las ´areas de inundaci´on se hace visible. Adem´as pudiera darse el caso que tuvi´esemos un hidrograma con esas caracter´ısticas que adem´as tuviese dos picos. Si dej´asemos entrar libremente el agua durante el primer pico de caudal y se llenan las zonas de inundaci´on, cuando llegue el segundo pico ya no tenemos margen de laminaci´on. Sin embargo fijando un nivel de referencia elevado para el primer pico (con cualquiera de los dos tipos de regulaci´on) y dejando actuar a las ´areas de inundaci´on controlada solamente en el segundo pico conseguimos una mejor gesti´on de este tipo de avenidas. 56 Cap´ıtulo 6 Conclusiones y trabajo futuro 6.1. Conclusiones La Hidr´aulica Computacional es uno de los campos de la ciencia en el cual la irrupci´on de los ordenadores exige una nueva manera de trabajar, aunque sin olvidar los desarrollos te´oricos y los m´etodos experimentales. En particular, la simulaci´on num´erica de ondas de avenida se esta convirtiendo en una herramienta muy potente y eficiente para la Administraci´on, aunque, cada vez m´as, la demanda de una mayor precisi´on en la descripci´on del modelo y en los resultados posteriores y un coste computacional no muy elevado se han convertido en exigencias muy frecuentes. Por ello, la investigaci´on en esta rama es y ser´a una parte esencial en el desarrollo de nuevas t´ecnicas para la gesti´on de las avenidas. Si atendemos al an´alisis dimensional realizado, podemos concluir: •A igual volumen del hidrograma e igual tiempo que permanece abierta la compuerta, el instante de apertura influye, de manera que si se abre la compuerta antes del pico de caudal o despu´es, el amortiguamiento no es ´optimo. Por otro lado, cuanto mayor sea el tiempo que permanece abierta la compuerta, mayor es el abanico de tiempos “validos” para un aprovechamiento ´optimo de las zonas de inundaci´on controlada. •Si el instante de apertura es fijo en el tiempo, al igual que el volumen del hidrograma, el tiempo que permanece abierta la compuerta es muy importante: cuanto m´as tiempo dejemos abierta la compuerta, m´as agua entrar´a al dep´osito, con lo que el amortiguamiento ser´a mayor. •La influencia que tiene el volumen del hidrograma en nuestro sistema es clara: si el tiempo que permanece abierta la compuerta es el mismo, y el instante 57 [40] G. Rosatti, J. Murillo, L. Fraccarollo. Generalized Roe schemes for 1D two-phase, free-surface flows over a mobile bed. Journal of Computational Physics 54, 543–590, 2007. [41] J.C. Rutherford. River Mixing. (Wiley, New York, 1994),p. 21, 1994. [42] P. Scotton y A. Armanini. Experimental investigation of roughness effects of debris flow channels. 6th Workshop on two-phase flow prediction. Erlangen, 1992. [43] J.J. Stoker. Water waves. Intersc. Pub. Inc., New York, 1957. [44] V.L. Streeter, E.B. Wylie y K.W. Bedford. Mec´anica de fluidos McGraw- Hill, 1999. [45] E.F. Toro. Riemann solvers and Numerical Methods for fluid dynamics: A practical introduction. Springer-Verlag, Berlin, 1997. [46] E.F. Toro. Shock-Capturing Methods for Free-Surface Shallow Flows. (Wiley, New York, 2001),p. 109, 2001. [47] M.E. V´ azquez-Cend´ on. Improved treatment of source terms in upwind schemes for the shallow water equations in channels with irregular geometry. Journal of Computational Physics 148, 497–498, 1999. [48] C.B. Vreugdenhil. Numerical methods for shallow-water flow. Kluwer Ac. Pub., Dordrecht, The Netherlands, 1994. [49] R.C. Ward y M. Robinson. Principles of hydrology. McGraw-Hill, 1990. [50] F.M. White. Mec´anica de Fluidos. McGraw-Hill, 1979. [51] Yun Li, K.H. Ang, G.K.H. Ang y G.C.Y. Chong PID Control System, Analysis and Design Glasgow ePrints Service 64