Latest revision as of 08:06, 18 April 2017

Resumen

La circulación del agua a través de medios porosos gruesos, como la escollera, la grava o incluso la arena de mayor tamaño, tiene la particularidad de exhibir una relación no lineal entre el gradiente hidráulico y la velocidad de filtración. Uno de los problemas de ingeniería en que se observa esto es cuando se produce un sobrevertido en una presa de escollera y el agua, pasando sobre su coronación, se infiltra en el espaldón y circula por su interior. La modelación numérica del problema se resuelve aquí mediante un esquema de diferencias finitas, que se implementa en el código MNLEE.RDS para su aplicación al sobrevertido en presas de escollera. Además, se proporciona un procedimiento para resolver el acoplamiento del movimiento del agua en lámina libre sobre el talud de la presa y de filtración por el interior del espaldón, mediante la introducción del concepto «ley de intercambio de caudales».

Abstract

Flow through coarse porous media, as rockfill, gravel or sand of big size, exhibit a non linear relation between hydraulic gradient and seepage velocity. This behavior is observed when water flows over the crest of a rockfill dam, seeps into the downstream shoulder and moves through the rockfill voids. A finite difference scheme is used here to solve this problem through numerical modeling. It is implemented in the code MNLEE.RDS for solving dam overtopping problems. In addition, a procedure is provided to solve the coupled problem of flow over and through the rockfill dam. For that purpose, the concept of «law of flow exchange»is introduced.

Palabras clave

Filtración ; Medios porosos ; Diferencias finitas ; Escollera ; Presa ; Sobrevertido

Keywords

Seepage ; Porous media ; Finite differences ; Rockfill ; Dam ; Overtopping

1. Introducción

Al realizar estudios de filtración del agua a través de medios porosos finos, como el núcleo de arcilla de una presa de materiales sueltos, se considera una relación lineal entre el gradiente de carga y la velocidad de filtración, de acuerdo con la conocida ley de Darcy. La experiencia ha demostrado que en medios granulares gruesos (escolleras, gravas y arenas gruesas) esta ley no es válida, pues dicha relación no es lineal. En este sentido podemos hablar de medios porosos no lineales. Como ejemplo de ello podemos citar la filtración en acuíferos de grano grueso o en materiales granulares utilizados como elemento de drenaje. También es el caso de la circulación de agua que se establece por el interior del espaldón de una presa de escollera cuando excepcionalmente, debido a la insuficiencia de capacidad de su aliviadero, durante una avenida el agua llega a pasar sobre su coronación. Actualmente, las distintas administraciones hidráulicas de gran número de países están abordando el estudio y mejora de la seguridad de sus presas, teniendo en este contexto una importancia singular el problema del sobrevertido en presas de materiales sueltos, pues ha sido causa de gran número de roturas en este tipo de presas. Los principales mecanismos involucrados en la rotura de una presa de escollera por sobrevertido son el deslizamiento en masa y el arrastre de partículas. Además, existe el peligro de que se produzca migración interna de partículas. Por tanto, para evaluar la seguridad o para diseñar una presa resistente a este fenómeno, es condición previa conocer las presiones intersticiales y los gradientes hidráulicos y, por tanto, resolver el problema de filtración en un medio granular no lineal. Para ello, es necesario conocer la fórmula de resistencia que expresa la antes referida relación gradiente/velocidad y, a partir de ella, establecer una ecuación de campo que caracterice el fenómeno de filtración en este tipo de medios granulares gruesos. Las condiciones de frontera particularizarán el problema a resolver, de modo que el problema quede adecuadamente planteado. La resolución por métodos numéricos implica la discretización del dominio de filtración y la resolución del sistema de ecuaciones resultante. Se expone en apartados sucesivos cada uno de los aspectos señalados. Debe llamarse la atención sobre la necesidad de que la relación entre el tamaño del dominio de filtración y el tamaño de partícula sea suficientemente alta para asegurar la aplicabilidad de los procedimientos diferenciales que aquí se emplean. Esta consideración, que siempre aplica, puede resultar condicionante en el caso de que el tamaño de grano del material poroso sea grande, como sucede en la escollera.

2. Fórmula de resistencia

Diversos autores han propuesto leyes alternativas para medios no lineales, siendo las más usadas las leyes cuadrática y exponencial.

