ESQUEMAS CONSERVATIVOS BASADOS EN EL MÉTODO DE SEPARACIÓN, PARA LA SIMULACIÓN NUMÉRICA DE VÓRTICES EN LA ATMÓSFERA

Yuri N. Skiba y Denis M. Filatov

Yuri N. Skiba. Doctor en Física y Matemáticas, Academia de Ciencias de USSR, Novosibirsk, Rusia. Profesor Investigador Titular, Centro de Ciencias de la Atmósfera, Universidad Nacional Autónoma de México (CCA-UNAM), México. Dirección: CCA-UNAM, C.P. 04510, México DF, México. e-mail: skiba@servidor.unam.mx

Denis M. Filatov. Doctor en Ciencias de la Computación, Centro de Investigación en Computación, Instituto Politécnico Nacional, México Profesor Investigador Asociado, CCA-UNAM. e-mail: denisfilatov@mail.ru

Resumen

Se propone un nuevo método para la construcción de esquemas en diferencias finitas para el modelo de aguas someras en un canal periódico sobre el plano y sobre la esfera. La ventaja principal de los esquemas consiste en que ellos conservan exactamente la masa y la energía total del sistema. El enfoque está basado en la separación del operador del modelo por coordenadas y por procesos físicos. Como resultado, la solución del problema original se reduce a la solución de problemas simples unidimensionales. Se propone un conjunto de tales esquemas, que son lineales o no lineales dependiendo de la especificación de parámetros del esquema. A diferencia de los métodos conocidos, este enfoque permite derivar esquemas numéricos conservadores después de hacer discreto el problema original tanto en tiempo como en espacio. Además, los algoritmos numéricos utilizados para calcular la solución son económicos desde el punto de vista computacional, ya que cada esquema se realiza rápida y fácilmente usando un método directo, sin iteraciones.

Summary

A new method for the construction of finite difference schemes for the shallow water equations studied in periodic plane and spherical channels is proposed. The main advantage of the schemes is that they conserve exactly the mass and the total energy of the system. The approach is based on the method of splitting the original operator by coordinates and by physical processes. Thereby, the solution to the original complex problem is reduced to the solution of simple one-dimensional problems. A set of such schemes is proposed, which are either linear or nonlinear, depending on the choice of parameters of the scheme. Unlike existing methods in the field, this approach allows deriving conservative finite difference schemes when the equations are considered to be discrete both in time and in space. Furthermore, the numeric algorithms used to compute the solution are computationally economic, because each scheme is easily implemented by using fast direct methods of linear algebra.

Resumo

Se propõe um novo método para a construção de esquemas em diferenças finitas para o modelo de águas superficiais em um canal periódico sobre o plano e sobre a esfera. A vantagem principal dos esquemas consiste em que eles conservam exatamente a massa e a energia total do sistema. O enfoque está baseado na separação do operador do modelo por coordenadas e por processos físicos. Como resultado, a solução do problema original se reduz à solução de problemas simples unidimensionais. Se propõe um conjunto de tais esquemas, que são lineares ou não lineares dependendo da especificação de parâmetros do esquema. A diferença dos métodos conhecidos, este enfoque permite derivar esquemas numéricos conservadores depois de fazer discreto o problema original tanto em tempo como em espaço. Além disso, os algoritmos numéricos utilizados para calcular a solução são econômicos desde o ponto de vista computacional, já que cada esquema se realiza rápida e facilmente usando um método direto, sem iterações.

La capa inferior de la atmósfera (tropósfera) tiene la altura promedio de 10km. Dado que las escalas horizontales de los procesos atmosféricos son de miles de km, la dinámica global de la atmósfera puede ser descrita aproximadamente por el modelo clásico de aguas someras (Arakawa y Lamb, 1981; Vreugdenhil, 1994). En coordenadas cartesianas las ecuaciones de aguas someras (EAS) son

donde u = u(x,y,t) y v = v(x,y,t) son los componentes de la velocidad de un fluido, H = H(x,y,t) es la profundidad del fluido, f = f(y) es la aceleración de Coriolis causada por la rotación de la Tierra, y h = h(x,y,t) = H(x,y,t) + h_r(x,y) es la altura de la superficie libre, donde h_r es la altura del relieve (Figura 1). Las EAS se resuelven con ciertas condiciones iniciales y de frontera. Se considerará el problema (1.1)-(1.3) en un canal con la condición periódica en x y con el componente normal de la velocidad nulo, (u,v)_n = 0 sobre la frontera lateral (Figura 2).

