( ƒj n + 1 + ƒj + 1n + 1 )                     ( ƒj n + ƒj + 1n ) ƒ (x,t )  ≈  θ ______________________   +  (1 - θ) ________________                                      2                                                2

(7)

∂                         ( ƒj + 1n + 1 - ƒj n + 1 )                       ( ƒj + 1n - ƒj n ) ____ ƒ (x,t )  ≈  θ ______________________   +  (1 - θ) ________________  ∂x                                   ∆ x                                               ∆ x

(8)

∂                        1                         _____ ƒ (x,t )  ≈  _____ [ ( ƒj n + 1 + ƒj + 1n + 1 ) - ( ƒj n + ƒj + 1n ) ]  ∂ t                      2∆t

(9)

en la cual ƒ = variable dependiente, y θ = factor de ponderación del esquema implícito.

La aplicación de las Ecs. 7-9 a las Ecs. 3 y 4 en forma discreta en el plano x - t produce un conjunto de ecuaciones algebraicas no lineales. Para resolver los valores de las variables dependientes en el nivel avanzado de tiempo sin recurrir a un proceso iterativo, se diseña un esquema de linearización. Este esquema permite la conversión de un sistema de ecuaciones no lineales en un sistema de ecuaciones lineales que se puede resolver directamente, sin iteración. El esquema de linearización genera errores adicionales a los ya impuestos por la solución discreta. Sin embargo, el tamaño y el crecimiento de estos errores se pueden controlar de manera efectiva permaneciendo dentro de los límites del método computacional, es decir, evitando cambios bruscos en los valores de las variables dependientes. En términos de efectividad, el esquema de solución linearizada se compara favorablemente con respecto a la solución no lineal.

El esquema de linearización tiene la siguiente forma:

∆ƒj  ( ƒj n + 1 )δ ≈ ( ƒj n )δ + δ( ƒj n )δ-1∆ƒj  (válido para _____ < < 1 )                                                                                    ƒj n

(10)

en la cual δ = número real.

La solución del sistema de ecuaciones lineales simultáneas se logra mediante el conocido algoritmo de doble barrida (6).

4.  CONDICIONES DE CONTORNO

En la frontera aguas arriba, se puede especificar el cambio en el nivel o la elevación del lecho con el tiempo. Alternativamente, se puede especificar la profundidad del flujo. Esto permite el modelado de un hidrograma de flujo de entrada de sedimentos, al relacionar éste con el cambio en la profundidad del flujo. Para diferentes descargas de agua y sedimentos, la condición de frontera aguas arriba se formula como sigue (7):

qs1 n                 q1 n + 1  ∆d1T = ∆d1S + ∆d1Q = d1 [ (__________ )1/m(__________ ) - 1 ]                                                   qs1 n + 1                q1 n

(11)

en las cuales ∆d_1T = cambio en la profundidad de flujo en la condición de frontera aguas arriba; ∆d_1S = cambio en la profundidad de flujo en la condición de frontera aguas arriba debido a la entrada de sedimentos en desequilibrio; y ∆d_1Q = cambio en la profundidad de flujo en la condición de frontera aguas arriba debido a la entrada de agua en desequilibrio. Para la descarga de agua constante en el tiempo, la Ec. 11 se reduce a lo siguiente:

qs1 n                 ∆d1S = d1 [ (_________ )1/m - 1 ]                         qs1 n + 1

(12)

Las Ecuaciones 11 y 12 se basan en la ecuación de transporte de material del lecho (Ec. 6) y son estrictamente de tipo no lineal. Una comparación entre esta formulación y otra basada en una linealización en la frontera [ver Cunge y Perdreau (3)] muestra que la primera es intrínsecamente precisa, mientras que la última puede resultar en errores graves en la frontera.

La condición de frontera aguas abajo se da generalmente como una relación caudal-tirante. Para el caso de caudal constante, el tirante aguas abajo generalmente se mantiene constante.

5.  ESTABILIDAD Y CONVERGENCIA

Todo modelo numérico debe satisfacer ciertos requisitos de estabilidad y convergencia. La estabilidad se refiere a la capacidad del esquema numérico para inhibir la generación y crecimiento de errores. La convergencia es una medida de la magnitud de los errores en la solución numérica, debido a una inadecuada discretización. En los esquemas implícitos, la estabilidad se relaciona con los errores de redondeo (9); la convergencia, con los errores de discretización. La estabilidad se ve afectada si la discretización se hace más fina, mientras que lo contrario se aplica desde el punto de vista de la convergencia. Ambas propiedades están en conflicto en esquemas numéricos: cuanto más estable es un esquema, menos convergente, y viceversa.

El primer autor, et al. (10) han presentado un análisis exhaustivo de las propiedades numéricas de los modelos transitorios de lecho implícito. Según este estudio, el comportamiento del modelo está controlado por tres parámetros numéricos y dos físicos. Los parámetros numéricos son: (1) La resolución espacial L/∆x; (2) la resolución temporal T/∆t; y (3) el factor de ponderación. Las resoluciones espaciales y temporales se refieren a la discretización del modelo en el tiempo y el espacio, en la cual L y T son la longitud característica y la duración de la onda en estudio, respectivamente. La relación entre la resolución espacial y la temporal es un valor adimensional denominado número de Courant (1). Para un L, T, y ∆x dados, el número de Courant es una medida del intervalo de tiempo ∆t.

Los parámetros físicos son: (1) Número de Froude del flujo de equilibrio F_o, definido como F_o = u_o / √ gd_o ; y (2) número de onda adimensional σ_*, definido como σ_* = (2π / L)(d_o / S_o), en las cuales u_o y d_o = velocidad y profundidad del flujo de equilibrio, respectivamente, y S_o = pendiente de fondo.

6.  EXPERIMENTOS NUMÉRICOS

Los experimentos numéricos tienen como objetivo evaluar el comportamiento general del modelo, medido por su capacidad para simular problemas transitorios hipotéticos, los cuales producen soluciones físicamente reales. Además, la experimentación numérica se utiliza para verificar el análisis teórico de estabilidad y convergencia.

Formación de ondas de arena. Las ondas de arena se definen como acumulaciones transitorias de sedimentos en el lecho de un río. Se sabe que varios fenómenos como inundaciones, deslizamientos de tierra y fallas de presas provocan la formación de ondas de arena. Una vez que se ha formado la onda, se traslada lentamente aguas abajo en condiciones de flujo subcrítico.

El modelo matemático presentado en este artículo simula efectivamente la formación de ondas de arena. Un hidrograma de carga de material del lecho que no está en equilibrio (en exceso o defecto de la velocidad de transporte de equilibrio) crea una onda de arena positiva (agradación) o negativa (degradación) a medida que ingresa al canal. El proceso de formación de ondas de arena está sujeto a una restricción numérica: el número de Courant C del lecho móvil debe ser igual o mayor que uno. De lo contrario, se presentan oscilaciones parasitarias en el modelo y, finalmente, la solución pierde sentido. Por lo tanto, al modelar la formación de ondas de arena, es de suma importancia tener un conocimiento previo de la rapidez de la perturbación del lecho (4). Con la celeridad conocida al menos aproximadamente, se verifica el número de Courant del lecho para asegurarse de que se cumpla la restricción numérica: C ≥ 1.

Migración de ondas de arena. Las ondas de arena migran aguas abajo en un flujo subcrítico y están sujetas a disipación. La cantidad de disipación observada en el modelo numérico suele ser la combinación de dispersión física y numérica. La cantidad de dispersión numérica es una indicación de la convergencia de la solución; la solución será más convergente a medida que se minimice la dispersión numérica. En general, esto se puede lograr utilizando una resolución espacial L/∆x tan alta como sea posible, un número Courant próximo a 1, y un valor del factor de ponderación en el rango 0.5 ≤ θ ≤ 0.75. Sin embargo, la estabilidad puede verse afectada si se busca la eliminación de toda la dispersión numérica; por lo tanto, en la práctica, puede ser necesario mantener una pequeña cantidad de dispersión numérica.

Ondas de lecho de erosión y deposición. El modelo matemático se puede utilizar para estudiar varios fenómenos transitorios de erosión y de deposición. La Figura 2 muestra el caso de una presa de retención de sedimentos. Este fenómeno transitorio se simula incrementando el tirante aguas abajo para reflejar la presencia de la presa. Esto provoca un aumento en el tirante en las secciones aguas arriba, reduciendo las velocidades de flujo y produciendo deposición. Una onda de arena se deposita en el extremo aguas arriba del tramo y migra aguas abajo lentamente, llenando eventualmente todo el volumen de retención.

Fig. 2  Deposición de sedimentos aguas arriba de una presa de retención.

