El número de Courant, Hidráulica Computational, Victor M. Ponce

Inundación en el Río Chané, departamento de Santa Cruz, Bolivia (15 de Enero de 1990).

EL NÚMERO DE COURANT

Victor M. Ponce

Profesor Emérito de Ingeniería Civil y Ambiental

Universidad Estatal de San Diego, San Diego, California

14 de junio 2024

RESUMEN. En este artículo tratamos el número de Courant utilizado en la hidráulica computacional, particularmente en la modelación del problema de convección pura. Arrojamos luz sobre las propiedades de estabilidad y convergencia de una clase de esquemas numéricos explícitos lineales de la ecuación de convección pura. Observamos que tanto la estabilidad condicional como la incondicional son posibles en esquemas explícitos. Además, hemos demostrado que la existencia de difusión numérica, aparentemente debido a los errores de primer orden del esquema numérico, afecta el cálculo, haciendo que la solución dependa del tamaño de la malla y sea claramente difusiva para todos los casos en el que el número de Courant C sea diferente de 1.

1. INTRODUCCIÓN

El número de Courant es un concepto fundamental en el campo de la dinámica de fluidos computacional (DFC) y, por extensión, de la hidráulica computacional. El número conecta efectivamente las propiedades físicas y numéricas de un esquema dado, con aplicación específica al problema de convección, es decir, la advección de la mecánica de fluidos. El número de Courant se define como la relación entre una velocidad física u y la "velocidad de la malla" Δx /Δt. En la práctica, el número de Courant comúnmente se expresa de la siguiente manera: C = u (Δt /Δx) (Ponce, 2014).

El nombre Courant reconoce la obra de Richard Courant, renombrado matemático alemán-estadounidense y uno de los más tempranos contribuyentes al campo de las matemáticas aplicadas. El concepto de número de Courant tiene sus raíces en la condición de Courant-Friedrichs-Lewy, la cual se conoce ampliamente como condición CFL (Courant et al., 1928).

En este artículo realizamos un análisis del número de Courant, arrojando nueva luz sobre su comportamiento. Nos centramos específicamente en la solución numérica del problema de convección, revelando propiedades que pueden haber estado ocultas hasta ahora en aplicaciones convencionales de la dinámica de fluidos e hidráulica computacional.

2. EL PROBLEMA DE CONVECCIÓN Y EL ESQUEMA GENERAL DE CUATRO PUNTOS

Comenzamos por conceptualizar el problema de convección de la mecánica de fluidos. Una cantidad física ζ se moverá, en el espacio x y en el tiempo t, con una cierta velocidad convectiva u siguiendo la ecuación diferencial parcial de primer orden (Ponce et al., 1979):

      ∂ζ              ∂ζ
   ^_____ + u ^_____ = 0
      ∂t              ∂x           (1)

Esta ecuación se puede resolver numéricamente de varias maneras. Una formulación lineal explícita de diferencias finitas está representada por el siguiente esquema (Fig. 1):

X (ζ_jⁿ⁺¹ - ζ _jⁿ) + (1-X) (ζ_j+1ⁿ⁺¹ - ζ_j+1ⁿ) Y (ζ_j+1ⁿ - ζ_jⁿ) + (1-Y)(ζ_j+1ⁿ⁺¹ - ζ_jⁿ⁺¹)
[ ^{_____________________________________}] + u [^{_____________________________________}] = 0
Δt Δx

(2)

Fig. 1 Discretización de la ecuación de convección (Ec. 1) siguiendo la Ec. 2.

En la Ecuación 2, resolviendo para el valor no conocido ζ_j+1ⁿ⁺¹:

ζ_j+1ⁿ⁺¹ = C₁ ζ_jⁿ + C₂ ζ_jⁿ⁺¹ + C₃ ζ_j+1ⁿ
(3)

Los coeficientes se definen como sigue:

                     (X + CY)
C₁ = ^{_____________________}
             (1 + C) - (X + CY) (4)

                C - (X + CY)
C₂ = ^{____________________}
            (1 + C) - (X + CY) (5)

                1 - (X + CY)
C₃ = ^{____________________}
            (1 + C) - (X + CY) (6)

en los cuales el numero de Courant C se define como sigue:

Δt
C = u ^_____
Δx (7)

Dado un número de Courant, habiéndose establecido previamente los factores de ponderación X e Y, y utilizando las Ecs. 2 a 7, el cálculo procede a calcular los valores de ζ para todos y cada uno de los puntos de la malla.