Es conocido que las EAS, como un sistema físico cerrado, poseen varias leyes de conservación. En particular, ellas conservan importantes características integrales como son la masa, la energía total, la enstrofía potencial, etc. (Vreugdenhil, 1994). El problema es que para la simulación numérica el sistema (1.1)-(1.3) se aproxima por un modelo discreto (o bien, por un esquema en diferencias finitas), y generalmente la discretización en espacio y en tiempo destruye todas o la mayoría de las leyes de conservación. En otras palabras, la masa, la energía total y la enstrofía en su forma discreta no se conservan al pasar de un momento de tiempo a otro. Sin embargo, es muy deseable conservar las leyes también en la forma discreta, ya que dichas leyes no sólo garantizan la estabilidad de los algoritmos numéricos aplicados, sino también la correcta distribución de energía de la solución entre diferentes frecuencias (Vreugdenhil, 1994) y, por lo tanto, los resultados fieles de la simulación.

Desde los años 60 han sido propuestos varios esquemas numéricos en diferencias finitas que conservan unas u otras características integrales de las EAS (Sadourny, 1975; Arakawa y Lamb, 1981; Takano y Wurtele, 1982; Kim, 1984; Ringler y Randall, 2002; Salmon, 2004). Por ejemplo, los algoritmos descritos en Ringler y Randall (2002) y Sadourny (1975) usan mallas geodésicas para hacer discreto el problema en espacio, mientras que el algoritmo propuesto en Salmon (2004) usa un patrón de malla bastante amplio, de 25 puntos. No obstante, a pesar de la indudable importancia de dichos algoritmos, cada uno de ellos posee leyes de conservación solo en la forma semidiscreta, cuando el problema ya es discreto en espacio y todavía es continuo en tiempo. Lamentablemente, la discretización en tiempo de un problema semidiscreto usando un esquema explícito siempre destruye las leyes de conservación (Marchuk, 1988). El único esquema que podría conservar las leyes después de discretizar el problema en tiempo es el esquema implícito de Crank-Nicolson (Crank y Nicolson, 1947; Marchuk, 1988). Sin embargo, su aplicación a menudo es demasiado difícil. Por ejemplo, al aplicar el esquema de Crank-Nicolson con cualquiera de los algoritmos desarrollados en Ringler y Randall (2002) y en Salmon (2004), la matriz del sistema lineal resulta densa (multidiagonal) y para resolverlo se requiere aplicar un método iterativo, lo que destruye de nuevo las leyes de conservación y aumenta drásticamente el costo computacional. Además, cabe enfatizar que las mallas geodésicas tienen aplicación limitada; son buenas solo para toda la esfera, pero es difícil usarlas si la frontera del dominio es de una forma arbitraria. Por lo anterior, es de interés desarrollar algoritmos discretos en espacio y en tiempo, que sean simples en la realización y posean las leyes de conservación.

En este trabajo se desarrolla un nuevo método para discretizar y resolver numéricamente las EAS. El modelo obtenido después de discretizar las EAS en espacio y en tiempo conserva exactamente la masa y la energía total. El enfoque emplea la técnica de separación del operador del modelo por coordenadas y por procesos físicos, con el fin de reducir la solución del problema bidimensional (2D) original a la solución de dos problemas unidimensionales (1D) de onda y de un problema simple (con dos ecuaciones diferenciales ordinarias) de rotación (Yanenko, 1971; Skiba, 1995). Cada uno de estos tres problemas continuos separados conserva la masa y la energía total. Los problemas discretos, derivados de éstos, también conservan exactamente la masa y la energía total. De hecho, se propone un conjunto de esquemas numéricos, los cuales son lineales o no lineales dependiendo de los parámetros del esquema (Skiba, 1995).

En la siguiente sección las EAS se transforman a la forma divergente y luego se aplica el método de separación. Después se derivan varios esquemas en diferencias finitas que conservan la masa y la energía total cuando el sistema se considera en un canal periódico sobre el plano. Posteriormente los resultados se extienden al caso esférico. Y, finalmente, se dan algunas pruebas numéricas.

