Simulación de flujo unidimensional en canal abierto (página 2)

El radio hidráulico () se utiliza en la ecuación (3) y en todas las demás ecuaciones; pero, cuando el canal tiene una configuración geométrica arbitraria (canal natural), se puede sustituir el radio hidráulico por la profundidad hidráulica. Esta aproximación () se asume válida para cuerpos de agua poco profundos, y se utiliza por la facilidad de calcular el ancho de la superficie líquida en lugar del perímetro hidráulico.

El coeficiente de Boussinesq , también llamado el coeficiente de momentum, está presente en la ecuación de movimiento para tomar en cuenta las distribuciones de velocidad no uniformes en las secciones transversales.

Las ecuaciones y describen, en general, el flujo no permanente en un canal de sección arbitraria teniendo áreas de transporte y de almacenamiento (o solamente de transporte). En su formulación, se asume que el agua es de densidad homogénea, que la presión hidrostática prevalece en todo el canal, que la pendiente de fondo del canal es pequeña y uniforme, que no hay procesos de transporte de sedimentos en el lecho del canal (no ocurre erosión ni sedimentación), que la geometría del tramo es suficientemente uniforme para permitir la aproximación unidimensional, y que la resistencia por fricción es la misma como en el flujo permanente, permitiendo el uso de la ecuación de Chézy o Manning.

1.3 Esquema Implícito de Diferencias Finitas

Existen numerosos métodos numéricos para producir soluciones aproximadas de las ecuaciones de flujo. En este trabajo, las ecuaciones de flujo serán discretizadas mediante el esquema de diferencias finitas implícito de Preissmann. Esta técnica, permite que el modelo utilice segmentos de diferentes longitudes y un esquema que va desde centrado hasta totalmente adelantado en el tiempo.

El método de solución implícito se emplea debido a su eficiencia inherente y propiedades de estabilidad superior. Es posible agregar un procedimiento de iteración opcional controlable por el usuario para mejorar la exactitud de los resultados.

Fig. 2. Grilla espacio-temporal Esquema Preissmann

El sistema de grilla espacio-temporal de la Fig. 2 muestra la región en que las ecuaciones de flujo son resueltas. Las derivadas temporal y espacial del valor funcional, ƒ, que representa la variable dependiente, nivel (elevación de la superficie líquida) o caudal, son discretizadas de la siguiente manera ((Abbott 1989) en Yzocupe 2006):

(6)

(7)

Donde, , y son factores de ponderación utilizados para especificar la posición temporal y espacial, respectivamente, dentro del incremento de tiempo e incremento de distancia en el cual la derivada y las funciones serán evaluadas.

Tomando , produce una derivada temporal en la posición espacial . Similarmente, cuando la derivada espacial esta centrada en la dirección temporal . Los errores de y , pero tomando se introducen errores de truncamiento que produce disipación numérica. Las derivadas temporales normalmente son calculadas con , aunque otros valores pueden ser ventajosos cuando se utilizan segmentos de longitudes desiguales ((Abbott 1989) en Yzocupe 2006).

De manera similar al tratamiento de la derivada espacial, el área de la sección transversal, el ancho de la superficie libre, el radio hidráulico, y las descargas en forma no derivativa, denotadas por ƒ (x,t), se discretizan como sigue ((Schafframek 1987) en Yzocupe 2006):

(8)

El factor de ponderación puede ser asignado en el rango . Así, estos valores funcionales pueden representarse en cualquier nivel de tiempo como las derivadas espaciales.

La determinación de valores apropiados para estos parámetros es importante porque ellos tienen efecto en la precisión, convergencia, y estabilidad del modelo. Tales valores son la determinación del incremento de tiempo (), la longitud de los segmentos del canal (), y la selección de los factores apropiados de ponderación del esquema de Preissmann.

Aproximación de las Ecuaciones de Flujo

Ecuación de Conservación de Masa

Ecuación de Conservación de Cantidad de Movimiento

Representación en diferencias finitas de las derivadas y coeficientes:

1.4 Discretización de las Ecuaciones Gobernantes

Las ecuaciones diferenciales parciales de flujo (1) y (5) son transformadas en expresiones discretas mediante la aplicación del esquema implícito de diferencias finitas de Preissmann, utilizando los operadores definidos en las ecuaciones (6), (7) y (8). Se utilizo la tilde (~ ) para denotar las cantidades tomadas como constantes locales, las que serán actualizadas a través de las iteraciones en el proceso de cálculo.

Rescribiendo la ecuación de conservación de masa :

(9)

(10)

(11)

(12)

(13)

(14)

Ecuación de Conservación de Cantidad de Movimiento