La Ecuación 3, junto con las Ecuaciones 4 a 7, constituye una familia de esquemas numéricos de la ecuación de convección, Ec. 1. Sus propiedades de estabilidad y convergencia dependen de los valores de X, Y, y C. Los dos primeros, denominados factores de ponderación, y C, el número de Courant, caracterizan el esquema elegido. Reiterando, el número de Courant es la relación entre la velocidad convectiva u y la relación de la malla, o "velocidad de la malla" Δx /Δt..

3. DIFUSIÓN Y DISPERSIÓN NUMÉRICA

La Ecuación 1 describe la convección, un proceso de primer orden; por lo tanto, no tiene en cuenta la difusividad, un proceso de segundo orden, y la dispersividad, un proceso de tercer orden (Ponce, 2023). Sin embargo, la Ec. 2, una solución numérica de la Ec. 1, normalmente mostrará un comportamiento difusivo y/o dispersivo debido al tamaño finito de la malla. Esto se atribuye a la difusión numérica (el error del término de primer orden) y a la dispersión numérica (el error del término de segundo orden).

Una medida de la cantidad de difusión y dispersión numéricas se puede obtener a partir del error de aproximación del esquema de diferencias finitas, la Ec. 2. Este error se obtiene expandiendo la función de la malla ζ(jΔx, nΔt) en serie de Taylor alrededor del punto (jΔx, nΔt), lo cual lleva a lo siguiente (Cunge, 1969; Ponce et al., 1979):

                     1                      1             ∂² ζ
R = u Δx [( ^____ - X ) + C (^____ - Y )]^_____
                     2                      2              ∂x²

                                1                        1                    ∂³ ζ
+ u Δx² { (1 - C ) [ ^____ ( X + CY) - ^____ (1 + C ) ] } ^_____ + ο(Δx³)
                                2                        3                      ∂x³ (8)

en la cual R = el error de aproximación. De la Ec. 8, observamos que para X = Y = 1/2, el esquema tiene una precisión de segundo orden. es decir, R = ο(Δx²). Además, para C = 1, la precisión es de tercer orden; es decir, la dispersión numérica desaparece. De hecho, para X = Y = 1/2 y C = 1 se obtiene la solución exacta de la ecuación de convección (Ec. 1).

En la Ecuación 8, el coeficiente del término de segundo orden ∂²ζ /∂x² se denomina coeficiente de difusión numérica y se define de la siguiente manera:

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

Los valores positivos de μ_n hacen que el esquema se difusione numéricamente; los valores negativos amplifican numéricamente. Mientras que en el primer caso la convergencia se ve perjudicada, en el segundo la estabilidad se ve afectada. Tanto la estabilidad como la convergencia mejoran cuando μ_n → 0. La Tabla 1 resume las propiedades de estabilidad y convergencia de tres esquemas lineales explícitos de primer orden.

Tabla 1. Propiedades numéricas de tres esquemas lineales explícitos de primer orden (Ec. 2) de la ecuación de convección pura (Ec. 1).
Descripción del esquema	X	Y	*μ_n*	Condición de estabilidad	Condición de convergencia (μ_n → 0)
(1)	(2)	(3)	(4)	(5)	(6)
I. Hacia adelante en el tiempo, hacia atrás en el espacio	0	1	(1/2) u Δx (1 - C)	C ≤ 1	C → 1
II. Hacia adelante en el espacio, hacia atrás en el tiempo	1	0	(1/2) u Δx (C - 1)	C ≥ 1	C → 1
III. Hacia atrás en el espacio, hacia atrás en el tiempo	0	0	(1/2) u Δx (1 + C)	Cualquier C	Ningún C

I. Esquema hacia adelante en el tiempo y hacia atrás en el espacio.

El esquema hacia adelante en el tiempo y hacia atrás en el espacio se centra en torno al nodo ζ_j+1ⁿ en la Fig. 1.
La derivada temporal se toma en el lado derecho del cuadro y la derivada espacial en el lado inferior.
El coeficiente de difusión numérica es: μ_n = (1/2) u Δx (1 - C).
El esquema es estable para C ≤ 1, e inestable para C > 1.
El esquema converge a la solución analítica de la ecuación de convección (Ec. 1) en cuanto C se aproxima a 1 (desde un valor menor que 1).
Cuanto más refinado sea el tamaño de la cuadrícula (conforme Δx ↠ 0), más rápida será la convergencia.
El esquema es condicionalmente estable, volviéndose inestable para C > 1.

II. Esquema hacia adelante en el espacio y hacia atrás en el tiempo.