Se parte de un cambio de variables (Skiba, 1995), definiendo

y se reescribe el sistema de las EAS (1.1)-(1.3) en la forma

En cada intervalo de tiempo [t_n,t_n + t], se aproxima el sistema (2.2)-(2.4) por tres sistemas más simples. Primero se consideran las ecuaciones en la dirección x dejando la coordenada y fija:

Después se intercambian las coordenadas y se considera el sistema con respecto a la dirección y bajo la x fija:

Finalmente, se trata separadamente la rotación (el término de Coriolis):

Los sistemas (2.5)-(2.7) y (2.8)-(2.10) son unidimensionales, mientras que el sistema (2.11)-(2.12) es un sistema simple de dos ecuaciones diferenciales ordinarias. Resolviendo consecutivamente cada uno de los tres sistemas en el mismo intervalo [t_n,t_n + t], finalmente se obtiene una solución que aproxima la solución del problema no separado original (2.2)-(2.4). La precisión de la aproximación se mejora cuando el paso t de la malla temporal tiende a cero.

Se construyen ahora (Skiba, 1995) los esquemas conservadores discretizando los sistemas 1D (2.5)-(2.7) y (2.8)-(2.10). Se introduce la notación t = t_n+1–t_n, Dx = x_k+1 –x_k, Dy = y_l+1 – y_l, Rⁿ_kl = R(x_k, y_l, t_n), , donde R representa cualquiera de las funciones u, v, U, V, H, h y z. Entonces, el primer sistema se aproxima como

omitiéndose el subíndice l para simplificar las denotaciones. A su vez, el segundo sistema discreto tiene la forma

omitiéndose también el subíndice k.

El último sistema, el (2.11)-(2.12), se aproxima como

Las ecuaciones (3.1)-(3.8) representan, de hecho, una familia de esquemas que aproximan el sistema original no lineal (2.5)-(2.12) de las EAS. Esta familia depende de los parámetros Por ejemplo, el esquema (3.1)-(3.8) es no lineal si

Es posible demostrar que el sistema separado (3.1)-(3.3) tiene dos leyes de conservación. Multiplicando (3.3) por DxDyt, sumando sobre todos los puntos interiores del dominio de interés, y tomando en cuenta las condiciones de frontera se tiene que

lo que muestra que la masa no se cambia en el tiempo. También, multiplicando (3.1) por DxDytU_k, (3.2) por DxDytV_k y (3.3) por DxDytgh_k, y después sumando sobre todos los puntos interiores, se obtiene

es decir, la energía total (la cinética más la potencial) se conserva. Similarmente se puede demostrar que los sistemas (3.4)-(3.6) y (3.7)-(3.8) también conservan las mismas características. Ya que cada uno de los esquemas separados (3.1)-(3.3), (3.4)-(3.6) y (3.7)-(3.8) conserva la masa y la energía total, el esquema completo (3.1)-(3.8) también conserva dichas características integrales (Skiba, 1995). El esquema posee las dos leyes de conservación, independientemente de la especificación de los parámetros

Otra importante ventaja del método de separación consiste en que éste permite reducir la solución del problema original complejo a la solución de problemas más simples, lo que se hace fácilmente usando algoritmos rápidos del álgebra lineal (Skiba y Filatov, 2005). Además, la implementación numérica de los esquemas se hace por métodos directos; es decir, no requiere el uso de métodos iterativos, y por lo tanto las leyes de conservación no se destruyen en la etapa de cálculos, a diferencia de Ringler y Randall (2002) y Salmon (2004), entre otros.

Los resultados obtenidos admiten una generalización al caso de la geometría esférica. Consideremos el modelo de aguas someras sobre la esfera. En las coordenadas esféricas las EAS (2.2)-(2.4) tienen la forma (Batchelor, 1970; White, 2002)

donde f = 2Wsin j, siendo W la velocidad angular de rotación de la esfera, a el radio de la esfera (de la Tierra), lÎ[0,2p] la longitud jÎ[-b,b] la latitud. Nótese que b< p/2 cuando se considera un canal periódico y b= p/2 para toda la esfera.

La separación del sistema (4.1)-(4.3) se hace de manera similar a la del sistema (2.2)-(2.4), y después de la discretización en espacio y en tiempo se obtienen los tres sistemas discretos siguientes:

en la dirección l;

en la dirección j ; y