Rescribiendo:

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22 a)

(22 b)

(23)

(24)

Las ecuaciones algebraicas lineales (13) y (19), que definen el flujo para el segmento , se pueden expresar también en la siguiente forma matricial ((Schafframek 1987) en Yzocupe 1992):

(25)

Entonces, se obtiene un par de ecuaciones algebraicas lineales por cada segmento de la grilla ;

(26)

Para puntos de cálculo se generará un sistema de (2ii-2) ecuaciones algebraicas lineales con 2ii incógnitas. Por lo que se necesita adicionar dos ecuaciones, las que provienen de las condiciones de frontera, para completar el número de ecuaciones necesarias. Luego de lo cual el sistema podrá ser resuelto.

1.5 Obtención y Solución de las Ecuaciones Algebraicas Lineales

Matriz Ejemplo para un Canal de Cuatro Segmentos

Para ilustrar el calculo matricial se aplica el sistema de ecuaciones obtenido en (26) para un canal de 4 segmentos (4 = ndx )y 5 puntos de grilla (ii = 5), se obtiene el siguiente sistema de ecuaciones lineales, , tal como se muestra en la ecuación (27):

Fig. 3. Grilla de cálculo para un canal de 4 segmentos

  (27)

Es decir, se obtiene un sistema matricial con una matriz de coeficientes penta diagonal,, de orden 2ii, un vector de incógnitas, x, de tamaño 2ii, y un vector de residuos, b, de tamaño 2ii. El sistema matricial puede resolverse, con las condiciones iniciales apropiadas, mediante el método de Eliminación de Gauss o también utilizando el método de Doble Barrido. Para nuestro caso se usó el método de Gauss con pivoteo parcial.

Después de resolver el sistema de ecuaciones se proporcionan los valores calculados a las variables del nivel de tiempo (la primera aproximación en el proceso de iteración fue al adoptar los valores del tiempo precedente). Luego, los coeficientes de las ecuaciones lineales pueden ser recalculados y actualizados y el sistema podrá ser resuelto de nuevo. Este proceso necesita de dos a tres iteraciones para obtener una muy buena aproximación, pero ya en la segunda iteración se obtiene una solución satisfactoria.

1.6 Condiciones Iniciales y de Frontera

Condiciones Iniciales

Para iniciar la solución del sistema de ecuaciones algebraicas lineales, se requieren los valores de las variables de flujo para el tiempo cero. Tales valores se pueden obtener de datos medidos o calculados de alguna otra fuente, tales como aproximaciones para un estado permanente, o resultados de alguna simulación anterior.

El uso sucesivo de los valores calculados como las nuevas condiciones iniciales permite que el proceso de cómputo proceda paso a paso hasta concluir la simulación. Una convergencia exitosa del cómputo a la solución correcta requiere que los valores iniciales sean razonablemente precisos; a menor precisión de los valores iniciales, mayor tiempo se tomará para disipar los errores iniciales y llegar a la solución correcta.

Condiciones de Frontera

La solución de las ecuaciones de flujo requiere que se especifiquen condiciones de frontera en los extremos del canal (Fig. 3) durante todo el tiempo de simulación para proveer el número suficiente de ecuaciones adicionales y satisfacer los requerimientos de la técnica de solución.

Estas condiciones de frontera pueden ser elaboradas a partir de registros históricos o calcularse mediante funciones especificadas por el usuario. Se tienen varias combinaciones de condiciones de frontera externas, estas pueden consistir de una descarga cero (por ejemplo, al final del canal), un caudal o nivel conocido en función del tiempo, o una curva de calibración conocida. Las condiciones de fronteras tipo series de tiempo, pueden ser leídas por el programa desde archivos de datos (Yzocupe 2004).

También se ha implementado una condición de frontera no reflectante, la cual consiste de una relación matemática que permite que las perturbaciones u ondas pasen libremente por la frontera y no se reflejen y regresen dentro del dominio de cálculo ((Vreugdenhil 1989) en Yzocupe 2006).

1. Algoritmo general para condiciones de fronteras aguas arriba:

2. (28)

   Si se proporciona

   

   

   Si se proporciona

   

   

3. Algoritmo general para condiciones de fronteras aguas abajo:

(29)

Si se proporciona

Si se proporciona

1.7 Diseño del Algoritmo de Solución

Caso de Flujo Permanente

Caso de Flujo No Permanente

CAPÍTULO 2. PRUEBA DE PERFORMANCE DEL MODELO

En el presente capitulo se realizan dos pruebas de simulación de flujo estático y permanente con el fin de evaluar la calidad y/o perfomance del modelo en un canal rectangular singular Fig. 5. Se han designado cinco puntos de control Fig. 4 para evaluar las variación de caudal y altura durantes el periodo de simulación (24 horas).

