1. INTRODUCCIÓN
La solución numérica de las ecuaciones de flujo no permanente en canales (ecuaciones de Saint Venant) es un tema establecido en la literatura de ingeniería hidráulica. Generalmente se prefieren los esquemas numéricos implícitos a los esquemas explícitos (9), porque permiten pasos de tiempo más grandes y, por lo tanto, reducen el esfuerzo computacional. Entre los esquemas implícitos, el esquema de cuatro puntos, también conocido como esquema de Preissmann (2, 5), ha recibido una amplia atención. Sin embargo, las propiedades de estabilidad y convergencia de este esquema aún no están claramente definidas.
Las ecuaciones de Saint Venant constituyen un conjunto de ecuaciones diferenciales parciales cuasilineales. En la actualidad, no existe un método para analizar la estabilidad y convergencia de soluciones numéricas de ecuaciones de este tipo. Por otro lado, existen varios métodos para analizar la estabilidad y convergencia de ecuaciones diferenciales parciales lineales (11), entre ellos, en particular, el enfoque heurístico de Von Neumann (7). Tradicionalmente, una salida a la dificultad matemática impuesta por la naturaleza cuasilineal de las ecuaciones de Saint Venant ha sido linealizarlas, asumiendo que el comportamiento del sistema lineal se aproxima al comportamiento más complejo del sistema cuasilineal. En la práctica, el análisis lineal sirve, si no como cuantitativo, al menos como una descripción cualitativa de los fenómenos no lineales.
2. ESTABILIDAD Y CONVERGENCIA
La característica esencial de un modelo numérico es la sustitución de una derivada tal como ∂ƒ/∂s por una relación de diferencias finitas tal como Δƒ/Δs. Al hacerlo, el modelo numérico debe satisfacer ciertos requisitos de estabilidad y convergencia. La estabilidad se refiere a la capacidad del esquema numérico de marchar en el tiempo sin generar una amplificación del error.
En la práctica, la estabilidad es una condición necesaria para el funcionamiento del modelo, ya que un modelo inestable tiene poca o ninguna utilidad. En los modelos implícitos, la estabilidad generalmente se logra a expensas de la convergencia. Por lo tanto, la convergencia es el tema central: Una vez que se demuestre que el modelo es estable, es necesario evaluar sus propiedades de convergencia.
Las propiedades numéricas de los esquemas implícitos de cuatro puntos han sido estudiadas por Liggett y Cunge (5), y Fread (3). Liggett y Cunge utilizaron un sistema simple de ecuaciones diferenciales parciales lineales para derivar conclusiones sobre las propiedades numéricas del esquema implícito de cuatro puntos. Si bien se dieron cuenta plenamente de las limitaciones de su análisis, señalaron, no obstante, que se esperaría que las ecuaciones de Saint Venant se comportaran de manera similar. Los cálculos reales han demostrado, sin embargo, que el valor "teóricamente más preciso" del factor de ponderación, θ = 0.5, no es práctico desde el punto de vista de la estabilidad. En la práctica puede ser necesario un valor de θ ≥ 0.6 (2).
Fread (3) ha analizado las propiedades numéricas del esquema implícito de cuatro puntos.
El presente trabajo tiene como objetivo determinar las propiedades de convergencia de los modelos numéricos implícitos de cuatro puntos del flujo no permanente en canales. Usando el método de Von Neumann, se comparan la celeridad y atenuación de la solución numérica con las de la solución analítica (8,7,11). Las conclusiones se refieren a la evaluación de la convergencia en función de parámetros físicos y numéricos.
3. ECUACIONES DE GOBIERNO
El flujo no permanente unidimensional en canales, por unidad de ancho, se puede expresar mediante un conjunto de las siguientes ecuaciones diferenciales parciales (4):
Ecuación de continuidad:
Ecuación de movimiento:
en las cuales u = velocidad media; d = profundidad de flujo; g = aceleración de la gravedad;
Las Ecuaciones 1 y 2 deben satisfacer el flujo no perturbado (equilibrio o flujo uniforme permanente) para el cual
4. SOLUCIÓN ANALÍTICA
Ponce y Simons han demostrado (8) que la solución analítica de las Ecs. 3 y 4 se caracteriza por la celeridad de propagación y la atenuación de las perturbaciones sinusoidales. Para la celeridad, eligieron la celeridad adimensional c* definida como sigue:
en la cual L = longitud de onda de la perturbación sinusoidal; y T = período de la onda. Para la atenuación, eligieron el decremento logarítmico δ definido de la siguiente forma:
en la cual ato and ato+T son las amplitudes de onda al principio y al final de un período T.
Ambos c* y δ son funciones del número de Froude del flujo de equilibrio Fo y el número de onda adimensional σ*, en el cual Fo = uo / ( gdo )1/2 y σ* = (2π/L) do /So (Apéndice I).
5. SOLUCIÓN NUMÉRICA
El esquema implícito de cuatro puntos se basa en las siguientes fórmulas (5), Fig.1:
en las cuales ƒ(x,t) es una variable dependiente; Δx y Δt son los incrementos de espacio y tiempo, respectivamente; y θ es el factor de ponderación del esquema implícito. El factor de ponderación θ se ingresa en la función y su derivada espacial para permitir una mayor flexibilidad en el funcionamiento actual del modelo. Mientras que el rango estable es 0.5 ≤ θ ≤ 1.0, un rango fuertemente estable es 0.6 ≤ θ ≤ 1.0. La sustitución de las Ecuaciones 7 y 9 en las Ecs. 3 y 4 conduce a las ecuaciones discretizadas en u' y d'. La solución se postula en la siguiente forma exponencial:
en las cuales x* = xSo /do; t* = tuoSo /do y β*R = 2πdo /TuoSo, y u* y d* son funciones de amplitud adimensionales. La sustitución de las Ecs. 10 y 11 en las ecuaciones discretizadas en u' y d' conduce a lo siguiente:
en las cuales i = √-1; ξ = [exp (-iβ* Δt*) - 1 ]; α = σ*Δx* /2; Δt* = Δt uo/ (do /So); y Δx* = Δx / (do /So).
Las Ecuaciones 12 y 13 constituyen un conjunto de dos ecuaciones homogéneas en u* y d*.
en la cual C /c* = Δt* /Δx*. El valor de C se conoce como el número de Courant, y se define como la relación entre la celeridad física c, y la "celeridad de la malla" Δx/Δt. La Ecuación 14 es una ecuación cuadrática compleja en ξ. Puede simplificarse a lo siguiente:
en la cual:
y
Llamando:
y resolviendo la ecuación cuadrática (Ec. 15), las dos raíces son:
en la cual:
y
Dado que ξ = [exp ( -i β* Δt* ) - 1], resulta que:
y
La celeridad adimensional de la solución numérica c*n se define como sigue (8):
lo cual se reduce a:
El decremento logarítmico de la solución numérica δn se define como sigue (8):
Dado que c* es una función de Fo y σ*, y M y N son funciones de Fo, σ*, α, C, y θ, resulta que c*n y δn son funciones de Fo, σ*, α, C y θ. Fo y σ* son los parámetros físicos, siendo Fo = número de Froude y σ* = número de onda adimensional σ* = (2π/L) do /So. Los valores de α, C, y θ son los parámeros numéricos, en los cuales α es la resolución espacial α = π /(L /Δx); C es el número de Courant C = c / (Δx / Δt); y θ es el factor de ponderación del esquema implícito.
6. ANÁLISIS DE CONVERGENCIA
El análisis de convergencia se basa en las siguientes relaciones:
El valor R1 es la relación de atenuación, y R2 es la relación de translación. Para R1 > 1, la atenuación numérica es menor que la atenuación física; para R1 < 1, la atenuación numérica es mayor que la atenuación física. Para R2 > 1, la translación numérica es mayor que la translación física, y para R2 < 1, la translación numérica es menor que la translación física.
La sensibilidad de las relaciones de convergencia a los parámetros físicos y numéricos se ha estudiado variando los parámetros dentro de un rango específico. Para el número de Froude Fo,
El análisis se ha realizado sólo para la onda primaria (8), tomando el signo menos asociado con D y E en las Ecs. 30 y 31. Las limitaciones de espacio impiden un análisis completo de los componentes primarios y secundarios de la propagación de la onda.
La Figura 2(a)-2(c) describe la relación de atenuación R1 como una función de Fo, σ*, α, C, y θ.
Las siguientes conclusiones se deducen de la Fig. 2:
El análisis de la relación de translación R2 muestra que esta relación se aproxima a la unidad para una amplia gama de parámetros Fo, σ*, α, y C. Como era de esperar, para números de onda adimensionales de rango medio (σ* = 10) y baja resolución espacial, las relaciones R2 se alejan de la unidad. Como ilustración de esta tendencia, la Fig.3 muestra los valores de R2 para diferentes Fo, σ*, C, y θ, y α = π / 10.
7. RESUMEN Y CONCLUSIONES
Se presenta un análisis teórico de la convergencia del modelo numérico implícito de cuatro puntos de ondas de aguas poco profundas. Se derivan la celeridad de propagación y el factor de atenuación de la solución numérica, y la convergencia se evalúa estableciendo las relaciones de atenuación y translación dada por las soluciones numéricas y analíticas. Se muestra que la convergencia es una función de la interacción compleja de cinco parámetros: (1) El número de Froude del flujo de equilibrio Fo; (2) el número de onda adimensional σ* = (2π/L) (do/So);
El análisis de convergencia se realiza para números de Froude en régimen subcrítico. Además, debido a las limitaciones de espacio, solo se analiza la convergencia de la onda primaria.
Se presentan las siguientes conclusiones:
La simulación es razonablemente apropiada para pequeñas ondas σ* (cinemáticas/difusivas) y para grandes ondas σ* (inercia-presión), siempre que se tenga el cuidado suficiente de minimizar la cantidad de atenuación numérica.
Para ondas intermedias σ* (ondas dinámicas, en las que predominan tanto la fricción como la inercia), la precisión de la simulación depende en gran medida del valor correcto del factor de ponderación θ. En la práctica, puede ser difícil determinar un valor óptimo de θ para asegurar tanto la estabilidad como la convergencia.
Los valores de θ < 0.5 siempre causan amplificación numérica; los valores de 0.5 ≤ θ < 1 puede causar amplificación o atenuación numérica; y θ = 1 siempre causa atenuación numérica.
8. LIMITACIONES Y APLICACIONES
Los resultados de este artículo se han obtenido utilizando una versión linealizada de las ecuaciones de Saint Venant. Por lo tanto, se hace énfasis en la naturaleza cualitativa de los resultados. En la actualidad, la bondad de un modelo no lineal sólo puede evaluarse mediante experimentos numéricos cuidadosamente diseñados. El análisis lineal presentado en este artículo provée un marco apropiado para el diseño de estos experimentos.
AGRADECIMIENTOS
El tema de investigación de este artículo fue apoyado en parte por la Estación Experimental de la Universidad Estatal de Colorado. Horst Indlekofer recibió una beca de la OTAN para hacer posible su visita a la Universidad Estatal de Colorado.
APÉNDICE Ι. ECUACIONES PARA LA CELERIDAD DE PROPAGACIÓN c* Y DECREMENTO LOGARITMICO δ DE PERTURBACIONES SINUSOIDALES EN EL FLUJO EN CANALES (8)
Estas ecuaciones son: c* = 1 ± D; δ = -2π (B ∓ E)/(|1 ± D |), en la cual A = (1/Fo2) - B2; B = 1/(σ* Fo2); C = (A2 + B2)1/2; D = [(C + A)/2]1/2; E = [(C - A)/2]1/2; Fo = uo / (gdo)1/2; y σ* = (2π / L)(do / So).
APÉNDICE ΙΙ. BIBLIOGRAFÍA
Chaudhry, Y. M., y D. N. Contractor. 1973. "Application of the Implicit Method to Surges in Open Channels," Water Resources Research, Vol. 9, No. 6, Dec., pp. 1605-1612.
Cunge, J. A., y M. Wegner. 1694. "Numerical Integration of the Saint Venant Equations by an Implicit Finite Difference Scheme," La Houille Blanche, Grenoble, France, Vol. 19, No. 1, Jan., pp. 33-39.
Fread, D. L. 1974. "Numerical Properties of Implicit Four-Point Finite Dfference Equations of Unsteady Flow," NOAA Technical Memorandum NWS HYDRO-18, National Oceanic and Atmospheric Administration, Mar.
Liggett, J. A. 1975. "Basic Equations of Unsteady Flow," Unsteady Flow in Open Channels, K. Mahmood and V. Yevjevich, eds., VoL. I, Water Resources Publications, Fort Collins, Colo.
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., Vol. I, Water Resources Publications, Fort Collins, Colo.
Lighthill, M. J., y G. B. Whitham. 1955. "On Kinematic Waves I, Flood Movements in Long Rivers," Proceedings, Royal Society of London, Vol. A 229, May, pp. 281-316.
O'Brien, G. G., M. A. Hyman, y S. Kaplan. 1950. "A Study of the Numerical Solution of Partial Differential Equations," Journal of Mathematics and Physics, Vol. 29, No. 4, pp. 223-251.
Ponce, V. M., y D. B. Simons. 1977. "Shallow Wave Propagation in Open Channel Flow," Journal of the Hydraulics Division, ASCE, Vol. 103, No. HY12, Proc. Paper 13392, Dec., pp. 1461-1476.
Price, R. K. 1974. "Comparison of Four Numerical Methods of Flood Routing," Journal of the Hydraulics Division, ASCE. Vol. 100, No. HY7, Proc. Paper 10659, July, 1974, pp. 879-899.
Quinn, F. H., y E. B. Wylie. 1972. "Transient Analysis of the Detroit River by the Implicit Method," Water Resources Research, Vol. 8, No. 6, Dec., pp. 1461-1469.
Roache, P. J. 1972. Computational Fluid Dynamics, Hermosa Publishers, Albuquerque. N. M.
APÉNDICE ΙΙΙ. SIMBOLOGÍA
Los siguiente símbolos son usados en este artículo:
A, B, C, D, E = parámetros, Eqs. 24-28;
a = amplitud de onda;
C = número de Courant;
c = celereridad de onda;
d = profundidad de flujo;
Fo = número de Froude;
ƒ = función;
g = aceleración de la gravedad;
i = √-1;
j = índice de discretización espacial;
L = longitud de onda;
M, Mo, M1, M2 = parámetros;
N, No, N1, N2 = parámetros;
n = índice de discretización temporal;
P, Q = parámetros;
R1 = relación de convergencia de atenuación;
R2 = relación de convergencia de translación;
Sƒ = pendiente de fricción;
So = pendiente de fondo;
T = período de onda;
t = variable tiempo;
to = tiempo inicial;
u = velocidad media;
x = variable de espacio;
α = parámetro de resolución espacial;
β* = factor de propagación complejo;
β*i = parte imaginaria de β*;
β*R = parte real de β*;
Δx = intervalo de espacio;
Δt = intervalo de tiempo;
δ = decremento logarítmico;
θ = factor de ponderación; y
σ* = número de onda adimensional.
Subíndices
o = flujo de equilibrio; y
n = numérico.
Superíndice
' = componente perturbada; y
* = variable adimensional.
|
211231 |
Documents in Portable Document Format (PDF) require Adobe Acrobat Reader 5.0 or higher to view; download Adobe Acrobat Reader. |