para la rotación de la esfera.

En las Figuras 3 y 4 se muestra la idea del método de separación por coordenadas para el caso esférico: primero se soluciona el sistema (4.4)-(4.6) en la dirección l, después se soluciona el sistema (4.7)-(4.9) en la dirección j. Con esto, es muy importante enfatizar que cuando entonces el denominador de la expresión es igual a cero, lo que resulta en que las ecuaciones no tienen sentido en los polos de la esfera. Por eso, para poder usar el sistema (4.1)-(4.3) y, por consiguiente, los sistemas (4.4)-(4.6) y (4.7)-(4.9) también, sobre toda la esfera, movemos la malla en la dirección j a medio paso Dj/2. Esto permite eliminar los puntos de polos de la consideración.

Es fácil verificar que el esquema (4.4)-(4.11) también conserva las dos características integrales: la masa

Como condiciones iniciales, para el parámetro h se consideró una función del tipo "sombrero", mientras para u y v se tomó un campo uniforme de velocidades perturbado por una estructura turbulenta. Tal especificación de los parámetros h, u y v cualitativamente describe un vórtice real observable en la atmósfera. A saber, en un ciclón atmosférico el parámetro h tiene la forma de función "sombrero" positiva, mientas que la estructura turbulenta es responsable de la circulación de éste.

El canal periódico fue el rectángulo [0,0; 1,0] ´ [0,0; 0,25], con los pasos de la malla Dx = 0,005, Dy = t = Dx/2.

En la Figura 5 se muestran los gráficos del comportamiento temporal de la masa, la energía total y la enstrofía potencial. Se puede ver que las primeras dos características se conservan durante todo el tiempo de modelación, lo que está en consonancia con el estudio teórico. En cuanto a la enstrofía potencial, ella no es constante, pero sus oscilaciones son muy pequeñas, alrededor de un valor promedio (aproximadamente 181,5), que es un resultado prometedor. Además, como se sigue de la Figura 6, la estructura turbulenta inicial del campo de velocidades (arriba) se conserva bastante bien en tiempo, al comparar con la Figura 6, abajo. Consideramos estos resultados como buenos, y ellos permiten concluir que el método desarrollado y los esquemas derivados son adecuados para la simulación (al menos, cualitativa) numérica de vórtices en la atmósfera. Los esquemas desarrollados pueden ser modificados en un trabajo posterior, con el fin de conservar la tercera importante característica integral de las EAS, la enstrofía potencial.

Como condiciones iniciales, para el parámetro h se consideró la misma función del tipo "sombrero", mientras para u y v se tomó un campo de velocidades uniforme. El canal periódico fue el conjunto [0,0; 2p] ´ [–p/4, p/4], los pasos de la malla fueron Dl = p/90, Dj = Dl/2, t = 0,001.

La Figura 7 muestra el comportamiento de la masa, la energía total, y la enstrofía potencial en tiempo. Puede verse que los gráficos son casi idénticos a los que se mostraron arriba para el caso cartesiano. Es un resultado previsto, ya que el sistema (4.1)-(4.3) no difiere mucho de las ecuaciones (2.2)-(2.4), y contiene solo el factor métrico adicional, a cos j.

Se desarrolló un nuevo método para la construcción de esquemas conservativos en diferencias finitas para el modelo de aguas someras. El modelo se consideró en un canal periódico sobre el plano y sobre la esfera. El enfoque presentado se basa en la separación del operador del modelo por coordenadas y por procesos físicos. La discretización en tiempo se realiza por el método de pasos fraccionados con el uso del esquema de Crank-Nicolson en cada paso, mientras que la discretización en espacio se realiza por el método de diferencias finitas. Se derivaron esquemas conservativos que son lineales o no lineales según la definición de sus parámetros. El método desarrollado tiene importantes ventajas: 1) A diferencia de los métodos conocidos que conservan la masa, la energía total y/o la enstrofía potencial solo en su forma semidiscreta, cada esquema conserva la masa y la energía total en su forma totalmente discreta, ambos en espacio y en tiempo. 2) La realización numérica de los esquemas es económica desde el punto de vista computacional, ya que cada esquema implica la solución de un sistema lineal con una matriz de cinco diagonales, lo que se resuelve rápidamente usando un método directo del álgebra lineal. El presente modelo, junto con los esquemas desarrollados, puede ser utilizado para la modelación numérica de vórtices atmosféricos (huracanes) permitiendo obtener resultados fieles y admite también la aplicación a los problemas oceánicos (tsunamis, ondas de marea, etc.)