Fig. 4 Ubicación de puntos de control en canal 1D (longitud 10 Km.)

Fig. 5 Sección transversal (canal rectangular)

1. El canal tiene una longitud total de 10 Kms y una pendiente de de fondo (Sf) igual a . La sección transversal es rectangular prismático y tiene un ancho de 100 m. El canal se divide en 10 segmentos , y 11 secciones transversales de las cuales se han elegido 5 secciones transversales como puntos de control como se muestra en la Fig. 4, el tiempo de simulación de es de 24 horas, con un paso de tiempo , lo que produce un 865 pasos de tiempo. El coeficiente de rugosidad de manning es constante e igual a 0.026 en todo el sistema, consideramos constante el valor de boussinesq y

2. Condiciones y Parámetros de Control

   Se realizan dos pruebas de flujo estático, en la primera se coloca valores de caudal en la frontera izquierda y derecha, en la segunda prueba se colocan valores de altura en la frontera izquierda y derecha, en ambas pruebas fueron cambiando los factores de ponderación del esquema de Preissmann

   Prueba estática (Q en las fronteras)

   Las condiciones iniciales y de frontera para la prueba estática del flujo en canal unidimensional se resumen en el siguiente cuadro:

   CONDICIONES INICIALES
   Caudal	Altura
   0	10

   CONDICIONES DE FRONTERA
   Frontera Aguas Arriba
   Frontera Aguas Abajo

   Los resultados de la simulación se han registrado durante 24 horas de simulación en 5 estaciones de control ubicadas a 1, 3, 5, 7, 9 kms. A continuación se muestran los hidrogramas de niveles y caudales para los diferentes factores de ponderación de Preissmann ()

   

   Fig. 6 Hidrograma de Niveles

   

   Fig. 7 Hidrograma de Caudales

   

   Fig. 8 Hidrograma de Niveles

   

   Fig. 9 Hidrograma de Caudales

   

   Fig. 10 Hidrograma de Niveles

   

   Fig. 11 Hidrograma de Caudales

   Los hidrogramas muestran que los niveles permanecen constantes durante todo el tiempo de simulación alrededor de 10 m para todos los factores de ponderación de Preissmann, en los hidrogramas de caudales se aprecian oscilaciones de amplitudes variadas en cada uno de los factores de ponderación de Preissmann, las fluctuaciones se presentan en el orden de lo que se puede considerar que el modelo reproduce las condiciones iniciales y de fronteras de , además se observa que las oscilaciones tienen un comportamiento en fase en los cinco puntos de control.

   Prueba estática (H en las fronteras)

   Las condiciones iniciales y de frontera para la prueba estática del flujo en canal unidimensional se resumen en el siguiente cuadro:

   CONDICIONES INICIALES
   Caudal	Altura
   0	10

   CONDICIONES DE FRONTERA
   Frontera Aguas Arriba
   Frontera Aguas Abajo

   Los resultados de la simulación se han registrado durante 3 horas de simulación en 5 estaciones de control ubicadas a 1, 3, 5, 7, 9 kms. A continuación se muestran los hidrogramas de niveles y caudales para los diferentes factores de ponderación de Preissmann ()

   

   Fig. 12 Hidrograma de Niveles

   En el caso de la prueba estática considerando valores de altura en las dos fronteras, tenemos que el modelo produce muy bien los valores de las fronteras durantes las primeras 4 horas de simulación como se aprecia en los hidrogramas de niveles y caudales de los cinco puntos de control en las figuras 12 y 13, Se realizaron pruebas con diversos factores de ponderación de Preissmann, siendo estable en un factores mayores a 0.7 en un periodo de simulación de 4 horas.

   

   Fig. 13 Hidrograma de Caudales

3. Prueba Estática
4. Prueba de estado permanente

Se realizan una pruebas de flujo permanente con valores de caudal en la frontera izquierda y de altura en la frontera derecha, en la pruebas de flujo permanente cambiando los factores de ponderación del esquema de Preissmann

Test de flujo permanente (Q y H)

Las condiciones iniciales y de frontera para la prueba de flujo permanente en canal unidimensional se resumen en el siguiente cuadro:

Los resultados de la simulación se han registrado durante 24 horas de simulación en 5 estaciones de control ubicadas a 1, 3, 5, 7, 9 kms. A continuación se muestran los hidrogramas de niveles y caudales para los diferentes factores de ponderación de Preissmann ()

Fig. 14 Hidrograma de Niveles

Fig. 15 Hidrograma de Caudales

