Ensayo comparativo entre modelos unidimensionales y bidimensionales en la modelización de la rotura de una balsa de materiales sueltos erosionables

Resumen

El desarrollo de ordenadores cada vez más potentes y las recientes investigaciones en el campo de la hidroinformática hace posible la modelización del flujo en lámina libre en dos dimensiones producido por una rotura de una presa de materiales sueltos. En el presente trabajo se comparan los resultados obtenidos mediante el modelo unidimensional HEC-RAS y el bidimensional CARPA. La comparativa se hace analizando el efecto que tienen las hipótesis simplificadoras implícitas en el uso de la plataforma HEC-RAS de unidimensionalidad, no infiltración y existencia de flujo inicial, sobre los hidrogramas de caudal en una sección testigo situada aguas abajo del eje drenante. Los modelos utilizados reproducen el flujo de agua generado por una posible rotura de la balsa número 1 del sector 5 del canal Segarra- Garrigues en la cuenca del río Ebro en España.

Abstract

The development of more and more potent computers and the recent research in the field of hidroinformatics makes possible the free surface flow modelling in two dimensions caused by earthen dam failures. In this paper, the results obtained by uni-dimensional model (HEC-RAS) and two-dimensional model (CARPA) are compared. The use of the HEC-RAS software assumes the hypothesis of unidimensionality to be true, no infiltration and existence of a minimal initial flow. The comparison is made by analyzing the effect of these hypothesis in the downstream flow hydrographs. The used models reproduce the water discharge generated by a possible failure of dam number 1 of the 5th District of the Segarra-Garrigues Irrigation Project in the Ebro river basin in Spain.

1

Introducción

La aplicación de la normativa vigente de clasificación de una presa de un gran embalse en función del daño que puede realizar su posible rotura, lo que se denomina «clasificación en función del riesgo potencial», en la clasificación de una pequeña balsa de regadío, puede llevar a la conclusión de que un agricultor deba implementar un Plan de Emergencia. De entrada, el sentido común nos dice que ello parece un tanto extraño.

Algunas de las causas de dicha situación pueden encontrarse en la propia modelización hidráulica de la avenida provocada por la rotura. Como siempre pasa en hidroinformática, cuanto más detallado es el modelo utilizado, mejor es el conocimiento que se obtiene pero más difícil es su construcción y mayor es su complejidad. La utilización de modelos excesivamente simplificados era obligada porque hace unas décadas no existían ordenadores suficientemente potentes para el uso de modelos sofisticados. Pero hoy en día, la tecnología informática disponible permite el uso de modelos cada vez más complejos y acordes con la realidad, de tal manera que el establecimiento de tales hipótesis simplificadoras resultan injustificables.

A lo largo de la historia han ido apareciendo hipótesis simplificadoras que en la actualidad no tienen por qué ser establecidas:

• –

  Régimen uniforme.

• –

  Régimen estacionario y suavemente variable.

• –

  Régimen no estacionario y suavemente variable.

• –

  Unidimensionalidad.

• –

  Bidimensionalidad.

• –

  Altamente variable (presión no hidrostática).

• –

  Sin pérdidas por infiltración.

• –

  Continuidad espacial del coeficiente de rugosidad.

• –

  Existencia de un flujo mínimo como condición inicial.

• –

  Sin transporte de sedimentos.

• –

  Propiedades reológicas del fluido constantes.

En este trabajo, se presenta un ensayo comparativo entre dos aproximaciones numéricas distintas de cálculo del flujo de agua en lámina libre que se produce ante una rotura de una balsa constituida de materiales sueltos. El objetivo es poner de manifiesto el efecto sobre la clasificación que tiene el establecimiento de algunas de las hipótesis simplificadoras. Concretamente, la unidimensionalidad, la inexistencia de infiltración y la existencia de un flujo mínimo.

A pesar de que los sedimentos generados por la rotura pueden producir cambios significativos en las propiedades reológicas del agua, en el presente ensayo comparativo no se realiza ningún tipo de análisis de sensibilidad a este parámetro ni se analiza el efecto que produciría sobre la clasificación la introducción de esta hipótesis.