El trabajo contó con el apoyo de los proyectos PAPIIT-UNAM IN105005 (México) y FOSEMARNAT-CONACyT 2004-01-160 (México), y por becas del Sistema Nacional de Investigadores (México). El segundo autor también agradece a la UNAM por una estancia y beca posdoctoral.

1. Arakawa A, Lamb VR (1981) A Potential Enstrophy and Energy Conserving Scheme for the Shallow-Water Equation. Month. Weather Rev. 109: 145-171.        [ Links ]

2. Batchelor GK (1970) An Introduction to Fluid Dynamics. Cambridge University Press. Cambridge, RU. 352 pp.        [ Links ]

3. Crank J, Nicolson P (1947) A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations of the Heat Conduction Type. Proc. Cambridge Philos. Soc. 43: 50-67.        [ Links ]

4. Kim VF (1984) Numerical Analysis of Some Conservative Schemes for Barotropic Atmosphere. Soviet Meteorol. Hydrol. 6: 251-261.        [ Links ]

5. Marchuk GI (1982) Methods of Numerical Mathematics. Springer. Berlín, Alemania. 532 pp.        [ Links ]

6. Ringler TD, Randall AD (2002) A Potential Enstrophy and Energy Conserving Numerical Scheme for Solution of the Shallow-Water Equations on a Geodesic Grid. Month. Weather Rev. 130: 1397-1410.        [ Links ]

7. Sadourny R (1975) The Dynamics of Finite-Difference Models of the Shallow-Water Equations, J. Atm. Sci. 32: 680-689.        [ Links ]

8. Salmon R (2004) Poisson-Bracket Approach to the Construction of Energy- and Potential-Enstrophy-Conserving Algorithms for the Shallow-Water Equations, J. Atm. Sci. 61: 2016-2036.        [ Links ]

9. Samarskii AA, Popov YuP (1969) Completely Conservative Difference Schemes. Zhurnal Vysshei Matematiki Matematicheskoi Fisiki 9: 953-958 (en ruso).        [ Links ]

10. Shokin YuI (1988) Completely Conservative Difference Schemes. En de Vahl Devis G, Fletcher C (Eds.) Computational Fluid Dynamics. Elsevier. Amsterdam, Holanda. pp. 135-155.        [ Links ]

11. Skiba YuN (1995) Finite-Difference Mass and Total Energy Conserving Schemes for Shallow-Water Equations. Russ. Meteorol. Hydrol. 2: 35-43.        [ Links ]

12. Skiba YuN, Filatov DM (2005) Finite Difference Splitting-based Schemes for Shallow-Water Flows Conserving the Mass and Total Energy. Memoria XI Congreso Latinoamericano e Ibérico de Meteorología y el XIV Congreso Mexicano de Meteorología. Cancún, México. Resumen r092.doc.        [ Links ]

13. Takano K, Wurtele MG (1982) A Four-Order Energy and Potential Enstrophy Conserving Difference Scheme. US Air Force Geophysics Laboratory Tech. Rep. AFGL-TR-82-0205. 85 pp.        [ Links ]

14. Vreugdenhil CB (1994) Numerical Methods for Shallow-Water Flow. Kluwer. Dordrecht, Holanda. 261 pp.        [ Links ]

15. White AA (2002) A View of the Equations of Meteorological Dynamics and Various Approximations. En Norbury J, Roulstone I (Eds.) Large-Scale Atmospheric Dynamics. Vol. I. Analytical Methods and Numerical Models Cambridge University Press. Cambridge, RU. pp. 1-100.        [ Links ]

16. Yanenko NN (1971) The Method of Fractional Steps. Springer. Berlín, Alemania. 429 pp.        [ Links ] ]]> 1 1981 109 145-171 2 1970 352 3 1947 43 50-67 4 1984 6 251-261 5 1982 532 6 2002 130 1397-1410 7 1975 32 680-689 8 2004 61 2016-2036 9 1969 9 953-958 10 1988 135-155 11 1995 2 35-43 12 2005 13 1982 85 14 1994 261 15 2002 1-100 16 1971 429