Con factor de ponderación de 0.5 se aprecia en los hidrogramas de niveles y caudales de los cinco puntos de control, tanto los caudales y los niveles presentan oscilaciones alrededor del de los valores de 10 metros en el caso de los niveles teniendo en cuenta las diferencias de nivel que se registran en cada punto de control por la consideración de la pendiente de fondo y el caudal las oscilaciones son alrededor de , presentando oscilaciones de mayor amplitud durante las primeras 14 horas de simulación

Fig. 16 Hidrograma de Niveles

Fig. 17 Hidrograma de Caudales

Con factor de ponderación de 0.6 se aprecia en los hidrogramas de niveles y caudales de los cinco puntos de control, tanto los caudales y los niveles presentan oscilaciones alrededor del de los valores de 10 metros en el caso de los niveles se muestra claramente en cada punto de control el efecto de la pendiente de fondo y el caudal las oscilaciones son alrededor de , presentando oscilaciones de mayor amplitud durante las primeras 10 horas de simulación

Fig. 18 Hidrograma de Niveles

Fig. 19 Hidrograma de Caudales

Con factor de ponderación de 0.7 se aprecia en los hidrogramas de niveles y caudales de los cinco puntos de control, tanto los caudales y los niveles presentan oscilaciones alrededor del de los valores de 10 metros en el caso de los niveles se muestra claramente en cada punto de control el efecto de la pendiente de fondo y el caudal las oscilaciones son alrededor de , presentando oscilaciones de mayor amplitud durante las primeras 8 horas de simulación

Fig. 20 Hidrograma de Niveles

Fig. 21 Hidrograma de Caudales

Con factor de ponderación de 0.8 se aprecia en los hidrogramas de niveles y caudales de los cinco puntos de control, tanto los caudales y los niveles presentan oscilaciones alrededor de los valores de 10 metros en el caso de los niveles se muestra claramente en cada punto de control el efecto de la pendiente de fondo y el caudal las oscilaciones son alrededor de , presentando oscilaciones de mayor amplitud durante las primeras 7 horas de simulación

Fig. 22 Hidrograma de Niveles

Fig. 23 Hidrograma de Caudales

Con factor de ponderación de 0.9 se aprecia en los hidrogramas de niveles y caudales de los cinco puntos de control, tanto los caudales y los niveles presentan oscilaciones alrededor del de los valores de 10 metros en el caso de los niveles se muestra claramente en cada punto de control el efecto de la pendiente de fondo y el caudal las oscilaciones son alrededor de , presentando oscilaciones de mayor amplitud durante las primeras 5 horas de simulación

CAPÍTULO 3. DISCUSIÓN DE RESULTADOS

El modelo de flujo estático y flujo permanente unidimensional en canales regulares abiertos muestra buenos resultados en las diferentes pruebas realizadas bajo ciertas condiciones iniciales y de frontera

Los precisión de los resultados son del orden de con diferentes valores de ponderación en el esquema numérico, en el caso de la prueba estática bajo condiciones de frontera de caudal cero se obtiene un flujo que permanece estable durante todo el periodo de simulación con valores aproximados a cero,

En las pruebas de sensibilidad en flujo permanentes del factor de ponderación se logro determinar que mediante el aumento del factor de ponderación se logra un menor tiempo de estabilidad en las curvas de evolución de los caudales y niveles del canal.

La técnica de diferencias finitas implícitas de Preissmann con factores de ponderación da una mayor flexibilidad en el manejo de las condiciones de fronteras en los canales

Cabe resaltar que el modelo es hidrodinámico y no contempla transporte de sedimentos ni cambios morfológicos en los tramos del canal, pero si considera una variación de pendiente de fondo constante.

REFERENCIAS

Abbott M. B.; Basco D. R. (1989): Computational Fluid Dynamics. An Introduction for Engineers. Logman Scientific, 425pp.

Becerra A.; Hernán H. (2006): Programa Computación de Simulación Hidráulica del Riego por Surco Usando el Modelo de Onda Cinemática, Dyna, Año 73, Nro. 149 pp 107-117.

Yzocupe, V. (1992): Simulación del Desplazamiento de Onda de Crecida (Flood Routing). Tesis de Ingeniero en Mecánica de Fluidos.

Yzocupe, V. (2006): Simulación de Flujo 1D en Canales Abiertos, Revista de Investigación de Física. Vol. 9 N° 1, pp. 1-11.

David Correa 1,

Juan Huamaní 1,

José Tezén 1,

José Mesías 1

ricardo_mesia[arroba]hotmail.com

1 Maestría en Ingeniería de Mecánica de Fluido Computacional, Universidad Nacional Mayor de San Marcos (Lima – Perú)