Después de esta introducción, se describe la balsa que ha servido de ejemplo para el estudio, se da la descripción matemática de los esquemas numéricos utilizados, se plantea el conjunto de ensayos para el estudio de sensibilidad a los factores de forma de la balsa que fueron planteados durante su diseño y que fueron considerados como criterios de proyecto, se muestra el comparativo de resultados cuando se supone cierta la hipótesis de unidimensionalidad y cuando no, y se muestra el efecto de introducción de algún tipo de infiltración.

1.1

Balsa objeto del estudio

La balsa objeto de este estudio de sensibilidad es la balsa de la margen derecha del Sector número 5 del sistema del Canal Segarra-Garrigues en la provincia de Lleida en la cuenca del río Ebro en España. En el momento de realización de este estudio de sensibilidad, la balsa se encontraba en fase de proyecto de manera que el resultado de la clasificación en función del riesgo potencial debía de ser considerada como criterio de diseño: el diseño tenía que ser tal que su clasificación diera de «Clase C». Este tipo de clase se corresponde con las presas cuya rotura o funcionamiento incorrecto puede producir daños materiales de moderada importancia y solo incidentalmente pérdida de vidas humanas [1]. La situación de la balsa se muestra en la figura 1, así como el conjunto de ejes de desagüe ante una posible rotura.

Figura 1.

Plano situación de la Balsa 1 de la margen derecha del Sector número 5 del Sistema del Canal Segarra-Garrigues en la provincia de Lleida en la cuenca del río Ebro en España. El plano también indica el conjunto de ejes y su nomenclatura que a primera vista tuvieron que ser analizados.

En la figura 1 puede verse la situación del núcleo urbano de Claravalls y los posibles ejes drenantes. El núcleo urbano es cruzado por el Eje A y se encuentra a unos 900 m de la balsa siguiendo el curso del Eje A.2. Como posteriormente se pudo comprobar, la sección en la entrada a la población que se corresponde con el punto de confluencia entre el Eje A.1 y el Eje A.2, y que denominaremos Sección Testigo, fue la sección crítica para la clasificación porque de las condiciones hidrodinámicas del flujo que se producían en este punto dependía el «daño» producido que decantaba la clasificación hacia la «Clase C» o a otro tipo.

El análisis detallado que genera un determinado hidrograma de caudal en términos de velocidad y calado a través de la geometría concreta de la población está fuera del alcance de este trabajo. Sirva decir solamente que el Eje A a su paso por la población se encuentra canalizado debido a la presencia de otra balsa existente de menor entidad que la analizada aquí.

Como puede verse en la figura 2, de forma cuadrangular, la balsa está ubicada en el extremo de una altiplanicie y se encuentra parcialmente enterrada en terrenos rocosos en profundidad.

Figura 2.

Detalle de la situación geomorfológica del entorno de la ubicación de la balsa objeto de este estudio.

2

Descripción de los esquemas numéricos

2.1

Esquemas utilizados

En el presente trabajo se utilizaron 2 herramientas de modelización distintas: el programa HEC-RAS desarrollado por el Hydrologic Engineering Center de los Estados Unidos [2] y el programa CARPA [3].

El primero es una herramienta de amplia difusión que resuelve las ecuaciones de Saint Venant en una dimensión por el clásico esquema de Preissmann [4], esquema que es inestable en situaciones de cambios de régimen o discontinuidades, como en el caso presente de un frente de onda, en cuyo caso se utiliza el método denominado LPI (Locar Partial Inertia) [5], desarrollado por Fread [6] afectando los términos de inercia de la ecuación del momento, en cierta manera añadiendo una difusión artificial, y por lo tanto alterando el comportamiento del fluido.

El programa CARPA resuelve las ecuaciones de Saint Venant en dos dimensiones que pueden escribirse de forma conservativa como:

donde U es el vector de variables de flujo, F es el tensor de flujo y H es el termino independiente o término fuente, que responden a las expresiones:

donde t es el tiempo, g la gravedad, h es el calado, u la velocidad en la dirección del eje de abcisas y v la velocidad en la de las ordenadas, S0x es la pendiente del fondo en la dirección de las abcisas y S0y en el de las ordenadas y Sfx es la pendiente de rozamiento en la dirección de las abcisas y Sfy en el de las ordenadas:

donde m es el coeficiente de rugosidad de Manning.

Siguiendo un discretización en volúmenes finitos explícita en el tiempo, que se ha escogido porque es especialmente adecuada a la hora de desarrollar métodos conservativos, se obtiene:

donde Fi,wl* es el flujo numérico normal a una pared de un volumen finito, wl representa el índice correspondiente a la l-ésima pared del elemento i, N_i el número de lados del elemento y S_i su superficie proyección en planta, el vector n es la normal exterior a la pared wl del elemento i, li,wl es su longitud y n representa el índice de evolución temporal a incrementos temporales de longitud ▵t que depende de la condición de Courant.

Para discretizar los términos de flujo se utilizan esquemas conservativos descentrados de tipo Godunov, concretamente el esquema descentrado de Roe en orden 2 de precisión en espacio. Esta formulación conservativa de las ecuaciones proporciona una buena resolución de los choques transónicos que se puedan producir en la solución, por lo que es un método recomendado para modelizar resaltos hidráulicos, rotura de presas y ondas de choque.

2.2

Esquema descentrado de Roe de orden 1

El flujo numérico se puede expresar como:

donde j indica el elemento que conecta con el i a través de la pared wl, donde los valores propios y vectores propios del jacobiano aproximado del tensor de flujo responden a:

con los coeficientes:

y los estados promediados en cada contorno según:

El esquema anterior es de orden 1 en espacio ya que el descentramiento del flujo convectivo es equivalente desde el punto de vista matemático a añadir un término de difusión.

2.3

Extensión del flujo numérico a orden 2

El esquema anterior es de orden 1 debido a la difusión numérica introducida en la discretización del flujo convectivo. A pesar de ello es utilizado de forma habitual en códigos de CFD como esquema por defecto, debido a su estabilidad numérica. Cuando se requiere un orden de precisión elevado con un tamaño de malla que no sea excesivamente fino, es necesario recurrir a esquemas de orden superior. En el módulo hidrodinámico de CARPA se consigue aumentar el orden de precisión del esquema de Roe mediante una extensión de orden 2 del flujo numérico, y una limitación TVD (Total Variation Diminishing) del mismo.

con

donde d_ij es la distancia existente entre los centroides de los elementos contiguos i y j. Se utiliza por defecto la función limitadora Minmod [7]:

En CARPA se utiliza de forma análoga una discretización centrada de todos los términos fuente excepto del término fuente pendiente del fondo. El principal motivo de utilizar una discretización descentrada de la pendiente del fondo frente a una discretización centrada es que se calcula de forma exacta la solución hidrostática con batimetría irregular, evitando de esta forma la aparición de oscilaciones espurias en la superficie libre del agua y en las velocidades. Estas oscilaciones son en general pequeñas, pero pueden llegar a ser de magnitud considerable en problemas con batimetrías irregulares.

CARPA ha sido verificado mediante un gran número de casos en [8] y [9], algunos de los cuales basados en soluciones analíticas. Cabe destacar también que CARPA presenta un novedoso tratamiento del término fuente que evita la difusión numérica del frente de onda evolucionando en un lecho seco. Algunos de los datos de tipo computacional de los ensayos numéricos han sido:

—Modelo unidimensional HEC-RAS:

• –

  Número de nodos del Eje 1 con una sección cada 2 m: 2.383.

• –

  Tiempo de resolución: algunos minutos en función del eje.

• –

  Tiempo de simulación: 3 horas.

• –

  Incrementos de tiempo: se ha limitado hasta el máximo permitido por el programa el paso de tiempo a 1 s.

—Modelo bidimensional CARPA:

• –

  Número de volúmenes finitos: 51.276.

• –

  Tiempo de resolución: unos 50 minutos (con el software y el hardware de 2006).

• –

  Tiempo de simulación: 3 horas.

• –

  Incrementos de tiempo: dependiendo de la condición de Courant y la homogeneidad de los volúmenes finitos «mojados», los incrementos de tiempo han sido del orden de 0, 01s.

3

Análisis de sensibilidad a los parámetros de forma

