1.  INTRODUCCIÓN

Las propiedades de estabilidad y convergencia de los esquemas numéricos explícitos de la ecuación de convección pura se conocen desde hace algún tiempo. Aquí se presenta un tratamiento teórico unificado. Los análisis de von Neumann (4) y Hirt (2) se utilizan para mostrar cómo la estabilidad incondicional y la precisión de segundo orden son posibles dentro del marco de un esquema explícito (6).

2.  ECUACIÓN DE CONVECCIÓN PURA

La ecuación de convección pura tiene la siguiente forma:

∂ζ              ∂ζ                   _____  +  u _____  =  0       ∂t              ∂x

(1)

en la cual ζ = cantidad física a ser convectada con la velocidad u. Un esquema general explícito de diferencias finitas de esta ecuación es el siguiente (Fig.1):

X (ζj n+1-ζ jn) + (1-X) (ζ j+1 n+1-ζj+1n)              Y (ζ j+1 n-ζ j n) + (1-Y)(ζj+1n+1-ζ jn+1) [ ____________________________________ ] + u [ ___________________________________]  =  0                              Δt                                                                   Δx

(2)

en la cual X y Y = factores de ponderación de las derivadas temporal y espacial, respectivamente; Δt y Δx = incrementos temporales y espaciales, respectivamente; y u = velocidad convectiva promedio en una ubicación particular de la malla.

Fig. 1. Esquema de discretización de la ecuación de convección pura.

En la Ecuación 2, resolviendo explícitamente el valor desconocido ζ _j+1 ⁿ⁺¹ :

ζ j+1n+1 = C1 ζjn + C2 ζj n+1 + C3 ζ j+1n

(3)

en la cual:

(X + CY)         C1 =   ___________________            (1 + C) - (X + CY)

(4)

C - (X + CY)         C2 =   ____________________             (1 + C) - (X + CY)

(5)

1 - (X + CY)         C3 =  ___________________            (1 + C) - (X + CY)

(6)

y

Δt         C =   u _____               Δx

(7)

es el número de Courant correspondiente a una ubicación particular de la malla.

3.  COEFICIENTE DE DIFUSIÓN NUMÉRICA

La Ecuación 1 es una ecuación no viscosa; por lo tanto, la viscosidad (o difusividad, según el caso) está ausente. Sin embargo, una solución numérica de la ecuación, a menudo muestra un comportamiento viscoso (difusivo). Esto debe atribuirse correctamente a la viscosidad artificial del propio esquema. El término "dispersión numérica" se utiliza a veces, aparentemente para indicar el hecho de que tanto los errores de amplitud como los de fase se agregan para producir un efecto de tipo dispersivo en la solución numérica.

Se puede obtener una medida de la cantidad de dispersión numérica a partir del error de aproximación del esquema de diferencias finitas. Esto se puede obtener expandiendo la función de malla ζ(jΔx, nΔt) en una serie de Taylor sobre el punto jΔx, nΔt, lo cual conduce a lo siguiente (1,2):

1                     1             ∂2 ζ R = u Δx [( ____ - X) + C ( ____ - Y )] _____                     2                     2              ∂x2                                1                        1                     ∂3 ζ + uΔx2 { (1 - C ) [ ____ ( X + CY) - ____ (1 + C ) ] } _____ + ο(Δx3)                                2                        3                     ∂x3

(8)

en la cual R = error de aproximación. A partir de la Ec. 8, para X = Y = 1/2, el esquema es de precisión de segundo orden, es decir, R = ο(Δx²). Además, para C = 1, la precisión es de tercer orden. De hecho, para estas condiciones se obtiene la solución exacta (1).

El coeficiente del término ∂²ζ /∂x² en la Ec. 8 se conoce como el coeficiente de difusión numérica, μ_n, definido como sigue:

1                          1 μn = u Δx [(____  - X )  +   C (____  - Y )]                     2                          2

(9)

Los valores positivos de μ_n hacen que en la amplitud de la onda haya difusión numérica; los valores negativos hacen que la amplitud de la onda se amplifique numéricamente. Mientras que en el primero la consistencia se ve afectada, en el segundo la estabilidad se ve afectada. Tanto la estabilidad como la consistencia mejoran a medida que μ_n tiende a cero. La Tabla 1 resume las propiedades de estabilidad y convergencia (consistencia en el límite) de tres esquemas usuales de primer orden.

4.  ANÁLISIS LINEAL DE VON NEUMANN

Se puede obtener una visión amplia de las propiedades numéricas de los diversos esquemas explícitos de cuatro puntos para el cálculo de convección pura mediante un análisis lineal de estabilidad y convergencia siguiendo el enfoque de von Neumann (3,4). De acuerdo con este método, se busca una solución de la Ec. 2 en la siguiente forma sinusoidal:

ζ = ζο exp [i ( σx - βt )]

(10)