El esquema hacia adelante en el espacio y hacia atrás en el tiempo se centra en torno al nodo ζ_jⁿ⁺¹ en la Fig. 1.
La derivada temporal se toma en el lado izquierdo del cuadro y la derivada espacial en el lado superior.
El coeficiente de difusión numérica es: μ_n = (1/2) u Δx (C - 1).
El esquema es estable para C ≥ 1, e inestable para C < 1.
El esquema converge a la solución analítica de la ecuación de convección (Ec. 1) en cuanto C se aproxima a 1 (desde un velor mayor que 1).
Cuanto más refinado sea el tamaño de la cuadrícula (conforme Δx ↠ 0), más rápida será la convergencia.
El esquema es condicionalmente estable, volviéndose inestable para C < 1.

III. Esquema hacia atrás en el espacio y hacia atrás en el tiempo.

El esquema hacia atrás en el espacio y hacia atrás en el tiempo se centra en torno al nodo ζ_j+1ⁿ⁺¹ en la Fig. 1.
La derivada temporal se toma en el lado derecho del cuadro y la derivada espacial en el lado superior.
El coeficiente de difusión numérica es: μ_n = (1/2) u Δx (1 + C).
El esquema es incondicionalmente estable, es decir, estable para cualquier valor de C.
Sin embargo, para cualquier valor de C, este esquema no converge a la solución analítica del problema de convección, inclusive cuando Δx → 0.

El presente artículo permite llegar a las siguientes conclusiones:

Los esquemas explícitos de la ecuación de convección pura (Ec. 1), tales como la Ec. 2, pueden ser condicional o incondicionalmente estables, dependiendo de la forma en que se formule el esquema. Por lo tanto, la afirmación de que todos los esquemas explícitos están sujetos a una condición de estabilidad es incorrecta.
Dado que el tamaño de la cuadrícula es siempre finito, los esquemas explícitos como los de la Ec. 2 introducen cantidades variables de difusión numérica. Sólo para el número de Courant C = 1 la difusión numérica desaparece; por lo tanto, en este caso la solución numérica del esquema representado por la Ec. 2 es la solución exacta de la ecuación de convección (Ec. 1).
En todos los demás casos, es decir, específicamente para C < 1 en el esquema hacia adelante en el tiempo y hacia atrás en el espacio, o C > 1 en el esquema hacia atrás en el tiempo y hacia adelante en el espacio, la solución numérica resulta ser efectivamente una onda de difusión, con cantidades variables de difusión numérica (Cunge, 1969; Ponce, 2023).

Al contrario de los esquemas explícitos, los esquemas implícitos se han visto favorecidos debido a sus percibidas mejores propiedades de estabilidad, supuestamente presentando lo que a menudo se ha denominado "estabilidad incondicional". La experiencia, sin embargo, demuestra lo contrario. En general, tanto los esquemas numéricos explícitos como los implícitos están sujetos a situaciones (condiciones) de estabilidad y convergencia. Se recomienda un análisis de von Neumann para mejorar la práctica de la modelación numérica (Ponce et al., 1978).

4. RESUMEN

En este artículo detallamos el papel del número de Courant en la hidráulica computacional, particularmente en el modelado del problema de convección pura. Arrojamos una luz renovada sobre las propiedades de estabilidad y convergencia de una clase de esquemas numéricos lineales explícitos de la ecuación de convección pura. Se ha demostrado que tanto la estabilidad condicional como la incondicional son posibles con esquemas explícitos. Además, se ha demostrado que la existencia de difusión numérica, atribuible a los errores de primer orden del esquema numérico, afecta el cálculo, haciendo que la solución dependa del tamaño de la malla y sea claramente difusiva para todos los casos distintos del número de Courant C = 1. Se recomienda el análisis de von Neumann para mejorar la práctica de modelado.

REFERENCIAS

Courant, R., K. Friedrichs, Y H. Lewy. 1928. Über die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen, vol. 100, p. 32-74.

Cunge, J. A. 1969. On the Subject of a Flood Propagation Computation Method (Muskingum Method). Journal of Hydraulic Research, Vol. 7, No. 2, 205-230.

Ponce, V. M., H. Indlekofer, y D. B. Simons. 1978. Convergence of four-point implicit water wave models. Journal of the Hydraulics Division, ASCE, Vol. 104, No. HY7, July, 947-958.

Ponce, V. M., Y. H. Chen, y D. B. Simons. 1979. Unconditional stability in convection computations. Journal of the Hydraulics Division, ASCE, Vol. 105, No. HY9, September, 1079-1086.

Ponce, V. M. 2014. Engineering Hydrology: Principles and Practices. Second edition, online.

Ponce, V. M. 2023. Is dispersion important in flood routing? Online paper.

Ponce, V. M. 2023. When is the diffusion wave applicable? Online paper.

240614