Otros ejemplos de formación y propagación de ondas, como el proceso de degradación debajo de una presa (Fig. 3) y el comportamiento de cortes de dragado y hoyos de préstamo (Fig. 4), se pueden estudiar de manera similar utilizando el modelo matemático aquí desarrollado.

Fig. 3  Degradación debajo de una presa.

Fig. 4  Migración de un hoyo de extracción de arena.

Modelado de caudal variable. Para modelar un flujo de entrada de caudal variable en la frontera aguas arriba, el hidrograma se discretiza en varios intervalos. Se hace que cada intervalo del hidrograma corresponda a uno o más intervalos de la onda del lecho. La simulación del caudal líquido se logra mediante un cálculo de remanso entre los intervalos de tiempo de la onda del lecho.

Sensibilidad a los parámetros numéricos. Se diseñó una serie de simulaciones de prueba con el propósito de evaluar la sensibilidad del modelo a los parámetros numéricos. El problema es el de rastrear la traslación y atenuación de una onda de arena sinusoidal a medida que viaja aguas abajo en un canal de ancho unitario. Los parámetros físicos se establecieron como sigue: F_o = 0.167, y σ_* = 132. Se especificó como condición inicial una onda de arena sinusoidal de 1 milla (1.6 km) de longitud y una amplitud de 1 pie (0.305 m). Los parámetros numéricos se variaron dentro de los siguientes rangos:

0.5 ≤ θ ≤ 1.0

(13)

L            10 ≤ _____  ≤ 100           ∆x

(14)

1 ≤ C ≤ 8

(15)

Las Figuras 5 (a) a 5 (c) muestran los resultados de la prueba numérica. La Figura 5 (a) muestra que cuanto mayor es el valor de θ, más estable es el modelo, es decir, mayor es la dispersión numérica. Un valor de θ = 0.6 muestra un patrón débil de inestabilidad. La Figura 5 (b) muestra que los valores bajos de L/∆x introducen la dispersión numérica de la misma manera que los valores altos de θ. La Fig. 5 (c) muestra que los valores altos de C también introducen dispersión numérica de una manera similar a un alto valor de θ. En general, la estabilidad mejora aumentando θ, incrementando C, y reduciendo L/∆x; lo contrario es cierto desde el punto de vista de convergencia. Por lo tanto, es necesaria una compensación entre las demandas conflictivas de estabilidad y convergencia. Un valor de compensación apropiado para θ parece ser θ = 0.75 (8). Los valores de L/∆x ≥ 40 son necesarios para minimizar la dispersión numérica introducida al tomar θ = 0.75. El valor apropiado del número de Courant C es una función de la resolución espacial L/∆x y la resolución temporal T/∆t. El segundo autor (7) ha demostrado que los valores de T/∆t > 20 son necesarios para permanecer dentro de límites razonables de error. Por lo tanto, para un cierto T/∆t, cuanto mayor sea el L/∆x , mayor será el número de Courant que se puede utilizar con éxito en el modelo.

Fig. 5   Efecto de: (a) factor de ponderación θ en la solución numérica; (b) resolución espacial L/∆x en solución numérica; y (c) número de Courant C del lecho en la solución numérica (los perfiles del lecho de arena se muestran a intervalos de 96 días).

Recomendaciones sobre estabilidad y convergencia. Valores en el rango 0.5 ≤ θ ≤ 0.6 pueden resultar en una solución débilmente estable o inestable. Como regla general, la estabilidad requiere que θ > 0.6. Se recomienda un valor de θ = 0.75, siempre que se utilice junto con una resolución espacial suficientemente alta (L/∆x ≥ 40) y un número de Courant, mayor que 1 lo más bajo posible. Un valor de C ≥ 1 es necesario para evitar oscilaciones debido a una mala posición de la condición de frontera aguas arriba.

El segundo autor (7) ha derivado un criterio de convergencia aproximado basado en simulaciones numéricas. Definió un parámetro de convergencia C_p como sigue:

L                         ____                 ∆x   Cp  =  _______                       C θ

(16)

Valores altos de C_p resultan en una buena convergencia; por el contrario, valores bajos de C_p son un indicativo de poca convergencia. Un límite inferior de C_p puede estar dictado por consideraciones de precisión, mientras que un límite superior puede ser necesario desde el punto de vista económico. La Ecuación 16 también se puede expresar de la siguiente manera:

T                          ____                  ∆t   Cp  =  ________                          θ

(17)

La Figura 6 muestra los resultados del análisis de convergencia. El grado de convergencia se define como la relación entre la celeridad numérica y la analítica, para soluciones numéricas en las cuales se varían T/∆t y θ. La Figura 6 resalta la importancia de la resolución temporal T /∆t como el parámetro de convergencia primario, a la vez que muestra el rol secundario del factor de ponderación θ. El rango recomendado para la resolución temporal es el siguiente: 20 ≤ T/∆t ≤ 100. Los valores de rango medio representarán una compensación entre precisión y limitaciones de costo.

Fig. 6  Grado de convergencia en función de la resolución temporal T/∆t y el factor de ponderación θ.

7.  RESUMEN Y CONCLUSIONES

Se desarrolla un modelo matemático para simular la formación y propagación de ondas transitorias del lecho en canales. El modelo es del tipo de caudal conocido, es decir, dentro de un intervalo de tiempo, el caudal líquido se considera constante. Esta técnica permite el modelado numérico implícito de la propagación de ondas transitorias del lecho, utilizando intervalos de tiempo mayores de lo que sería posible utilizando el método de enrutamiento secuencial convencional.

La inestabilidad debida a la mala posición de la función de transporte de sedimentos se evita utilizando una función potencial entre la descarga del material del lecho y la velocidad media del flujo. Los errores en la conservación de la masa debido a una linealización en la frontera se eliminan mediante el uso de una formulación no lineal de la condición de frontera aguas arriba. El modelo puede tener en cuenta la descarga de agua que varía en el tiempo mediante la interfaz de un cálculo de remanso entre uno o varios intervalos de tiempo. Esta característica permite una gran flexibilidad para modelar la superficie del agua a largo plazo y las ondas transitorias del nivel del lecho.

La estabilidad y la convergencia son analizadas, estableciéndose criterios para asegurar la mejor operación del modelo. En particular, se muestra que la convergencia es principalmente una función de la resolución temporal T/∆t, con el factor de ponderación θ desempeñando un papel secundario.

APÉNDICE I.  BIBLIOGRAFÍA

1. Abbott, M. A., A. Damsgaard, y G. S. Rodenbuis. 1973. "System 21, 'Jupiter,' A Design System for Two-Dimensional Nearly-Horisontal Flows," Journal of Hydraulle Research, Vol. 11, No. 1, 1-28.

2. Chen, Y. H. 1973. "Mathematical Modeling of Water and Sediment Routing in Natural Channels," thesis presented to Colorado State University, at Fort Collins, Colo., in partial fulfillment of the requirements for the degree of Doctor of Philosophy.

3. Cunge, J. A., y N. Perdreau. 1973. "Mobile Bed Fluvial Mathematical Models," La Beadle Blanche, Grenoble, France, Vol. 28, No. 7, 561-580.

4. De Vries, M., "River Bed Variations-Aggradation and Degradation," Publication No. 107, Delft Hydraulics Laboratory, Delft, The Netherlands, Apr., 1973.

5. "HEC-6, Scour and Deposition in River and Reservoirs", Report 723-G2-L2470 , U.S. Army Corps of Engineers, Hydrologic Engineering Center, Davis, Calif., Mar., 1976.

6. Liggett, J. A., y J. A. Cunge.1975. "Numerical Methods of Solution of the Unsteady Flow Equations," Unsteady Flow In Open Channels, K. Mahmood and V. Yevjevich, eds., Water Resources Publications, Fort Collins, Colo., 89-182.

7. Lopez-Garcia, J. 1977. "Mathematical Modeling of Alluvial Bed Transients," thesis presented to Colorado State University, at Fort Collins, Colo., in partial fulfillment of the requirements for the degree of Master of Science.

8. Mahmood, K., y V. M. Ponce. 1976. "Mathematical Modeling of Sedimentation Transients in Sand-Bed Channels," Report CER75-76-KM-VMP28, Engineering Research Crater, Colorado State University, Fort collins. Colo., Apr.

9. O'Brien, G. G., M. A. Hyman y S. Kaplan. 1951. "A Study of the Numerical Solution of Partial Differential Equations," Journal of Mathematics and Physics, Vol. 29, No. 4, 1951, 223-251.

10. Ponce, V. M., H. Indlekofer y D. B. Simons. 1979. "The Convergence of Implicit Bed Transient Models," Journal of the Hydraulics Division, ASCE, Vol. 105, Nº HY4, Apr.