en la cual ζ_ο = función de amplitud; σ = número de onda de la solución sinusoidal; y β = factor complejo de propagación, de tal manera que β = β_R + i β_I. La cantidad β_R = frecuencia de onda; y β_I = factor de amplificación (o atenuación). La sustitución de la Ec. 10 en la Ec. 2 conduce a lo siguiente:

[ (1 - X) exp ( i σ Δx) + X ] - CY [exp (i σ Δx) - 1] exp (- i β Δt ) = _______________________________________________________                            [(1 - X) exp ( i σ Δx) + X ] + C(1 - Y) [exp ( i σ Δx ) - 1]

(11)

Sustituyendo exp(i σ Δx) = p + i q, en la cual p = cos (σΔx) y q = sin (σΔx) en la Ec. 11 conduce a:

[p (1 - R ) + R ] + iq (1 - R ) exp (- i β Δt ) = _______________________________                               [p (1 - S ) + S] + iq (1 - S )

(12)

en la cual:

R = X + CY

(13)

S = R - C

(14)

Operando en la Ec. 12 se obtiene:

M + i N exp (- i β Δt ) = 1 - __________                                     P

(15)

en la cual:

M = (1 - p) (1 - S ) C

(16)

N = q C

(17)

y

P = 1 - 2 ( 1- p) (1 - S ) S

(18)

Separando la Ec. 15 en sus componentes reales e imaginarios:

N             tan (- βR Δt ) =   __________                              P - M

(19)

[( P - M)2 + N 2 ] 1/2            exp (- βI Δt ) = _____________________                                        P

(20)

Dado que la solución analítica de la Ec. 1 no exhibe atenuación, se deduce que la Ec. 20 es una medida de la atenuación (o amplificación) numérica de la solución numérica, después de un incremento de tiempo Δt. Una relación de convergencia con respecto a la amplitud de onda se define de la siguiente manera:

R1 = exp (βI Δt)

(21)

De la Ecuación 19, la velocidad de fase de la solución numérica es:

βR           1                        N un = ______  = ______ tan -1   ( _________ )             σ          σΔt                   P - M

(22)

Una relación de convergencia con respecto a la fase de onda se define de la siguiente manera:

un             1                        N R2 = ______  = _______ tan -1   ( _________)             u          σuΔt                   P - M

(23)

o igualmente:

1                            N R2  = _________   tan -1   ( _______ )           C (σΔx)                   P - M

(24)

en la cual σΔx = resolución espacial. Debido a σ = 2π/L, en la cual L es la longitud de onda sinusoidal, se deduce lo siguiente:

2 π σΔx = _______                L              ____               Δx

(25)

en la cual L / Δx = número de intervalos discretos de espacio por longitud de onda sinusoidal. Ambos R₁ y R₂ son funciones de X, Y, C, y L / Δx.

La Figura 2 muestra la variación de R₁ y R₂ como una función de X, Y, C, y L / Δx. Se pueden obtener las siguientes conclusiones:

1. Para C = 1, la coincidencia exacta de soluciones numéricas y analíticas R₁ = R₂ = 1 se da para valores de X e Y a lo largo de la diagonal que une X = 0, Y = 1 , y X = 1, Y = 0.

2. Para C ≠ 1, no hay combinación de X y Y para el cual R₁ = R₂ = 1.

Fig. 2. Razones de convergencia R1 y R2 como una función de X, Y y C: (a) L / Δx = 10; (b) L / Δx = 40.

La inestabilidad, R₁ > 1, del Esquema I (Tabla I) para números de Courant > 1, y del Esquema II para números de Courant < 1, es confirmada. Asimismo, se confirma la estabilidad incondicional y una mayor dispersión numérica del Esquema III.

5.  PROPIEDADES DE ESTABILIDAD Y CONVERGENCIA

Ecuación lineal con coeficiente constante. En este caso, la velocidad convectiva u es constante en el espacio y el tiempo. La solución numérica usando las Ecs. 3-6 con X = Y = 1/2 y C = u Δt / Δx = 1 reproducirá exactamente la solución analítica.

Ecuación lineal con coeficiente variable. En este caso, la velocidad convectiva u varía en el espacio y el tiempo. Por lo tanto, el número de Courant varía de una celda a la siguiente. El análisis de estabilidad lineal de von Neumann muestra que la relación de convergencia (amplitud de onda) R₁ es igual a 1 para X = Y = 1/2 y para todos los valores de C. Sin embargo, la relación de convergencia (fase de onda) R₂ es igual a 1 solo para C = 1. Para números de Courant diferentes de uno, siempre habrá una cierta cantidad de error de fase.

La estabilidad del esquema requiere que μ ≥ 0; por otra parte, la convergencia de la amplitud de la onda mejora a medida que μ_n → 0. Usando la Ec. 9, se pueden formular varios esquemas explícitos de segundo orden, incondicionalmente estables pero precisos. Las propiedades numéricas de tres de estos esquemas se dan en la Tabla 2. Los Esquemas IV y V pueden considerarse como generalizaciones de los Esquemas I y II, respectivamente, para el caso de números de Courant variables. Sin embargo, a diferencia de los esquemas I y II, los esquemas IV y V no difusionan numéricamente la amplitud de onda. El esquema VI, un esquema centrado en el espacio y el tiempo, tiene una precisión de segundo orden para todos los valores de C.

