Un método de elementos finitos incondicionalmente estable en norma uniforme para resolver las ecuaciones de Euler 2D

Revista Internacional de Métodos Numéricos para Cálculo y Diseño en Ingeniería. Vol. 9, 3, 231-258( 1993) UN METODO DE ELEMENTOS FINITOS INCONDICIONALMENTE ESTABLE EN NORMA UNIFORME PARA RESOLVER LAS ECUACIONES DE EULER 2D TOMAS CHACON REBOLLO e IBRAHIM BLESS RANERO* Departamento de Análisis Matemático, Universidad de Sevilla RESUMEN Este trabajo presenta dos algoritmos de tipo transporte e interpolación con Elementos Finitos para la resolución numérica de las ecuaciones de Euler para flujos incompresibles bidimensionales en todo el espacio IR'. En la primera versión de nuestro algoritmo, la vorticidad es discretizada mediante elementos finitos triangulares de primer grado. En la segunda, mediante elementos finitos triangulares de segundo grado. La velocidad se obtiene en ambos casos calculando exactamente el producto de convolución del núcleo de Biot-Savart con una aproximación lineal a trozos sobre cada triángulo de la vorticidad discreta. Se prueba que el primer algoritmo es incondicionalmente estable y convergente con una precisión de primer orden, en norma uniforme. Sin embargo, en la práctica la precisión alcanzada resulta escasa, debido a la difusión numérica introducida en el paso de interpolación. En el caso del segundo algoritmo, los ensayos numéricos muestran un notable incremento de la precisión, incluso para tiempos largos. Sin embargo, en este caso el algoritmo deja de ser estable en norma uniforme. SUMMARY We introduce two Finite Element transport-interpolation algorithms to solve the twodimensional Euler equations in the wole R2. In the first of these algorithms, the vorticity is discretized with triangular finite elements of degree one, and of degree two in the second one. The velocity is computed by convolution of the Biot-Savart kernel with a piecewise affine interpolate of the vorticity. We prove that the first algorithm is unconditionally stable in uniform norm, with first order accuracy. However, in practice its precision is rather low, due to the numerical diffusion introduced in the interpolation step. The second algorithm is shown numerically to produce a remarkable increase of precision, even for long integration times. However, in this case the algorithm is no longer stable in uniform norm. * Investigación financiada parcialmente por Proyecto DGICYT 9B91-0619 Recibido: Octubre 1991 OUniversitat Politecnica de Catalunya (España) ISSN 0213-1315 232 T. CHACON REBOLLO E 1. BLESS RANERO INTRODUCCION En este trabajo, nos interesaremos por la resolución numérica de fluidos bidimensionales incompresibles y no viscosos. Estos fluidos están gobernados por las ecuaciones de Euler, que se pueden escribir como una ecuación de convección pura, para la vorticidad del flujo. Los problemas de convección no lineal aparecen frecuentemente en ingeniería. Por ello, resulta de interés el desarrollo de algoritmos de resolución numérica de los mismos. La ecuación de convección se reduce a una familia de Ecuaciones Diferenciales Ordinarias a lo largo de las curvas características del fluido. Esto ha dado origen a los algoritmos "Lagrangianos" para las ecuaciones de Euler 2D, basados en la idea de calcular la vorticidad a lo largo de las curvas características del fluido. El primero de estos métodos fue el de "punto-vórtice", introducido por Rosenhead [18], que se basa en una discretización de la vorticidad como suma de masas de Dirac y tiene como principal ventaja, la de ser no disipativo, aunque no es estable para largos períodos de tiempo (cf.g). Al principio de los años 70, Chorin7, Kuwahara y Takarni introdujeron la idea de discretizar la vorticidad, sustituyendo las masas de Dirac por funciones regulares que las aproximen, para estabilizar el método. Esto dio origen a métodos con mayor precisión, no disipativos y estables, que han sido estudiados ampliamente en los pasados 15 años (Cf.13~14~15*3,8~17). Se trata de los métodos llamados "de Burbuja -Vórtice7'. Por otra parte, también es posible construir algoritmos lagrangianos con Elementos Finitos. Bardos, Bercovier y Pironneau introdujeron en [2] un algoritmo para resolver las ecuaciones de Euler formuladas en términos de la función de corriente asociada a la vorticidad, basándose en una discretización constante a trozos de la vorticidad sobre una triangulación. En este algoritmo, la función de corriente se aproxima mediante elementos finitos conformes de primer grado. De esta forma, la velocidad discreta es constante a trozos, pero la componente normal es continua a través de los lados de la triangulación. Este método posee la ventaja fundamental de los métodos de vórtices, al ser no disipativo. Además, es uniformemente convergente y permite la manipulación de condiciones de contorno. Sin embargo, el cálculo de las curvas características, presenta algunas dificultades computacionales, que se deben al hecho de que la velocidad no es globalmente continua. Un algoritmo puramente lagrangiano con Elementos Finitos fue a continuación introducido en [6] por Chacón y Hou. Este algoritmo está basado en una discretización afín a trozos de la vorticidad, calculando la velocidad directamente, mediante la convolución de la vorticidad discreta con el núcleo de Biot-Savart y trasladando los vértices de la triangulación a lo largo de las curvas características del fluido, en cada paso de tiempo. Este algoritmo es uniformemente estable, para un cierto período de tiempo, si las curvas características son calculadas con un algoritmo de segundo orden. Sin embargo, para obtener la estabilidad en largos períodos de tiempo, es necesario el uso de técnicas de remallado o un paso de tiempo muy pequeño. Nuestro propósito en este trabajo es el de introducir una versión de tipo "Vortex- METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 233 in-Cell" del algoritmo Chacón-Hou citado arriba. La malla móvil utilizada en éste para discretizar la vorticidad es reemplazada por una malla fija, siendo necesario interpolar la vorticidad tras cada paso de tiempo. Esto produce una notable mejora de las propiedades de estabilidad del método, producida fundamentalmente por la introducción del paso de interpolación. En la sección titulada Descripción de los Algoritmos introduciremos un primer algoritmo que utiliza esta técnica. La velocidad es calculada de la misma manera que en [6], mediante la convolución de la vorticidad discreta con el núcleo de Biot-Savart. Tambíen se introduce una variante de este algoritmo, en la cual la velocidad es calculada mediante la convolución del núcleo de Biot-Savart con un interpolado constante a trozos de la vorticidad. En la sección tilulada Análisis de Convergencia analizaremos las propiedades de convergencia y estabilidad de nuestros algoritmos. Probaremos que bajo suposiciones razonables para el operador de interpolación, ambas versiones del primer algoritmo son incondicionalmente estables y convergentes, con exactitud de primer orden en norma uniforme, para todo intervalo de tiempo finito. La sección titulada Operadores de Interpolación está dedicada a describir algunos operadores de interpolación que cumplen las hipótesis requeridas para asegurar la convergencia. Probaremos que las interpolaciones "puntual" y "promediada" son buenas elecciones para este propósito. En la sección titulada Propiedades de Conservación analizaremos las propiedades de conservación de nuestros algoritmos. El paso de interpolación que introducimos para calcular la vorticidad, produce un cierto aumento de la difusión numérica, debido a que los algoritmos de transporte e interpolación que conocemos no conservan el área. Daremos una demostración directa de que nuestros algoritmos conservan el área, con un error de segundo orden. Esta demostración se basa en el hecho de que las velocidades discretas son exactamente de divergencia nula. En la sección titulada Ejemplos Numéricos, mostraremos algunos ensayos numéricos en un problema con solución analítica conocida. Este ejemplo muestra una buena correspondencia entre las predicciones teóricas y los resultados numéricos. Sin embargo, existe un alto nivel de difusión numérica. Por último, en la sección titulada Una versión de Segundo Orden describiremos la segunda versión de nuestro algoritmo. En ella, la vorticidad se aproxima mediante elementos finitos triangulares de segundo grado. Las curvas características se discretizan hacia atrás en tiempo mediante un método no estándar de segundo orden. La velocidad se calcula convolucionando el núcleo de Biot-Savart con una aproximación afín a trozos de la vorticidad. El incremento de cálculo requerido se ve compensado por una técnica de cálculo rápido de la velocidad, que es la etapa más costosa de nuestros algoritmos. Mostraremos finalmente un ejemplo numérico en el que se aprecia una convergencia de segundo orden, manteniéndose el error prácticamente constante incluso para tiempos largos. 234 T. CHACON REBOLLO E 1. BLESS RANERO PLANTEAMIENTO DEL PROBLEMA Nuestro propósito es resolver numéricamente las ecuaciones de Euler 2D para flujos incompresibles en todo el espacio R2 u,t + u-Vu + Vp = 0, V-u = O enlR2x]0,T[, u(x, 0) = uo(x) en IEt2, lim u(x,t) = 0. IxI-)~ 1 Aquí, u(x,t) y p(x,t) representan el campo de velocidad y la presión del fluido respectivamente, en el punto x E IR2 y en el instante t, uO(x) es un campo de velocidad inicial dado. Además, [O,T] es el intervalo de tiempo durante el que analizamos el comportamiento del fluido y 1 1 denota la norma b" en IR2. Es conocido que el problema anterior es equivalente a la formulación "velocidadvorticidad" de las ecuaciones de Euler con ausencia de frontera finitas: La función K(x, t) anterior es conocida como el núcleo 2D de Biot-Savart. Este núcleo tiene una singularidad en el origen, aunque es localmente integrable. Más aún, la convolución con K es un operador acotado de L~~,(IR~) en ~l!p(R~), 1 5 p < +m. En lo que sigue, utilizaremos ampliamente esta propiedad. Se conocen resultados de existencia y regularidad de las soluciones del problema 1, por ejemplo para wo regular12. Consideraremos el caso wo E C2(lR2) con soporte compacto. Entonces existe una solución w E c2(BI2 x [O, TI), para cualquier T > 0. Además, esta solución tiene soporte compacto en un tiempo cualquiera, es decir, existe una constante RT > O tal que Observemos, que si se da el campo de velocidad u, entonces la ecuación de transporte para w en (1) puede ser integrada explícitamente. Para cada S E [O, T] fijado, consideremos la ecuación de las curvas características X(t; S, x) asociadas al campo de velocidad u con origen en S: dX -(t; S, x) = u(X(t; S, x), t) para t E [O,T], X(s; S, x) = x. dt (3) entonces w es constante a lo largo de la curva t E [O,T] -t (X(t;s,x),t) E IR2: w(X(t; S, x), t) = W(X, S). (4) METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 235 En particular, resulta que Utilizaremos este hecho en lo que sigue para construir nuestros algoritmos. DESCRIPCION DE LOS ALGORITMOS Sea una triangulación de Bt2, donde h es el diámetro máximo de los triángulos de Th. Denotemos por {a;); 2 1 los vértices de los triángulos de Th. Definamos el espacio de Elementos Finitos lineales a trozos donde Pk, para k 2 O entero, es el espacio de los polinomios sobre Bt2 de grado 5 k. Una función vh E Vh, está únicamente determinada por los valores u; = vh(cri),Vi > 1 (C f.5). Consideremos un operador de interpolación lineal Supondremos que existe X 2 O independiente de h, tal que si sop(v) c B(0, R), entonces sop(vh) c B(0, R + Xh). (6) Definamos ahora dos versiones de un primer algoritmo para discretizar las ecuaciones de Euler, basados en la descripción lagrangiana (4). Consideremos N 2 1, y At = TIN. Llamemos t, = nAt, O < n < N, y denotemos por GK a una aproximación de w(., t,). Algoritmo Al 1. Inicialización: 2. Dada 6; E Vh con soporte compacto, definimos (a) La velocidad discreta üh por (b) La característica discreta í?; retrocediendo con el método de Euler 236 T. CHACON REBOLLO E 1. BLESS RANERO (c) La vorticidad discreta 6;+' en el instante tn+i mediante transporte $ interpolación, Observaciones 1. e es obtenida discretizando (3) hacia atrás en tiempo mediante el método de Euler, con s = tn+i. 2. El cálculo de la velocidad discreta mediante (7) se puede hacer de la siguiente manera: Denotemos por {q;) la base canónica de Vh7 dada por q;(aj) = fijj. Dado un triángulo T E Th7 denotemos por I(T) el conjunto de índices i E IN tales que a; es un vértice de T. Como las únicas integrales a calcular son - x')qi(x')dxl, para i E I(T) estas integrales pueden ser expresadas analíticamente como funciones de x y programadas directamente. Esta definición es consistente: Sea R, = in f {R >_ Olsop(Wz) c B(O, R)). A partir de la definición de la velocidad discreta (7) tenemos entonces If'Z(x)l 2 1x1 - At RnlGfilmDe aquí que sop(6; O F;) c B(0, RE), donde RE = Rn(l + AtlGhloo) METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 237 R,+i 5 Ri + Xh, a partir de la hipótesis hecha sobre rh. De esta manera, zú;+' tiene soporte compacto y fiX+' puede ser correctamente definida, ya que K E [L/,,(R~)]~. La segunda versión de nuestro primer algoritmo permite hacer un cálculo más rápido de la velocidad discreta sin pérdida de orden de convergencia, como veremos. Consideremos el espacio de las funciones constantes a trozos sobre 7: Hh = {vh : Dt2 + IR, vhl~ E PO,VT E z} Consideremos también un operador de interpolación lineal Sh : c0(Dt2) + Hh. Supondremos que sh verifica las dos propiedades siguientes: si sop(v) c B(0, R), entonces sop(shv) c B(0, R + Xh), (11) para alguna constante X > O; y si v E L~(R~), entonces shv E L~(R~). (12) ahora podemos describir la segunda versión de nuestro primer algoritmo: Algoritmo A2 Todo como en el Algoritmo Al, con (7) reemplazado por Esta definición es consistente: A partir de las propiedades (11) y (12), obtenemos C entonces, como en el Algoritmo Al, esto implica que 27th" tiene soporte compacto. Observación Para calcular la velocidad mediante (13), la única integral que se necesita es que se puede obtener analíticamente, de la misma forma que en el caso de la interpolación constante a trozos. 238 T. CHACON REBOLLO E 1. BLESS RANERO ANALISIS DE CONVERGENCIA Probaremos ahora la convergencia uniforme de los algoritmos Al y A2 a la solución de (1)) bajo algunas suposiciones para los operadores de interpolación rh y sh. TEOREMA 1 Supongamos que vh verifica las siguientes propiedades: si v E L~(R~), entonces rh v E L"(R~) ylrhvl, 5 Ivloo, (14) donde r; es una constante numérica. Supongamos también que h = OAt, donde O es una constante numérica positiva del orden de la unidad. Entonces, el Algoritmo Al converge uniformemente a la solución de (1). Además, la discretización es de primer orden y se cumple la siguiente estimación para el error de discretización: donde C es una constante que depende únicamente de T y wo. Demostración Analizaremos separadamente la estabilidad y la consistencia. Para simplificar la notación, definamos - n+l Observemos que la sucesión {wh Inzl es obtenida cuando un paso del Algoritmo Al es aplicado a la solución exacta de (1). Estabilidad: A partir de la propiedad (15) se tiene que ['lb;" - '6;+'lm = Irh(wn 0 Pn - '6; 0Ft)lm 5 lwn 0Tn - u; 0F2loo Entonces ~q+l - ti$+llm 5 Iwn oTn - w; oFJm + l(wn - 'LO?) Oel, Además, a partir de la propiedad (14) tenemos METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 239 17U;+lloo 5 12o; 0 F,yoo 5 IG;loo < l~:loo < Iwolw Entonces, de (10) se tiene que Rn+l I Rn(1 + IwolooAt) + Xh. y de aquí se obtiene que (4 max Rn _< RT ; independientemente de At, O<n<N donde Ahora bien, de (9) sigue que si tomamos, por ejemplo, At < 1, 19(x)l > 1x1 - ~T)/~olrn~t (4 > 1x1 - RT lwoloo. Entonces, donde RT viene dado por (2). Además, como en (8) obtenemos Iu(.,t)loo < RTlwoloo; o 5 t < T. De esta manera Irn(x)l > 1x1 - RTlwolooAt. Consiguientemente, sop(wn o Tn) c B(0, RZT), O 5 tn < T, con RZT = RT(~ + lwoloo) Llamemos ahora RJT = max(RIT, R2~) y definamos BT = B(O,RJT), QT = BT X [O,T]. Entonces 246 T. CHACON REBOLLO E 1. BLESS RANERO Por tanto, bastará con encontrar cotas uniformes en Lp para Vüh. 1. Recordemos que K es un operador acotado de L~~,(IR~) en w'J)(~R~), 1 < p < +m. Entonces, en el caso del algoritmo A2, donde la constante cp puede depender únicamente de T y wo. A partir de (ll), (12), (24) y (27), obtenemos IVühl, < cplshWhloo < ~~lchl~ < C~IWOI~. (31) En el caso del algoritmo Al, la segunda desigualdad en (31) se satisface directamente. 2. Sea RIT definido por (17) y definamos las constanes Demostraremos recursivamente que si O < T < To y At es suficientemente pequeño, entonces IVÜhloo < Cl(C0 + 1) Y IVchloo < Co + 1. (32) Consideremos primero que, a partir de (7) y del Lema 1, Si n = O, entonces (32) se cumple inmediatamente. Supongamos que (32) es cierto para O, 1,. . . , n. Demostraremos que lo sigue siendo para n + 1. Consideremos que Supongamos ahora que Entonces De aquí que, 5?~ sea globalmente inversible, y (f';)-' E c1(1EL2). Definamos METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 247 donde int(T) denota al interior topológico del triángulo T Entonces, aOT = (T;)-'(aT) y 1 1 dy 5 dy = 0. = LT 1 det [(~e)(?(y))] 1 1 - 21Viii1&(At)2 Como V(WAOT~)~~, E cO(R~), IV(Wi o f;)lw = max sup lV(6f o p;)(x)l. ,E7h sEf2~ Consideremos ahora que Entonces IV'.újkll, < (1 + 2ClAtlVGhl,)lVWfl,. (36) Esta estimación permite acotar IVWf lw durante un intervalo limitado de tiempo. En efecto, definamos la sucesión & = IVG;I,, paran = 0,1, ..., N. - Como (36) es también verdadera para O, 1,. . . , n, entonces la sucesión {(k}i=l está acotada término a término por la sucesión {~k}~=l definida por Esta sucesión se obtiene aplicando el método de Euler progresivo a la ecuación diferencial ('(t) = 2Cit2(t), ((0) = to. La solución de esta ecuación diferencial es 248 T. CHACON REBOLLO E 1. BLESS RANERO que es una función creciente no-negativa, que tiende a infinito en el instante t=- l >To. 2C1€0 - Podemos probar ahora que Para ello, observemos en primer lugar que ello es cierto para (0. Supongámoslo cierto para k = 0,1,. . . , n. Definamos Pk = [(tk) - tk; O 5 k 5 N. Entonces, de (37) y (38) obtenemos y de aquí que I~k+ll I [l + (2Co + i)~t]lpkl + 2~:~,3(~t)~, O 5 k 5 n. Con esta desigualdad y (34) llegamos a la cota deseada: Entonces obtenemos nuestro resultado: Observación El hecho de que la aplicación discreta de flujo sea invertible con inversa uniformemente acotada globalmente puede utilizarse para dar una demostración sencilla de la convergencia del algoritmo Al con rh en normas LP para cortos intervalos de tiempo. METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 249 EJEMPLOS NUMERICOS En esta sección, analizaremos los resultados de la ejecución práctica de los algoritmos anteriores en un caso test. La convergencia de primer orden predicha por el teorema 1 es confirmada numéricamente en velocidad y vorticidad. En nuestros tests utilizamos una solución estacionaria (u, w) de las ecuaciones de Euler 2D (1). Esta solución es con T = \/x: + x;, Las funciones u y w son muy regulares. El soporte de w es el círculo unidad. El flujo correspondiente es radialmente simétrico y rota alrededor del origen, produciendo un gran gradiente normal en la velocidad tangencia1 en una banda estrecha alrededor del círculo T = 0.4. Esta solución de las ecuaciones de Euler 2D ha sido ampliamente utilizada para mostrar las propiedades de convergencia de métodos de vórtices. De esta forma, tenemos un buen test para mostrar los resultados obtenidos por la ejecución de nuestro algoritmo. Para discretizar la vorticidad hemos utilizado triangulaciones uniformes, con una talla h = 0.15, h = 0.10 y h = 0.05. Si los triángulos se enumeran con cuidado, es posible encontrar el triángulo que contiene un punto x E 1Et2 dado en un número limitado de operaciones. Entonces, el paso de interpolación requerido por la actualización de la vorticidad resulta poco costoso computacionalmente. Para triangulaciones generales, este costo se ve incrementado, aunque puede ser mantenido en un número de operaciones del orden de si es programado con suficiente cuidado. La principal dificultad práctica que encontramos cuando programamos las dos primeras versiones de nuestro algoritmo es el cálculo del soporte de la vorticidad discreta. La estimación (10) puede ser utilizada para dar un círculo Cn+' que contenga sop(2ún+l), a partir de sop(2ún). Sin embargo, calcular 2ún+' en todo Cn+l requeriría una gran cantidad de trabajo computacional inútil. Para obviar esta dificultad, en nuestros cálculos hemos utilizado una técnica especial para "predecir" más precisamente sop(~Z"+~). Esta técnica debe ser descrita como sigue: (a) Inicialización: Aproximar l?: = sop(2ú0) mediante una línea quebrada de vértices {t!)E1. (b) Paso de tiempo: 250 T. CHACON REBOLLO E 1. BLESS RANERO i. Dado el conjunto de puntos {(r)gl que describe la frontera computacional I'g del sop(Wn), definimos los puntos (:';n+' = (r + A t üh((F) (Método de Euler progresivo) ; ii. Definimos la frontera aproximada I'z+' del sop(zZn+') como la línea n+l M quebrada de vértices {ti );=,, con las mismas conexiones de I'g. iii. Definimos el soporte computacional de Wn+' como la unión de todos los triángulos de que quedan dentro de I';+', más los que se intercepten con I';ltl. Aunque esta técnica introduce una cantidad adicional de error de discretización, nuestros cálculos muestran que éste sigue siendo de primer orden. Para medir los errores hemos utilizado una seminorma de tipo L2. Para una función v : R2 +- IR2, esta seminorma se define como sigue: donde al es el j-ésimo vértice del triángulo T. El error normalizado en velocidad es definido por, Con esta definición de errores, el orden de convergencia, por ejemplo en velocidad pu(t), ha sido estimado utilizando dos valores consecutivos de h: Aunque hemos demostrado la convergencia de primer orden en norma uniforme, por razones prácticas es más conveniente calcular los errores utilizando normas 12; esto permite obtener curvas de error más regulares. En cambio, las normas 1, producen grandes oscilaciones en las curvas que muestran los órdenes de convergencia, que desaparecen cuando utilizamos normas 12. Las Figuras 1 y 2 muestran las curvas de error en velocidad y vorticidad respectivamente, correspondientes al Algoritmo A2, con At = h. Se observa un crecimiento exponencial de los errores, aunque decrecen con h. En nuestro caso, este crecimiento es más rápido que en los métodos de burbuja vórtice, debido a la etapa de interpolación. Como la velocidad es más regular que la vorticidad, podemos esperar que los errores en velocidad sean más pequeños que en vorticidad. Sin embargo, vemos que sucede lo contrario. Este comportamiento aparentemente anormal se debe a pérdida de precisión cuando reinterpolamos la velocidad mediante una función constante a trozos, en el cálculo de la vorticidad con (13). Las Figuras 3 y 4 muestran la evolución de los órdenes de convergencia calculados en velocidad y vorticidad, respectivamente. La curva correspondiente a METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D EXTREMA l Y'S MIS I : am 0.9 1.m 1.50 2.- 2.50 iw 1se4.m Figura 1 C.30 4.2 EXTREMA l Y'S MIS 1 . Figura 2 T. CHACON REBOLLO E 1. BLESS RANERO UODULEF : chaco" EU1 12/10/88 raterfi.dbj NUMBER OF CURVES : 3 EXTREMA l X'S AXIS 1 : o. 4.2 EXTREMA l Y'S AXIS I : 0.83 2.6 - : REFERENCE .......... : h-.l5 vr hm.10 . ---. - - : h..10 va h-.O5 CURVES'S DRRWN m Figura, 3 MODULEF : chactin eul 12/10/88 1 CMRGEKE (ROER (vmrxcxrv) I ratea--0.dbj 2. 68 NUUBER OF CURVES : 3 2. 48 EXTREMA ( Y'S AXIC I : O. 4.2 2. 28 EXTREUA Y'S AXIS 1 : 0.82 1.2 2. m - : REFERENCE .......... 1. m : h'.lS "S hm.10 .-S -. . - : hm.10 vs h-.O5 1. 68 l. 48 l. 20 1. BB CURVES'S DRAWN .................................................. .................... 8. 88 8. 68 0.- LY 1.- 1.- 1.- 2.Y 03.Y4.W m Figura 4 METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 253 la velocidad muestra una exactitud de segundo orden en los momentos iniciales, debida probablemente a ciertas simetrías en las triangulaciones que producen cancelaciones. Se observa una transición suave al primer orden con el paso del tiempo. Finalmente esta exactitud se pierde progresivamente cuando el tiempo crece, como es usual en la resolución con Elementos Finitos de muchos problemas de evolución. Además, la curva correspondiente a la vorticidad muestra un comportamiento más estable en este caso. La estimación de orden en el tiempo t = O ha sido omitida, ya que nuestro error calculado es cero. Observemos que en ambos casos las curvas correspondientes a pequeños valores de h están dentro del orden teórico p = 1. La ejecución del Algoritmo Al en este caso test es similar, aunque todos los errores son algo más pequeños, en concordancia con la mayor exactitud en el cálculo de velocidades. 1 UNA VERSION DE SEGUNDO ORDEN Ahora nuestro propósito es describir una versión del algoritmo con segundo orden de consistencia. Si tenemos en cuenta que T = N At, para anular los efectos de la propagación de error al cabo de un tiempo T, convendrá exigirle a nuestro algoritmo un error de consistencia local de tercer orden en tiempo y en espacio. Para ello, es necesario en primer lugar construir un esquema con error de consistencia O(At3) para discretizar la ecuación de las curvas características (3). Sin embargo, la condición inicial está dada en el instante t,+l y la velocidad únicamente es conocida en instantes anteriores. Ello hace necesario construir un esquema numérico no estándar, de tipo mixto Runge-Kutta y multipaso: l LEMA 2. Sea 2n+l = xn+l + At u(xn+l , t,). (40) 1 Entonces, l es tal que lxn - ~(t,,t,+~,x~+')l 5 cat3 (42) O Por otra parte, el error local en la interpolación de la vorticidad debe ser también de tercer orden en h. Ello se consigue interpolando ésta localmente mediante polinomios de segundo grado. Concretamente, la vorticidad ha sido discretizada sobre el siguiente espacio de elementos finitos: l Wh = {wh E C'(B~)IW~,~ E P2,VT E E}. (43) Denotemos por {/.Lj}jEIN los puntos medios de los lados de los triángulos de Th. Entonces, una función wh E Wh queda unívocamente determinada por los valores 254 T. CHACON REBOLLO E 1. BLESS RANERO wh(~i)y Vi E IN y wh(pj), Vj E IN. (44) Denotaremos por ph : C0(Et2) + Wh al operador que interpola los valores dados en (44). Nuestro algoritmo puede ahora ser descrito como sigue: Algoritmo B Supongamos que wo E C2(Et2) y tiene soporte compacto. Aproximaremos w(.,t,) mediante una función U; de Wh. (a) Inicialización: i. Vorticidad: 6' = Ph 0 Wg7 ii. Velocidad: a0 = K * (ri O w:), (b) Iteración en tiempo: Supongamos dados tú:, tú;, . . ., W; E Wh. i. Velocidad discreta: ii. Característica discreta: pr: iii. Dado x E IEt2, donde 2 = x + Atú;(x). iv. Vorticidad: METODO DE ELEMENTOS FINITOS PARA ECUACIONES DE EULER 2D 255 Observación Es necesario ejecutar el primer paso de este algoritmo con un método de segundo orden, para conservar su orden de precisión. Observemos que en el Algoritmo B la velocidad discreta no está calculada mediante convolución del núcleo deBiot-Savart con la vorticidad discreta, sino con el interpolado afín a trozos de ésta. De esta forma se consigue una notable reducción en el número de operaciones necesarias para calcular la velocidad, sin que por ello el error de truncamiento local deje de ser de tercer orden. Esta afirmación se formaliza como sigue: LEMA 3. Supongamos que wo E C3(JR2) con soporte compacto. Definamos 6; = ~h(~(.,t?)), .u.; = K * (rhOW;). Sean q+' E Wh y %+' la velocidad y vorticidad calculadas mediante un paso de tiempo del Algoritmo B, a partir de los valores Wh, 6; y 6;-l. Entonces existe una constante C que depende únicamente de wo y T, tal que la+' Wh - w(-,tn+l)lL,,(R2) < C(h3 + At3); "h < c(h3 + At3); 14+' - ,+IL( - o De esta forma, el Algoritmo B tiene efectivamente un error de consistencia local de tercer orden. Sin embargo, este algoritmo no es ahora estable en norma uniforme. En efecto, el operador de interpolación ph, en general incrementa la norma uniforme de la función interpolada. Por otra parte, ya hemos comentado anteriormente, que en el algoritmo Al, el paso más costoso computacionalmente es el cálculo de la velocidad en los N vértices de la triangulación, que requiere en total O(N2) operaciones. Como se puede apreciar, en el algoritmo B este paso es aún más costoso, para calcular la velocidad discreta serán necesarias del orden de ocho veces este número de operaciones. Por esta razón hemos utilizado una técnica de cálculo rápido de la velocidad, que reduce este número de operaciones de O(N2) a O(N log2(¡V)) (Cf. [4] y [lo]). Un ejemplo numérico En esta sección, analizaremos los resultados de la ejecución práctica del algoritmo B en el caso test empleado en la Sección 6. En él, se confirma numéricamente la convergencia de segundo orden, previsible a partir de las estimaciones de consistencia. Para discretizar la vorticidad hemos utilizado triangulaciones uniformes, con tallas h=L h=' 8 7 12 Y h = h. Para medir los errores hemos utilizado la misma seminorma que en la sección referente a Ejemplos Numéricos. Las Figuras 5 y 6 muestran las curvas de error en velocidad y vorticidad respectivamente, con At = 9. Se observa una notable disminución de la difusión numérica, ya que el error permanece prácticamente constante.