Como criterios de diseño de partida se tenía:

• –

  El volumen de almacenamiento: 497.548 m³.

• –

  La forma en planta: cuadrangular con talud de 2:1 (H:V), esto es, no adaptada a la forma del terreno.

• –

  Se permitía ajustar la posición de la balsa de manera que se podía jugar con la profundidad de agua en la balsa enclavada en el terreno rocoso.

Así pues, se decidió llevar a cabo un análisis de sensibilidad a los 2 parámetros que determinan el hidrograma de caudal generado por la rotura, que son la altura máxima de la lámina H sobre cimientos y el volumen V cercado que es capaz de generar escorrentía, ambos correlacionados a través del factor de forma de la balsa.

Estos parámetros establecen el modelo de formación de la brecha en balsas de materiales sueltos erosionables según se especifica en las recomendaciones de la guía técnica [1]. Esta guía recomienda un modelo lineal de progresión de la brecha que se supone de forma trapezoidal. El modelo de progresión recomendado se define a partir de 2 variables que dependen de los parámetros anteriormente citados que son el tiempo de formación y el ancho medio de la brecha. Estas variables se definen de la siguiente manera:

donde T(h) es el tiempo de formación de la brecha, V(Hm³) es el volumen cercado en terraplén, H(m) es la altura de la lámina de agua en la balsa hasta cimientos y B(m) es el ancho medio de la brecha.

Como ya se ha dicho anteriormente, la balsa se ubica en una zona rocosa que requiere llevar a cabo voladuras para poderla enclavar. Por consiguiente, tuvo que hacerse un análisis de sensibilidad a los parámetros de forma. Concretamente, se analizaron los casos mostrados en la tabla 1.

Tabla 1.

Resumen de los ensayos realizados en el análisis de sensibilidad a los parámetros de forma

V	H	T	B
0,003	0,1	5,42	2,25
0,035	0,5	1,63	7,44
0,066	1,0	1,17	10,26
0,096	1,5	0,96	12,44
0,127	2,0	0,83	14,28
0,153	2,5	0,75	15,74
0,185	3,0	0,68	17,35
0,214	3,5	0,63	18,68
0,242	4,0	0,58	19,91
0,270	4,5	0,55	21,07

Mediante un modelo unidimensional, sin contemplar la infiltración y con un caudal mínimo de generación de flujo se ensayaron todos los casos indicados en la tabla 1 a lo largo de todos los ejes a estudiar. Algunos de los resultados obtenidos en el ensayo para el Eje A.2 pueden verse en la figura 3.

Figura 3.

Hidrogramas de caudal de paso por la Sección Testigo obtenidos con un Modelo Unidimensional en tres casos de los ensayados para el Eje A.2 cuando H = 4,5 m (la brecha llega a cota 403), H = 3,5 m (la brecha llega a cota 404) y H = 2,5 m (la brecha llega a cota 405).

La respuesta de la propia geometría en el punto de conexión entre el Eje A.1 y el Eje A.2 en la Sección Testigo a la entrada de la población hace que un caudal inferior a los 30 m³/s genere una clasificación de «Clase C». Por lo tanto, se considera como óptima una profundidad de agua en la balsa de 2, 5 m (que se corresponede con la llegada de la brecha hasta la cota 405 m), según se desprende de la gráfica de la figura 3, con lo que H = 2, 5 m fue el valor adoptado en el diseño de la balsa. Ello se consiguió ubicando la balsa en un lugar dónde el terreno natural del perímetro perforado de la balsa estuviera como mínimo a dicha cota.

Llegados a este punto, ya estamos en disposición de analizar la influencia de algunas de las hipótesis anteriormente mencionadas.

4

Hipótesis de unidimensionalidad y flujo mínimo inicial

La hipótesis de que el flujo es en una dimensión, supone que todo el caudal generado por la rotura de la balsa debe seguir un solo camino, el de los denominados Ejes A.1, A.2, A.3 y B de la figura 1. Entonces, siendo coherentes con esta hipótesis, se definió una canalización artificial que conducía todo el caudal desde la balsa a la entrada del cada uno de los ejes.