Ley cuadrática:

Ley exponencial:

Siendo i el gradiente hidráulico y V_d la velocidad de filtración. Esta última se define como la relación entre el caudal filtrado y la sección de paso total, incluidos huecos y esqueleto sólido.

George y Hansen (1992)[1] obtuvieron unas fórmulas que permiten determinar los coeficientes de la fórmula cuadrática a partir del coeficiente y el exponente de la fórmula exponencial y viceversa.

Conversión de ${\textstyle i=c\quad V_{d}^{m}}$ a ${\textstyle i=a\quad V_{d}+b\quad V_{d}^{2}}$

Siendo i = 0 para V_d = 0 y V_{d ,max} el valor máximo de la velocidad de filtración en el intervalo considerado.

Conversión de ${\textstyle i=a\quad V_{d}+b\quad V_{d}^{2}}$ a ${\textstyle i=c\quad V_{d}^{m}}$

La modelación se ha realizado partiendo de una fórmula de resistencia de tipo exponencial. No obstante, si se dispone de una fórmula de tipo cuadrático para un material determinado, basta con la conversión mediante las anteriores fórmulas de George y Hansen. A partir de resultados obtenidos empíricamente, Parkin, Trollope y Lawson (1966)[2] definieron un nomograma que permite obtener el coeficiente c de la fórmula exponencial a partir del índice de huecos y la superficie específica del medio granular considerado. Así mismo, llegaron a la conclusión de que el exponente de dicha fórmula varía poco en torno al valor 1,85 próximo al valor 2 correspondiente a movimiento totalmente turbulento, y alejado del valor 1 correspondiente a movimiento laminar. La causa de que el exponente de la fórmula exponencial tenga un valor próximo a 2, pero sin llegar al mismo, ha sido objeto de controversia. Con independencia de cuál sea la explicación, es un hecho comprobado que el valor de dicho exponente se encuentra, en efecto, en dicho rango de variación.

3. Movimiento del agua. Planteamiento general del problema

Una vez conocida la fórmula de resistencia, puede plantearse el estudio del movimiento del agua a través de la escollera. El análisis diferencial conduce a un modelo matemático diferencial que se plantea como un problema de campo. El análisis se limitará al caso de movimiento permanente, bidimensional de un fluido viscoso e incompresible, el agua, en un medio poroso homogéneo e isótropo. En primer lugar es preciso definir el problema. Para ello debe conocerse la ecuación de campo que caracteriza el fenómeno físico: movimiento del agua a través de un medio poroso determinado, y las condiciones de frontera, que particularizan el problema para un caso concreto. Como resultado final del cálculo, se obtiene un campo de presiones y un campo de velocidades, si bien la ecuación de campo puede plantearse en términos de otras magnitudes, típicamente su potencial. En cualquier caso, una vez resuelto el problema, es posible conocer la presión y la velocidad en cualquier punto del dominio de filtración y, a través de la velocidad, el caudal de agua filtrado. Por tanto, el planteamiento del problema consiste en la definición de la ecuación de campo y de las condiciones de frontera. Cuando hablamos de resolver el problema, nos referimos a la determinación de las presiones y velocidades en cada punto del medio. Se dice que el problema está bien planteado cuando la solución existe y es única. Por otra parte, al plantear un problema de campo es necesario garantizar que la ecuación de campo se cumple en todo el dominio. El incumplimiento de esta condición invalida el planteamiento del problema.

4. Ecuación de campo

Consideremos una partícula de medio poroso de dimensiones suficientemente reducidas en relación con el dominio de filtración para que los procedimientos diferenciales sean aplicables: dxdy , y unos ejes de coordenadas ortogonales (x , y ). Sea el vector velocidad de filtración en un punto:

siendo x e y vectores unitarios de dirección la de los ejes de coordenadas x e y respectivamente y sentido el definido como positivo. ( fig. 1 )

Fig. 1. Partícula diferencial del medio poroso.

El caudal de agua que entra en la partícula de medio poroso es:

El caudal que sale es:

Admitiendo que ni el agua ni los elementos sólidos sufren deformación volumétrica y que la partícula está saturada, el caudal de agua que entra en ella tiene que ser igual que el que sale (ley de conservación de la masa o principio de continuidad), de donde se obtiene:

4.1. Medios porosos lineales

Si se cumple la ley de Darcy, tendremos que:

siendo K el coeficiente de permeabilidad e i el gradiente hidráulico. Por tanto:

siendo:

Donde: z es la altura sobre el plano de referencia y p es la presión. En realidad, h es la carga hidráulica total, pero la carga de velocidad tiene en la práctica un valor muy bajo, despreciable frente al de la carga piezométrica. Sustituyendo (12) en (10) tenemos que:

En un medio homogéneo e isótropo el coeficiente de permeabilidad es constante en todo el dominio. Por tanto:

Esta es la conocida ecuación de Laplace, que gobierna el movimiento del agua en un medio poroso cuando se cumple la ley de Darcy. Se trata, por tanto, de la ecuación de campo para movimiento lineal.

Si consideramos el operador laplaciano:

la ecuación de Laplace puede expresarse:

4.2. Medios porosos no lineales

Diversos autores han desarrollado ecuaciones de campo para medios porosos no lineales basándose en una fórmula de resistencia de tipo exponencial o cuadrático.

Engelund (1953)[3] y otros, considerando una fórmula de resistencia cuadrática con término lineal nulo, han llegado a la siguiente ecuación de campo.

siendo: h carga piezométrica, V la velocidad de filtración y α la inclinación del vector velocidad. A continuación veremos el desarrollo completo de las ecuaciones de campo de Parkin y de Volker para movimiento no lineal.

4.2.1. Ecuación de campo de Parkin

La ecuación de campo de Parkin (1971)[4] se basa en una ley exponencial:

admitiendo que los coeficientes c y m son constantes en el dominio de filtración. Estos coeficientes no varían de forma significativa en el rango de velocidades que se dan en la práctica, por lo que la hipótesis resulta aceptable. De (18) se obtiene que:

Tanto la velocidad como la pérdida unitaria de energía se refieren a la dirección de avance del agua en el punto considerado (línea de corriente del movimiento en el medio poroso). En general, la relación puede expresarse en forma vectorial. Teniendo en cuenta que:

donde, despreciando la carga de velocidad queda:

Admitiendo que el coeficiente c es constante, puede expresarse:

y llamando:

Se tiene:

Por definición del operador gradiente:

Siendo x e y los vectores unitarios en las direcciones de los ejes cartesianos x e y respectivamente.

Sustituyendo en (24) :

Por otra parte, tenemos que:

Por tanto, sustituyendo en (26) :

En consecuencia, las componentes del vector velocidad

son:

Como hemos visto en el apartado anterior, la ecuación de continuidad aplicada a una partícula de medio poroso indeformable conduce a:

Derivando las expresiones (30) y (31) respecto a x e y respectivamente, obtenemos:

Sustituyendo las expresiones (33) y (34) en la ecuación de continuidad (32) , obtenemos la ecuación de campo de Parkin.

donde recordemos que:

Parkin, Trollope y Lawson (1966)[2] hacen notar que, supuesto constante el coeficiente c , la anterior ecuación de campo es equivalente a otra en la que ϕ se sustituya por h , para obtener directamente la carga piezométrica. En consecuencia, la carga piezométrica (h ), y por tanto la presión (p ), dependen únicamente del exponente (m ) de la fórmula de resistencia. La ecuación de campo (35) se reduce a la ecuación de Laplace, válida para movimiento lineal, cuando m  =1, o lo que es lo mismo, cuando se cumple la ley de Darcy.

5. Condiciones de frontera

El fenómeno de filtración del agua a través del espaldón de una presa de escollera sometida a sobrevertido da origen a un problema de frontera mixto con superficie libre. Consideremos en primer lugar el movimiento del agua a través de una presa homogénea (fig. 2 ). Las condiciones de frontera son:

