Document
-
Rating
-
Date
September 2018 -
Size
16.6MB -
Views
2,055 -
Categories
Transcript
- 25 - EL METODO ESPECTRAL PARA LA RESOLUCION DE ECUACIONES DIFERENCIALES EN DERIVADAS PARCIAl ES. SU APLICAC ION A LA ECUAC ION BAROTROPI CA DE LA VORTICIDAD Jos é Antonio Garcia-Moya Zapata. Instituto Nacional de Meteorología Abstract The Spectral Method is showned as one of the better methods to resolve P artial Diferencial Equation Systems. T his is because it improves the n umeric solutio n with respect to analytic solution. Moreover, the numeric solutio n mant ains exactly the integral invariants of the analytic solution. In this paper, it is explained the m atematical underlying principies of the method. At the end, it is applied to the reso lution of the Barotropic Vorticity Equatio n fo r both local and spheric coordinares. Resumen E n la resolución de sistemas de ecuaciones diferenciales en derivadas parciales, el método espectral se ha mostrado como uno de los más eficaces. Esto no sólo es debido a que mejora las soluciones numéricas obtenidas con respecto a la solución analítica, sin o también a que la solución numérica conserva exactamente los invariantes integrales de la solución analítica. Se explica en este trabajo detalladamente el fundamento matem ático del método. Después, y como aplicación, se resuelve la ecuación barotrópica de la verticidad tanto en coordenadas locales como esféricas. l. 1nt roducción Desde los primeros experimentos en predicción numérica del tiempo, a finales de los años cuarenta, los métodos de puntos de rejilla han Revista de Mc:teorología, A.~. E. ·Diciembre 1985 sido los más usados en los modelos numéricos del flujo atmosférico a escala planetaria. Sin embargo, se han empezado a usar métodos alternativos. U no de estos métodos es el espectral. E l método espectral fue in troducido en los modelos meteorológicos por Silberman ( 19 54), que estudió la ecuación de la verticidad para una atmósfera barotróp ica en coo rdenadas esféricas. Estudios posteriores demostraron algunas propiedades importantes del método. El trabajo de Ellsaesser (1966) demost ró que para un modelo barotrópico equilibrado el método espectral podría competir, en cuanto a ejecución y eficacia, con el de puntos de rejilla usado por el C.S. National Meteorological Center. Sin embargo, aunque el método espectral fu ncionaba bien para modelos simples en baja resolución, no se le consideraba como una alternativa real a los métodos de puntos de rejilla para integraciones de modelos no-adiabáticos más complicados. La razón de esto era que el cálculo de los llamados " coeficientes de interacción" provenientes de los términos no lienales de las ecuaciones necesitaba dem asiado tiempo de ordenador cuando se quería conservar un n úmero de compone ntes ond ulatorias en los desarrollos en serie, suficientes para obten er una adecuada resolución espacial, que incorporara procesos físicos locales, como por ejemplo, el ajuste con vectivo. La situación cambió comp letamente con la introducción del método de la transformada, desarrollado independientemente por Eliasen y otros (1970) y Orszag (1970). En este método no están implicados los coeficientes de interacción -26 - y la memoria de ordenador y las operaciones aritméticas requeridas se reducen sustancialmente. Además el método implica un estado en cada paso de tiempo, en el que lo valores puntuales de las variables se calculan en una rejilla auxiliar en el espacio físico . Como han señalado Eliasen y otros (1970) ésta hace que se puedan incluir en el modelo los efectos no-adiabáticos locales, de manera similar a lo que se hace en los modelos de puntos de rejilla. como tales, pueden ser consideradas como elementos de algún Espacio Vectorial H (por ejemplo, el de todas las funciones continuamente diferenciables ). A continuación se describen brevemente algunas generalidades sobre el método espectral y las propiedades más importantes de los armónicos esféricos y, en particular, de los que utilizan las funciones asociadas de Legendre. donde g* es el complejo conjugado de g. También podemos definir una norma como: D espués se aplica el m étodo para resolver la ecuación barotrópica de la vo rticidad, tanto en coordenadas locales (tomando como funciones básicas las de Fourier) como en coordenadas esféricas (tomando los armónicos esféricos con las funciones asociadas de Legendre). 11. Si " f' y " g" son dos funciones de H podemos definir un producto escalar de la forma: (f, g) = Lf· g*dx Con este producto escalar y esta norma, H se convierte en un E spacio de Hilbert. Si suponemos ahora que em(x) es una base ortonormal de H (o sea, ( em, en) = O para n =1- m y (en, en) = 1) cualquier elemento F de H puede expresarse como: Armónicos esféricos. Funciones asociadas de Legendre Las ecuaciones en predicción numenca del tiempo pueden escribirse, de forma general, como: Sin embargo, en general esto no es una base ya que M no es finito. Los pmson las proyecciones ortogonales de F sobre el subespacio generado por em. lo que es lo mismo, o j = 1.. . J [II. l ] donde J es el número de variables independientes que aparecen en las ecuaciones y que son función de las tres coordenadas espaciales "x" y del tiempo "t", o sea, F¡ = F¡ (x, t). Para todas las aplicaciones meteorológicas se supon e que F¡ es tan continua, o sea, diferenciable, como se necesite. Los L¡ son operadores, generalmente no lineales, formados por derivadas parciales espaciales y, en algunos casos, integrales espaciales. Entonces, para cualquier instante de tiempo t, los campos F¡ son funciones continuas de "x" y, (II.2] Si suponemos el caso de una sola variable F = F(x, t) el sistema (II.l] se reduce a: ;;} t =L (F) con la condición inicial F(x, O) [II.3] = f(x). Ahora bien, como no podemos trabajar con un número infinito de componentes, tenemos que proyectar F sobre un subespacio finito de H, al que llamaremos H (procedimiento de truncación). Re.visrl. de. Mereorologfa, A.M.E. - Diciembre t985 j_ -27- R(F) el error cometido al truncar la serie infini- De este modo F se aproxima por: ta de F, o sea: M F(x, t) = 1: Fm(t) ern(x) [II.4] m=l R(F) El siguie nte paso es revolver la ecuación [II. 3] para F, o sea: ~F -----=- L(F) F(x, O) ;}t = F(x, O) Sin embargo, esto no es posible ya que F pertenece aH (el subespacio truncado y L(F) no tiene por qué. Entonces tenemos que reemplazar esta ecuación por: - ;;, F _ d onde L(F) es la expresión truncada de L(F) que ya sí pertenece a H. Derivando parcialmente con respecto al tiempo la ecuación [II.2] resulta: ~ t d pm d t f ~F - - e *m(x) dx S ~ t = LL(F)e*m(x) = dx O sea, que: De esta manera hemos obtenido un sistema de ecuaciones diferenciales ordinarias. Entonces, si podemos evaluar con algún procedimiento el segundo miembro, podremos conocer la evolución temporal de los coeficientes f m, y con ellos, la evolución temporal de la variable meteorológica F con la ecuación [II.4]. E ste tipo de truncaciones constituyen el "método G alerkin" en el que las funciones em(x) son las llam adas funciones test o funciones básicas. Como ya hemos dicho, L(F) no tiene por qué pertenecer a H por lo que podemos definir como Revista de Metcorologfa, A.M. E. . Oíciembre 1985 - L(F) - = L(F) - - L (F) =1= O Ahora bien, puesto que los em(x) son ortogonales, R(F) pertenece al complementario de H (o sea H - H ) y, por lo tanto es ortogonal a todos los elementos de H y, en particular, a F, o sea: ), R(F)F dx =O (II. 5] Esta propiedad es muy importante, como veremos después. Entonces la serie de coeficientes pmforma_n el "espectro" de F, y de ahí el nombre de " espectral" dado a este método. - =L(F) ~ t ---= ~F =~- Puesto que la mayor parte de las veces se trabaja en predicción n umérica con problemas hemisféricos o globales, es conveniente elegir como funciones para los desarrollos en serie las que tengan simetría esférica. Además, también es conveniente elegir funciones que simplifiquen el operador L de la ecuación (II.3]. Como hemos dicho antes, este operador incluye operadores diferenciales espaciales (derivadas verticales y horizontales, laplacianos horizontales, etc.), té rminos no lineales y, posiblemente, algunas integrales verticales. Esta simplificación se obtiene si las funciones elegidas son funciones propias de algunos suboperadores de L. Este es el caso de las funciones armónicos esféricos: donde, como siempre, Aes la longitud,n [II.14b.] - -29- [II.14c.) La propiedad [II.14b] implica que [II.12) puede escribirse como: co F= I [11. 15] m=-c:c n=-co La propiedad [II.14c] prueba la sim etría o antisimetría de las funciones del desarrollo respecto al Ecuador. Hay algunas fórmulas de recurrencia para calcular los valores de P~, la usada en el Centro E uropeo de Predicción a Plazo Medio (CEPPM) es: Ahora bien, sabemos que tenemos que truncar el desarrollo en serie [II.15] para tener un número finito de componentes. Antes, hay que tener en cuenta que cuando m y n aumentan es como si disminuyeran las escalas horizontales, o sea, aumenta la resolución espacial del modelo. De manera que usando los armónicos esféricos (Y~) tenemos un control directo sobre las escalas que queremos despreciar, de la misma manera que los métodos de puntos de rejilla lo hacen al elegir el intervalo de la rejilla. La forma del desarrollo truncado es: lvf F= I pm(11 ) = cmpm-2(11 ) - ¿mn p n-m-2(11 )+ n ,..... n n-2 r 1 r +e~ P~_ 1 (J.!) n= - M I n 1m 1 [II.16) El hecho de que m vaya de -M a M asegura que F es real, ya que de [11.13] y (II.14a] se deduce que F;;m = F~· , y po r lo tanto F;;m • y~m + siendo: e~ = [ 2n + 1 m + n-1 F~ m+n 2n + 3 n-]3 1/2 - m+ m+n 2n + 1 m+n - 1 2n- 1 m+n m+n-2 e~ = [ 2n + 1 . n - m] v z 2n- 1 n +m para m > O También hay otra fórmula de recurrencia para calcular las derivadas meridionales (ya que los armónicos esféricos no son funciones propias d · Y~ es real. Como se puede demostrar que F converge a F Garraud & Simmons 1983), la principal cuestión que queda por resolver es la elección de N(m). Las dos truncaciones más comunes usadas en los modelos de predicción numé rica son las llamadas "triangular" y "romboidal", que pueden representarse en el plano (m, n) como se muestra en la figura 1. n-m+ 1 de - N (m) - ): La razón para elegir una determinada truncación se relaciona con el mejor uso de los grados de libertad (es decir, de las com ponentes espectrales retenidas), o sea, que trata de obtener mejores resultados con el mismo número de grados de libertad. dj.! e;:1 =n con R~vista e~= [ P:;+l - (n + 1) nz- mz ] 4nz - 1 de Meteorologfa, A.M. E. · Diciembre 1985 112 ~ P:;_ 1 [Il17) Aunque al principio, los modelos espectrales usaban truncaciones romboidales, los más re cientes usan truncaciones triangulares, ya que parece que dan mejores resultados y son más eficientes desde el p unto de vista computacional. - 30 - n n X X X X X X X X X X X X X X X X X X X X X. X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X X m Truncación triangular Truncación romboidal Figura l.-Representación en el plano (m, n) de las trun caciones romboidal y triangular. 111. l a ecuación barotrópica de la verticidad en coordenadas cartesianas o locales y esféricas En este punto nos vamos a limitar a obtener las expresiones más adecuadas de la ecuación barotrópica de la vorticidad para la aplicación del método espectral, para ello partimos de la conocida expresión de dicha ecuación en el nivel de no divergencia: - d -(/;+ f) dt =o Desarrollando el operador dj d t como ;}j ;} t + y teniendo en cuenta que el parámetro de Coriolis "f' sólo depende de la latitud
, A y radio de la Tierra "a") tenemos que transformar el segundo miembro de la ecuación y el operador laplaciano. Para ello, tenemos que tener en cuenta las relaciones siguientes: +1 --2-"'-~-22-J - !l ""/\. Una de las propiedades más importantes que cumple la solución analítica de esta ecuación es la de conservar tres magnitudes m edias, éstas son, el momento angular medio (M), la energía cinética medía (K ) y la vorticidad cuadrática media o "enstrofía" (E). O sea, que la solución analítica cumple: dM dk = o; dt dE =0· - - = 0 ' dt dt [III.S] donde las barras indican un valor medio global definido por: f f¡ 2 ( 5 =1- " 41t o ( )d'Ad¡.t = -1 [III.6] = - 1- dx = a · cos q> · d'A ( ( ) ds Js S dy = a · dq> y donde M = u · a · cos q>, K = (V 2)/ 2 y E = = 2)/ 2. Para proba r la invariancia del momento angular medio, deduciremos para él una expresión alternativa. (s Y el operador laplaciano resulta: 1 + - -- - - cos2 q> 1 M= - ;;)2 g A2 Si hacemos ahora el cambio ¡.t = sen q> y agrupamos términos, resulta la ecuación en coordenadas esféricas y el operador laplaciano de la siguiente manera: v 1 - - v 2\jf- - g t a2 2 \ji - ~').._ ___:_y_ • gcp _a fu S u a cos q> ds = S f( u~-V f az ~!l ~ L S ds = S ;;)J..l. - az ) ds= acoscpg/.... 1 acosg/.... cosq> [III. 3] ds =a2S J..l. S ;:}V ;;¡').._ ~q> S f ¡.t S ~(¡.tcosq>) a ~q> Sds S Y, por tanto: (1- ¡.t2) _g_] ;:}¡.t + [III.4] Revisca de. ;\l[ereorologfa. .\ ..M.E . . Diciembre 1985 dM d t = S f S g s J..l. - - ·ds d t [III.7] - 32- 1 [ -(k 2 Teniendo en cuenta la expresión de la ecuación no divergente de la verticidad resulta: dM dt = ~ ( S J, = !-! [- v '\7 es + f) J ds = [ ~ ( [~ V(
+ f) = V (S+ f) V - (<¡> + f) · · V ·V, pero como se trata de f1ujo no divergente el segundo sumando es cero. Entonces podemos transformar el integrando de la siguente manera:
- - v [~Cs + f)
Z1
A
vz) = lvl lv
\ji
'JI
1
1 .
- d- 1v 1 dt
=
-
g
v-d t
-33-
d t
de manera que el teorema de Gauss implica:
:;, K
e \ji J,
-1
= -
S
:;,t
____:__l.
ds
[III.8]
:;,r
De la misma forma, obtenemos que:
:;, E
= _1 e 1; ____:__l. ds S
:;, t
J,
[III. 9]
dt
U sando la ecuación de la vorticidad, los integrandos de ambas magnitudes pueden escribirse como:
=-\jf~~ + DV + \jf~+ D~ ·V =
= - 9'[\jf(s + D · vJ + Cs + D V · V "' = = - ~ [\jf(s + DVJ + Cs + D(f V\jf) · "' = = - v[\jf[s + D· V]; = Cs + D
- f
di;
IV.
Ecuación de la verticidad en coordenadas cartesianas. Desarrollo de la función de corriente en serie de Fourier
Partimos, en este caso, de la ecuación [III. 1] que es la ecuación barotrópica de la vorticidad en forma vectorial y en coordenadas cartesianas:
-
d
~
-
d t
d t
d
2n
:;,¡;
~
k
d t
1;
(/; + D2 Y' V - f - d t
=-
~
[+
Y'
Cs + D2 • v ] - 2 n 11
d t
Y aplicando el teorema de Gauss a las integrales resulta que:
dR dt
=Oy -
dE
- =0.
dt
Revtsta de Meteorología, A.1LE. · Diciembre 1985
2n -
y
+-
= \jf(x,
y, t,)
\jf(x+ - - ,
--=-Y'
d
d
\ji X
~~
Y' (/; + D- f - - = - V \l
-
~
( Y' 2\ji)+ (kA \l \ji). Y' (\7 2 \ji) + ~.-=o
Suponemos que la función de corrientes es periódica en las coordenadas "x" e "y", o sea:
d t
d t
~~
= - (/; + D V
Estas tres invariancias son muy importantes ya que constituyen lo único que sabemos de la solución analítica de la ecuación baratrópica de la vorticidad, por tanto, cuando obtengamos una solución espectral aproximada (truncada) lo primero que deberemos hacer es exigirle que cumpla estas tres invariancias, para así estar seguro de que la solución truncada tenderá a la analítica cuando el número de términos retenidos tienda a infinito y que, por tanto, se la pueda considerar una solución aproximada de la ecuación, consistente con la solución analítica.
1
'
t) =
Entonces para este caso del plano-~ y esta condición de periodicidad, las funciones ortogonales básicas apropiadas son las de Fourier, o sea: q>m, n
(x, y)=
e i(mla+nly)
Que, como se sabe, son las soluciones propias de la ecuación de Poisson:
-
Siendo b = (m 2 corriente es: \!f(x, y, t)
•
k2
+ n 2 • F). Y la función de
~ ~ ~Cm,
n (t) ei(mkx +nlrl
[IV. 1]
34 -
[, o sea, m · k, donde k es ahora el número de onda según el eje "X".
El producto
(k
\l \!') • \l ( v 2 \!') queda:
A
Además, para que \!f sea real los coeficientes del desarrollo han de satisfacer la cond ición:
Para simplificar la notación, introducimos el vector número de onda M = m · k · T+ n · e ·) y el radio vector R = x · T+ y · y entonces queda:
T
\!'(X, y, t)
= L
b=M· M
cm (t)eiM · R
[IV.2)
Ya que + Lr H x
[k
A HJ
· L = k( H AL)=
Ahora, multiplicamos la ecuación de la vorticidad por el conjugado de las funciones básicas exp (-i · M · R), integramos sobre el dominio periódico y obligamos a que el error sea o rtogonal con las funciones básicas y resulta:
Con este desarrollo en serie la expresión de la vorticid ad es:
-
-
(L · L)
C M (t) , \]2 , iMR =
+ i~ I
dCl
e'( L- ~!) R +
dt
L,CLei(L-M)R+
L
e H V é iR =
L iHCH eiHR [IV.4)
H
Sustituyendo en la ecuación de la vorticidad queda:
(:L · L)
L
eH eiHR ] . [
dCL dt
e'LR +
[I
i(k A
[ "!K [
e-'M R d y dx
H
[IV.S)
-I
L
· dydx =
Y además resulta que:
= L
I I(L · L) H
~(HA L) CHCLei(H+L - M)TJ.
[IV .3)
\] '!'
11
en
=O
Esto se hace, lógicamente para cada M del desarrollo en serie incial de \!f. Como las funciones básicas de Fourier son o rtogonales, cada integral se anulará salvo que el exponente de la exponencial sea O, o sea, que para cada M queda:
I
- CM · M)
H
H)
H
-~CL(L. L) i E e'LR) ] + i ~ LL L e L eiLR=eN
Teniendo en cuenta que H A(M - H) = H A la ecuación diferencial ordinaria de cada coeficiente CM queda:
M - A A H,
X
donde k es el vector unitario según el eje vertical, eN es el error cometido al sustituir la función de corriente por el desarrollo en serie y L, es la componente " x" del vector número de onda
Revista de Mereo rología, A.M. E. -Dic iembre 1985
-Lx HY+
dCM
i ~m k CM
dt
M· M
- - = ---,----- +
I H
-35-
Entonces, si N es el número de términos del desarrollo de \jf, obtenemos así un sistema de N ecuaciones diferenciales ordinarias para los coeficientes del desarrollo. El último término de la ecuación surge de la interacción entre diferentes ondas que provienen del término de advecciónno lineal de la ecuación de la vorticidad. En particular, la onda M está afectada por la interacción de todas las ondas H con la M - H. Cuando se desprecia este término de interacción resulta un sistema de ecuaciones lineales desacopladas que tiene como solución la onda de Rossby.
nales, vamos a demostrar que en este caso se siguen conservando con este esquema la energía cinética media y la verticidad cuadrática media o ''"enstrofía". La energía cinética media para la región es:
-
K=
r 27t/k 21t/l
~
-
ro -V·V - - dy dx=
o
2
(k A V\ji) (k A V \ji) 2
Pmk
Y a que, como el v \ji es un vector horizontal, el ángulo que form a con el vector vertical k es 90° , luego lk A V\j/1 2 ~ [ lk 1 · 1V\jf 1· sen 90°)2 = 1'\7 \j/1 2 = V\jl. Sustituyendo el valor de V\jl d ado por [IV.4), resulta:
'v\v ·
V elocidad de fase: _ _.......:....P_m_k_ __ mzkz + n ZF
-
r
r 27t/k .2 27t;l
K = o
H
o
Como vemos, la velocidad de fase de la onda resultante es efectivamente la de la onda de Rossby. Entonces cada una de las ecuaciones del sistema [IV.6) puede resolverse con un esquema de diferencias finitas centradas. Si despreciamos los términos de orden superior al primero resulta:
c~~r = c~-~r + 2~t •
[
d~M J r
[IV. 7)
1
-2 rr M
D e la condición de ortogonalidad de las funciones básicas deducimos que la integral será cero, salvo cuando A = - ~f, y e ntonces quedará:
1 - K= - 2 L M . M CM c_M = M . [IV.8]
Enton ces con un proceso iterativo para las ecuaciones [IV. 6) y [IV. 7) puede calcularse la evolución temporal de los coeficientes del desarrollo de la función de corrientes y sustituyendo estos valores en las ecuaciones [IV .2), [IV.3), [IV.4) y [IV.S] pueden calcularse los valores previstos de la vorticidad y la v elocidad del viento. Además, puesto que el procedimiento de Galerkin conserva los invariantes que surgen de las no linealidades cuadráticas de las ecuaciones origi-
Revista de Meteorologfa, A. .M.E. - Diciembre 1985
- 21
1: M
IMI 2 le
M
12
Ya que, según vimos antes, CM = c~M" y ésta es la expresió n de la energía cinética para todo el dominio. Para el caso de la "enstrofía" resulta:
[21tjk [21t¡1 E =
Jo
t /;
2
dy dx
=
-
36 -
Con lo que con las funciones básicas utilizadas y la aproximación de Galerkin, la condición [IV. 1O] se transforma en:
1: i~M M
=:E:E
1C
12
M
X
+ 1: 1: - M H k (H
A
H M
í27t/k í27t/l L Jo
e•