La problemática presentada en el párrafo anterior puede clarificarse aún más si se observa el sitio concreto donde se ha construido la balsa y que se muestra en la figura 4.

Figura 4.

Zona elegida para la construcción de la Balsa.

Si se observa con detalle la fotografía de la figura 4, se hace difícil elegir el camino que tomaría el agua en caso de rotura. La única forma de saberlo es la utilización de un modelo de flujo en lámina libre en dos dimensiones cuyos caminos son elegidos por el propio flujo y son solución del modelo y no una imposición al modelo. Además, la dirección tomada por el agua depende claramente del lugar de la balsa donde se produce la rotura. Analizada la información geométrica de la balsa encajada en el territorio, de manera que la profundidad de agua en la balsa máxima fuera de 2, 5 m, se encontraron tres posibles puntos de rotura denominados aquí Rotura Norte (figura 5), Rotura Oeste (figura 6) y Rotura Sur (figura 7).

Figura 5.

Rotura Norte.

Figura 6.

Rotura Oeste.

Figura 7.

Rotura Sur.

Para cada hipótesis de punto de rotura se estudió el flujo producido por la misma rotura mediante el modelo 2D. En el «Apéndice: Resultados de las simulaciones 2D»adjuntado al final de este trabajo, se dan las secuencias del avance de la lámina de agua de cada caso y en ellas puede verse como (figs. 11–32):

• –

  En los tres casos, el caudal generado por la rotura se reparte entre los 3 ejes posibles, dependiendo la proporción de reparto, del lugar de rotura.

• –

  En los tres casos, el Eje B es el que recibe la menor proporción de caudal y al tratarse de un eje que no era crítico en cuanto al posible daño, se descartó continuar su análisis.

• –

  Gracias a la aportación de la Rotura Oeste y de la Rotura Sur al Eje B de parte de caudal, dichas roturas pasaron a ser secundarias en cuanto a su daño potencial en comparación con la Rotura Norte.

• –

  Como consecuencia de los puntos anteriores, el flujo producido por la Rotura Norte a través de los ejes A.1 y A.2 es el que puede infringir el máximo daño a la población.

• –

  Justo antes de la entrada del caudal en la población, los dos Ejes A.1 y A.2 vuelven a unirse y lo hacen justo aguas arriba de la población en la Sección Testigo. Con ello se ve reducido el efecto de laminación del caudal punta y su efecto sobre la clasificación.

Desde el punto de vista de la clasificación, uno de los resultados que pueden resumir mejor el efecto de la aplicación de la hipótesis de unidimenisonalidad es el comparativo de hidrogramas de llegada a la Sección Testigo a la entrada del núcleo urbano (figura 8).

Figura 8.

Comparativo de hidrogramas de caudal en la Sección Testigo situada a la entrada del encauzamiento en el núcleo urbano de Claravalls. Se compara el Modelo Unidimensional (Eje A.2) con el Bidimensional (Rotura Norte). Ambos modelos comparten las mismas hipótesis, principalmente el mismo volumen y mismo grosor de la lámina almacenada, pero no la dimensionalidad y el caudal inicial.

Por otro lado, existe una gran sensibilidad al punto de rotura de la balsa como indican los hidrogramas de la figura 9.

Figura 9.

Comparativo de hidrogramas de caudal generados para el caso de la Rotura Norte y la Rotura Oeste.

A la vista de las gráficas de las figuras 8 y 9, cabe destacar los siguientes puntos:

• –

  Una reducción del 17 % del caudal punta. Como ya se ha dicho anteriormente, si la Sección Testigo se hubiera situado aguas arriba de la confluencia entre los Ejes A.1 y A.2, este efecto laminador hubiera sido mucho más evidente, de manera que, muy posiblemente, la Rotura Norte hubiera dejado de ser la crítica.

• –

  La existencia de cierto caudal inicial totalmente necesario desde el punto de vista numérico en la aplicación de la mayoría de modelos unidimensionales provoca que la llegada del «frente de onda» a la población sea mucho más rápido (curva con más pendiente, más vertical) porque, como es sabido, el agua fluye más facilmente en un lecho con agua que en un lecho seco.