• Tramo AB: la carga piezométrica es constante e igual a la altura de agua sobre la base (H ). Se trata de una condición de tipo Dirichlet.
• Tramo DE: análogamente, la carga piezométrica es en este caso H ′. También es una condición de tipo Dirichlet.
• Tramo AE: admitiendo que la base es impermeable, el caudal que circula a través de este tramo de frontera es nulo. Por tanto, se trata de una condición tipo Neumann.
• Tramo BC: se trata de una superficie libre, cuya situación se desconoce a priori. La presión en todos sus puntos debe ser nula (referida a la presión atmosférica) y no puede circular agua a través de la misma. Se trata de una condición tipo Cauchy.
• Tramo CD: la carga piezométrica viene dada por la altura sobre el plano de referencia, ya que la presión es nula por encontrarse en contacto con la atmósfera. La condición es de tipo Dirichlet. Si el agua supera la cota de coronación se producirá el sobrevertido. Cuando la altura de lámina vertiente (H_L ) sea pequeña, la circulación superficial o exterior, sobre la coronación y el talud de aguas abajo, se reducirá a unos pocos metros. A partir de un determinado punto solo se producirá circulación interna a través del cuerpo de la presa (fig. 3 a). A medida que aumenta la altura de lámina vertiente, el punto extremo de infiltración se desplazará hacia aguas abajo sobre el talud de la presa y el D lo hará hacia aguas arriba. Cuando la lámina sea suficientemente elevada, por ser el caudal de llegada al embalse netamente superior al que puede circular por el interior de la presa, se producirá movimiento superficial sobre toda la longitud del talud de aguas abajo ( fig. 3 b). La condición de frontera en la coronación y el paramento de aguas abajo puede plantearse como de potencial impuesto, pero desconocido a priori, que es preciso determinar mediante tanteos sucesivos de modo semejante a como se procede con los problemas de superficie libre. En realidad, se trata de determinar la geometría de la superficie libre del agua que circula sobre el talud de la presa. Por tanto, para determinar la situación de una superficie libre se sigue el siguiente proceso:
• 1.° Se estima la situación de la superficie libre.
• 2.° Se calcula la solución al problema, obteniéndose la presión en la superficie libre.
• 3.° Si la presión no es nula en todos los puntos de la misma, se modifica su geometría, elevándola en la zona donde es positiva y reduciendo su cota donde es negativa.
• 4.° Se vuelve a realizar el cálculo y el proceso se repite hasta que se consigue que la presión sea nula en todos los puntos de la superficie libre. En el caso de que exista circulación sobre el talud de la presa, la condición que debe cumplirse es la que se explica a continuación. Se consideran dos secciones verticales separadas una cierta distancia Δx , siendo:
• q_{1(2),e (i )}  : caudal que atraviesa la sección 1 (2) circulando exteriormente (interiormente) sobre el talud de la presa (a través del cuerpo de la presa).
• Δq_1,2 : caudal que pasa del exterior al interior de la presa (o viceversa, si el signo es negativo) entre las secciones 1 y 2. Aplicando el principio de continuidad al recinto de medio poroso RSVU (fig. 4 ), tenemos:

El movimiento del agua sobre el talud es de caudal variable, dado que existe un intercambio de caudales con el interior de la presa.

Fig. 2. Presa homogénea sin elemento impermeable.

Fig. 3. Fases de movimiento del agua sobre la presa a) Fase inicial b) Fase final.

Fig. 4. Caudales circulantes a través de las secciones paralelas.

El movimiento del agua sobre el talud y a través del cuerpo de la presa son fenómenos acoplados, no siendo posible en pura teoría calcular uno de ellos con independencia del otro. En efecto, para resolver el movimiento del agua a través de la presa es preciso conocer la carga piezométrica en el talud de aguas abajo, y por tanto el perfil del agua que circula sobre el mismo. Pero para determinar este perfil tenemos que conocer el caudal que se añade o sustrae en cada punto del talud, que se determina mediante la expresión (36) , previa resolución del movimiento del agua a través de la presa.

En consecuencia, debe seguirse un proceso iterativo.