Para valores fijos de X y Y, el error de segundo orden de los Esquemas IV a VI será directamente proporcional a | 1 - C |. Asimismo, para valores fijos de X, Y y C (C ≠ 1), el error de segundo orden será directamente proporcional a ∂³ζ /∂x³. La equivalencia de los Esquemas IV a VI se puede demostrar sustituyendo los valores de X y Y de la Tabla 2 en las Ecs. 4-6. Esto lleva a:

C1 = 1

(26)

1 - C C2 =     ________               1 + C

(27)

y

1 - C C3 =     ________               1 + C

(28)

para los tres esquemas. Además, para C = 1, C₂ = C₃ = 0, y ζ _{j + 1} ^{n + 1} = ζ_j ⁿ, es decir, la cantidad ζ está siendo sometida a convección pura. Para C ≠ 1 se introduce un efecto dispersivo porque el cálculo de ζ _{j + 1} ^{n + 1} se basa no solamente en ζ_j ⁿ sino también en ζ_j ^{n + 1} y ζ_{j + 1} ⁿ .

6.  EJEMPLO ILUSTRATIVO

Se estableció un ejemplo hipotético para demostrar las propiedades numéricas de los esquemas descritos en la Tabla 2. Un canal de 20,000 pies de largo se dividió en 20 tramos de Δx = 1,000 pies (304.8 m) cada uno. Una distribución sinusoidal ζ, la cual varía de 0 a 100 unidades, es especificada en el eje de tiempo como condición de frontera aguas arriba. La distribución se especifica en 20 intervalos de Δt = 1,000 segundos cada uno, y la simulación se lleva a cabo durante un total de 50 intervalos de tiempo, es decir, 50,000 segundos. Una velocidad convectiva que varía linealmente en el espacio entre u = 0.5 pies/s (0.15 m/s) en el límite aguas arriba, y u = 1.5 pies/s (0.45 m/s) en el límite aguas abajo.

Los resultados de la simulación se muestran en la Fig. 3. La estabilidad, la convergencia de la amplitud de la onda y la conservación de la masa se mantienen para todos los números de Courant. Como era de esperar, se detectó una pequeña cantidad de ruido numérico.

Fig. 3. Resultados de un ejemplo hipotético.

Esto es causado por los errores de segundo orden y tiende a concentrarse en las regiones de cambios bruscos en ∂³ζ /∂x³ (o ∂³ζ /∂t ³). El error numérico es una función de | 1 - C | y desaparece conforme C → 1.

7.  RESUMEN

Se presenta un tratamiento teórico unificado de las propiedades numéricas de una clase de esquemas explícitos de la ecuación de convección pura. Para los casos de números de Courant que varían lentamente, la teoría se usa para mostrar que la estabilidad incondicional y la precisión de segundo orden pueden lograrse dentro del marco de una formulación explícita.

APÉNDICE I.  BIBLIOGRAFÍA

1. Cunge, J. A. 1969. "On the subject of a flood propagation computation method (Muskingum method)," Journal of Hydraulic Research, Vol. 7, No. 2, pp. 205-230.

2. Hirt, C. W. 1963. "Heuristic stability theory for finite difference equations," Journal of Computational Physics, Vol. 2, pp. 339-355.

3. Leendertse, J. J. 1967. "Aspects of a computational model for long-period water-wave propagation," Memorandum RM-5294-PR, May, The Rand Corporation, Santa Monica, Calif.

4. O'Brien, G. G., M. A. Hyman, and S. Kaplan. 1950. "A study of the numerical solution of partial differential equations," Journal of Mathematics and Physics, Vol. 29, No. 4, pp. 231-251.

APÉNDICE II.  SIMBOLOGÍA

C = número de Courant;

C₁ = coeficiente, Ec. 4;

C₂ = coeficiente, Ec. 5;

C₃ = coeficiente, Ec. 6;

j = índice de discretización espacial;

L = longitud de onda;

M = parámetro, Ec. 16;

N = parámetro, Ec. 17;

n = índice de discretización de tiempo;

P = parámetro, Ec. 18;

p = cos (σΔx);

q = sin (σΔx);

R = error de aproximación; Además, parámetro, Ec. 13;

R₁ = relación de convergencia de amplitud de onda;

R₂ = relación de convergencia de la fase de onda;

S = parámetro, Ec. 14;

u = velocidad media;

μ_n = velocidad de fase de la solución numérica;

X = factor de ponderación de la derivada temporal;

Y = factor de ponderación de la derivada espacial;

β_I = factor de amplificación (o atenuación);

β_R = frecuencia de onda;

Δt = incremento temporal;

Δx = incremento espacial;

ζ = alguna cantidad física a ser convectada;

μ_n = coeficiente de difusión numérica, Ec. 9; y

σ = número de onda.