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