• 1.° Estimar una ley de variación del caudal circulante exteriormente.
• 2.° Modelar el movimiento del agua sobre la presa, determinando la carga piezométrica y el potencial a lo largo de la frontera.
• 3.° Modelar el movimiento del agua a través de la presa. Conocido el campo de velocidades puede determinarse la ley de intercambio de caudales con el exterior y la ley de variación del caudal circulante exteriormente.
• 4.° Comparar esta ley con la estimada en el primer paso.
• 5.° Si coinciden, el problema está resuelto. En caso contrario, debe modificarse la ley y reiterar el proceso. En definitiva, vemos que la condición de frontera correspondiente a la coronación y el paramento de aguas abajo obliga, en pura teoría, a resolver el problema mediante tanteos sucesivos. Desde un punto de vista práctico, debe tenerse en cuenta que en general la altura de lámina vertiente sobre la presa será pequeña, muy inferior a las alturas de presión existentes en el dominio de filtración, por lo que la hipótesis de presión nula a lo largo del talud de aguas abajo resulta en general aceptable, salvo en el caso de que la lámina vertiente sea importante. Esta hipótesis de presión nula simplifica notablemente la modelación, ya que hace innecesaria la iteración.

6. Solución del problema mediante un esquema de diferencias finitas

El método de las diferencias finitas resulta más intuitivo y conduce a una formulación más sencilla que el método de los elementos finitos. En cambio, este último tiende a ser más rápido y permite un mejor tratamiento de contornos de geometría complicada (Volker, 1969; Parkin, 1971)[5]  and [4] . Para el desarrollo de la formulación en diferentes finitas se utilizará la siguiente notación simplificada:

y denominaremos N  = 1/m , quedando la ecuación diferencial de campo de Parkin:

6.1. Formulación general en un nudo interior

Se considera una malla rectangular de lados b_x y b_y en las direcciones ortogonales de los ejes x e y respectivamente. Denominamos nudo interior al que está rodeado por otros 8 nudos interiores al dominio de filtración. ( fig. 5 )

Fig. 5. Nudo interior.

Las derivadas en un nudo interior pueden aproximarse mediante las siguientes expresiones:

Por tanto:

Por otra parte:

y por tanto:

Finalmente:

de donde:

Por tanto, las derivadas de la función potencial ϕ (x , y ) en un nudo interior O (i , j ) pueden aproximarse así:

Para obtener una formulación más sencilla, llamaremos:

De este modo queda:

Por tanto, en el punto O (i , j ), omitiendo los subíndices que en todos los casos son (i , j ):

y operando:

Para simplificar nuevamente la formulación, llamaremos:

Sustituyendo en la ecuación diferencial de campo tenemos:

Denominando:

queda:

Llamando ahora:

se obtiene:

Si llamamos:

tenemos:

de donde:

y finalmente, denominando:

queda la expresión simple:

en la que puede despejarse:

Deshaciendo ahora los cambios de variable realizados, salvo los de K_x y K_y , para no complicar la expresión, se obtiene fácilmente la expresión del potencial ϕ (i , j ) en un nudo interior O (i , j ) en función de las derivadas ϕ_x (i , j ), ϕ_y (i , j ), ϕ_xy (i , j ) y de los valores de K_x (i , j ) y K_y (i , j ), que pueden todos ponerse en función de los valores de ϕ en los nodos adyacentes, como se vio con anterioridad:

siendo:

6.2. Formulación en un nudo de contorno impermeable horizontal (base del espaldón)

En un nudo situado en un contorno impermeable horizontal, como la base del espaldón, no se dispone de nudos situados bajo el mismo, pero la situación puede resolverse imponiendo la condición de que la variación del potencial de velocidad en la dirección vertical tiene que ser nula

Esto debe cumplirse, ya que en caso contrario existiría circulación de agua en sentido vertical, lo que no es posible siendo el contorno horizontal impermeable. Imponiendo esta condición, se llega a la siguiente expresión:

siendo:

6.3. Formulación en un nudo de contorno impermeable vertical

En este caso debe ser:

Procediendo de forma análoga obtenemos la expresión del potencial en un nudo de contorno impermeable vertical:

siendo:

6.4. Formulación en un nudo afectado por el contorno inclinado del talud

Únicamente resulta problemática la aproximación de la derivada cruzada:

De acuerdo con la notación de la figura 6 :

El valor de ϕ_x en el punto C puede aproximarse de diversas formas. Se ha optado por hacer una extrapolación lineal a partir de los valores de ϕ_x con los puntos D y L .

Extrapolando linealmente:

Fig. 6. Nudo afectado por el contorno inclinado del talud.

6.5. Punto singular intersección de dos contornos impermeables, uno horizontal y el otro vertical