• –

  La gran mayoría de balsas a modelar no tienen cuenca, con lo que, en estos casos, el caudal mínimo inicial es una hipótesis que no debería establecerse.

• –

  La Rotura Norte presenta un caudal punta superior con lo que pasa a ser la crítica.

Como conclusión parcial de este apartado puede decirse que, la aplicación a ultranza de la hipótesis de unidimenisonalidad para la defensa de un determinado modelo, magnifica el daño potencial y su transposición a la clasificación resulta más pesimista. Desde nuestro punto de vista, la utilización de modelos en 2D hoy en día no reviste ningún inconveniente, tanto desde el punto de vista de la capacidad actual de los ordenadores ni desde la complejidad de su utilización. Por lo tanto, recomendaríamos utilizar solamente este tipo de modelos para la clasificación de pequeños embalses.

5

Hipótesis de inexistencia de pérdidas

La introducción en el modelo bidimensional de cierta capacidad de infiltración en toda la superficie de forma homogénea en todo el dominio también tiene un efecto laminador del hidrograma de caudal circulante por la Sección Testigo a la entrada del núcleo urbano. En la figura 10 se muestra dicho efecto. Concretamente, se introdujo con CARPA el modelo de infiltración lineal que después de los 15 mm primeros de infiltración acumulados, la tasa de infiltración permanecía constante a lo largo del tiempo a razón de 5 mm/h.

Figura 10.

Comparativo de hidrogramas de caudal de la Rotura Norte en la Sección Testigo.

A parte del efecto laminador de la infiltración, también reduce la longitud de los ejes a estudiar puesto que si un modelo es unidimensional y no tiene infiltración, su caudal puede generar una superficie mojada enormemente larga. Nuestra experiencia nos dice que puede llegar a tener decenas de kilómetros, cosa faltada de sentido común.

A la vista de las curvas de la figura 10 cabe destacar que la utilización de un modelo bidimensional con cierta capacidad de infiltración produce una reducción del 30% del caudal punta generado por el modelo unidimensional sin infiltración.

6

Conclusiones

• 1.

  La profundidad de agua en la balsa hasta cimientos H es el factor más sensible en la clasificación presentada aquí. Dado que la «Clasificación de tipo C»se ha tomado como un criterio de diseño de la balsa, dicha profundidad ha tenido que ser optimizada en función del daño potencial.

• 2.

  Con un Modelo 2D se puede simular mejor el flujo en lámina libre y por lo tanto, se puede ajustar mejor la clasificación, reduciendo costes económicos: en el ejemplo en estudio se ha reducido un 17% el caudal punta de paso por la Sección Testigo a la entrada del núcleo urbano, con lo que se hubiera podido aumentar el calado de la balsa y reducir la profundidad del enclavamiento y el consiguiente volumen de voladuras.

• 3.

  En geometrías planas, plantear la hipótesis de unidimensionalidad es falsa y puede falsear enormemente los resultados de una clasificación. Con un modelo bidimensional, se puede encontrar el comportamiento del flujo en términos de reparto de caudales entre barrancos cuando la balsa se sitúa en zonas planas y en un máximo relativo.

• 4.

  Como consecuencia del punto anterior, recomendaríamos el uso exclusivo de modelos bidimensionales para el estudio de balsas de materiales sueltos y erosionables puesto que nuestra experiencia nos dice que la unidimensionalidad es la excepción y no la regla.

• 5.

  La infiltración es un factor clave de laminación de hidrogramas porque afecta a la generación de un lecho de agua inicial por donde el frente de onda evoluciona más rápidamente y por consiguiente, se lamina menos. Con los modelos bidimensionales se puede introducir algún tipo de infiltración a partir de cobertura SIG.

• 6.

  Dada la importancia del factor de infiltración se recomendaría un estudio previo de los parámetros que caracterizan la infiltración mediante infiltrómetros o cualquier otra técnica.

Financiación

Este estudio ha sido financiado en parte por la empresa Ingeniería E. Inalba ubicada en la población de Mollerussa en la provincia de Lleida.

Agradecimientos

Agradecemos también al personal de dicha empresa su predisposición y apoyo en la elaboración de estudios como el presente.