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.

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.

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.

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:

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.

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):

  (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.

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.

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:

(28)

Si se proporciona





Si se proporciona





2. Algoritmo general para condiciones de fronteras aguas abajo:

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).

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

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 ()

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

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

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

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

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

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.

David Correa ¹,