En este punto singular el potencial de velocidad se puede aproximar considerando el valor medio de las extrapolaciones lineales realizadas a partir de los dos nudos más próximos en fila y columna, es decir, en el contorno impermeable horizontal y vertical cuya intersección origina el punto singular. Considerando dicho punto singular como origen de filas y columnas, tendrá i = 0 y j = 0. La extrapolación a partir de los dos nodos más próximos de la fila 1, contorno horizontal impermeable, será:

y la extrapolación a partir de los dos nodos más próximos de la columna 1, contorno vertical impermeable:

Por tanto:

Esta formulación en diferencias finitas se ha implementado en el código que se ha denominado Movimiento No Lineal en Espaldón de Escollera. Red de filtración (MNLEE.RDF).

7. Resultados del modelo numérico

Una vez resuelto el problema definido por la ecuación diferencial de campo y las condiciones de contorno (mediante el código MNLEE.RDF), el potencial de velocidad (ϕ ) es conocido en todos los nudos del dominio de filtración. A partir de ellos puede obtenerse también en cada nudo la velocidad de filtración y el gradiente hidráulico, que definen sendos campos vectoriales, y el campo escalar de presiones del agua. También puede obtenerse el caudal filtrado a través de la escollera.

7.1. Presión del agua

El potencial de velocidad en un punto del dominio de filtración es:

siendo «h »la carga hidráulica o carga piezométrica, ya que la carga de velocidad es despreciable frente a la suma de altura y presión.

y c la constante de la fórmula de resistencia:

Por tanto, conocido ϕ (i , j ) en cada nudo del dominio, una vez resuelto el problema, de forma inmediata puede obtenerse la carga h (i , j ) correspondiente.

y la presión del agua

Si bien la presión es fluctuante en las proximidades del talud, debido a la turbulencia, la fluctuación se atenúa rápidamente hacia el interior del medio poroso (Detert et al., 2010)[6] .

7.2. Velocidad de filtración

Como se vio al presentar la ecuación de campo de Parkin, las componentes del vector velocidad de filtración son:

A partir de estas ecuaciones puede obtenerse el vector velocidad de filtración en cada nudo del dominio. Su módulo será:

y el ángulo de avance respecto a la horizontal

con sentido antihorario positivo.

7.3. Gradiente hidráulico

El vector gradiente hidráulico en un punto lleva la dirección de la velocidad de filtración, o sea, es tangente a la línea de filtración que pasa por ese punto, y su módulo puede obtenerse directamente a partir de su propia definición o bien aplicando la fórmula de resistencia, una vez conocida la velocidad de filtración en el punto considerado. Por la propia definición de gradiente hidráulico:

Teniendo en cuenta que:

donde c es una constante, tenemos:

Las coordenadas de ${\textstyle {\boldsymbol {grad}}{\mbox{ }}\phi }$ son:

Por tanto, el módulo del gradiente hidráulico puede obtenerse en el nudo (i , j ):

En cada nudo los valores de las derivadas pueden obtenerse mediante las expresiones deducidas con anterioridad. También puede aplicarse la fórmula de resistencia, una vez calculada la velocidad en cada nudo como se vio anteriormente. El gradiente hidráulico tendrá el siguiente módulo:

Obviamente, el valor obtenido por medio de cualquiera de las dos expresiones debe ser idéntico en todos los nudos, como se ha comprobado mediante el código MNLEE.RDF.

7.4. Caudales filtrados

Admitiendo que el espaldón de aguas abajo está totalmente saturado por el agua que vierte sobre la coronación de la presa, el caudal que circula a través de distintas secciones verticales es distinto, ya que entre cada 2 secciones se produce una aportación o detracción de agua a través del talud (fig. 7 ).

Fig. 7. Intercambio de caudales entre el interior y el exterior del espaldón de escollera.

Denominaremos: Q_c (j ) el caudal que circula a través de la sección definida por la columna j , Q_f (j ) caudal que se filtra a través del talud entre las secciones definidas por las columnas «j »y «j  + 1». El caudal Q_c (j ), y análogamente Q_c (j  + 1), puede obtenerse por suma de los caudales circulantes entre cada dos nudos. Conocida la componente horizontal de la velocidad en cada nudo, puede considerarse una velocidad media en el interespacio entre cada dos nudos consecutivos en la misma vertical. Por tanto: