Algunos elementos para la construcción de un Método de Multiescala Variacional “a posteriori”

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Algunos elementos para la construcci´on de un M´etodo de Multiescala Variacional “a posteriori”. Tom´ as Chac´ on Rebollo1, Antonio Dom´ ınguez Delgado 2 1Dpto. E.D.A.N., Universidad de Sevilla, Aptdo. 1160, E-41080 Sevilla. E-mail: [email protected]. 2Dpto. de Matem´atica Aplicada I, Universidad de Sevilla, Sevilla. E-mail: [email protected]. Palabras clave: M´etodos Multiescala, SUPG, convecci´on dominante, estabilizaci´on. Resumen El presente trabajo est´a dedicado a la resoluci´on de problemas de la Mec´anica de Fluidos incompresibles en r´egimen de convecci´on dominante. Nuestro prop´osito es presentar una t´ecnica de post-proceso de una soluci´on Galerkin con E.F. P1 inicial oscilante, que permite recuperar a partir de la misma, otra soluci´on no oscilatoria, y que adem´as es una aproximaci´on de segundo orden de la soluci´on continua. Presentamos tambi´en, una extensi´on de esta t´ecnica al caso no lineal, proporcionando un m´etodo eficaz para la resoluci´on de choques. Como extensi´on del M´etodo Multiescala Variacional debido a Hughes, estamos introduciendo un M´etodo Multiescala Variacional “a posteriori”. 1 Introducci´on Consideramos el problema de convecci´on-difusi´on en una dimensi´on: (a u0−ν u00 =fen [0,1] u(0) = 0, u(1) = 1 (1) donde la constante arepresenta la velocidad de convecci´on del fluido y νes el coeficiente de difusi´on. Supondremos a6= 0, ν > 0. 1 Tom´as Chac´on, Antonio Dom´ınguez Consideramos una discretizaci´on de (1) mediante una aproximaci´on Galerkin con elementos finitos, usando funciones afines a trozos sobre una partici´on 0 = x0< x1< . . . < xN= 1 del intervalo [0,1] de tama˜no uniforme h= 1/N. El espacio aproximante ser´a entonces: ˜ Xh=nxh∈C0[0,1] |xh|]xj−1,xj[∈P1,1≤j≤No. Denotemos por {ϕi}N i=0 la base can´onica del espacio considerado, es decir cada ϕies af´ın en cada subintervalo ]xj−1, xj[,1≤j≤Nyϕi(xj) = δij,0≤i, j ≤N. Aproximamos la soluci´on ude (1) mediante, uh=PN i=0 uiϕi∈˜ Xh,con u0= 0, uN= 1 y ui=uh(xi). Denotemos por Peh=|a|h 2νel llamado n´umero de P´eclet de malla asociado a la discretizaci´on. Es bien conocido que cuando Peh<1 la sucesi´on {ui}N i=0 satisface el principio del m´aximo discreto, justo lo contrario que sucede cuando Peh>1, y por ello en este caso, es decir, en situaciones de convecci´on dominante, aparecen oscilaciones esp´ureas en la soluci´on discreta de de (1). El m´etodo SUPG (Streamline Upwind Petrov-Galerkin) proporciona una t´ecnica para eliminar las oscilaciones esp´ureas antes descritas en la discretizaci´on Galerkin de las ecuaciones de convecci´on-difusi´on 1D. Tambi´en permite amortiguarlas en problemas de convecci´on multidimensional. Recordamos que la estabilizaci´on propuesta por el m´etodo SUPG, [4], equivale a discretizar por diferencias finitas centradas la ecuaci´on a u0−(ν+νnum)u00 =f , (2) con νnum =ν Pehξ(Peh), donde ξ(α) = coth(α)−1 α. La elecci´on indicada de νnum, se deduce imponiendo que sobre los nodos de la malla, la soluci´on apr´oximada uhque proporciona el m´etodo de Galerkin aplicado a (2) coincida con la soluci´on exacta de (1), es decir uh(xi) = u(xi),0≤i≤N. Por tanto, con la t´ecnica empleada en la construcci´on del m´etodo SUPG, se busca conseguir un m´etodo m´as estable y de precisi´on ´optima. En [8], se demuestra que el esquema num´erico (2) proporcionado por el m´etodo SUPG es equivalente a a˜nadir al espacio ˜ Xh, una burbuja “´optima”por elemento y ”condensarla”despu´es. Por “burbuja”sobre el elemento ]xi−1, xi[ por ejemplo, entendemos cualquier funci´on de H1 0( ]xi−1, xi[ ), extendida por cero a todo [0,1]. ”Condensar”las burbujas significa eliminar expl´ıcitamente los grados de libertad ”burbuja”y obtener la ecuaci´on equivalente satisfecha por la componente P1 de la soluci´on. 2 Algunos elementos para un M´etodo de Multiescala Variacional “a posteriori” Este m´etodo es un caso particular del M´etodo Multiescala Variacional introducido por Hughes en 1995 , [7]: la componente P1 de la velocidad constituye las ”escalas resueltas”, mientras que la componente ”burbuja”forma las escalas ”no resueltas”, cuyo efecto es incorporado a la ecuaci´on satisfecha por las escalas resueltas. Este efecto es una viscosidad num´erica que en el caso ´optimo es la precisa para que la soluci´on proporcionada por el m´etodo SUPG coincida con la soluci´on exacta en los nodos de la malla. 2 Autoestabilizaci´on ´optima: caso de coeficientes constantes. Consideremos una descomposici´on de Xhen la forma Xh=YhLZh, donde YhyZhson subespacion de Xh, siendo Yhde elementos finitos, con YhTZh={0}. Si fes constante, es posible determinar Zhde forma que la aproximaci´on filtrada de la soluci´on oscilante uhcoincida en los nodos de la partici´on de paso 2hcon la soluci´on exacta udel problema de convecci´on difusi´on con coeficientes constantes (1), de modo an´alogo al m´etodo SUPG. Tomamos como subespacio Yh, el espacio de las funciones afines a trozos sobre la partici´on 0 = x0< x2< . . . < xN−2< xN= 1 del intervalo [0,1], y de tama˜no uniforme 2h. Construimos el espacio Zhde base {ϕ1, ϕ3, . . . , ϕ2M−1}de modo que se tenga la descomposici´on uh= M X i=1 yiφi+ M X i=1 z2i−1ϕ2i−1,con yi=u(x2i) para 1 ≤i≤M. (3) Para ello, el soporte m´ınimo que debe tener ϕ2i−1,2≤i≤M, es el intervalo [x2i−3, x2i], y el soporte m´ınimo de ϕ1debe ser [x0, x2]. En este caso, (3) puede ser resuelto anal´ıticamente de forma ´unica. Definimos las siguientes funciones de base, ϕ2i−1(x) =            −1 h(x−x2i−1) + 1,si x∈[x2i−1, x2i], a2i−1(x−x2i−1)+1,si x∈[x2i−2, x2i−1[, (1 h−a2i−1)(x−x2i−3),si x∈[x2i−3, x2i−2[ y i≥2, 0 en otro caso. (4) Imponiendo la condici´on de que la soluci´on discreta coincida en los nodos pares con la soluci´on exacta, obtenemos z1=uh(x1)−1 2u(x2), a1=1 h, y para i= 2,3, . . . , M, z2i−1=uh(x2i−1)−1 2(u(x2i−2) + u(x2i)), a2i−1=−1 h·uh(x2i−2)−u(x2i−2) z2i−1 −1¸. 3 Tom´as Chac´on, Antonio Dom´ınguez ||||| |||| - - ±1 1 2n2n–12n–22n–32i2i–12i–22i–3210 xxxxxxxxxxx 2n–1 ϕ 2i–1 ϕ 1 ϕ Figura 1: Base ´optima de Zhpara la ecuaci´on de Convecci´on-Difusi´on (1), para Peh= 10. En la Figura 1 aparecen representadas las funciones de base de Zh. Podemos observar c´omo las pendientes a2i−1se incrementan al acercarse el punto x2i−1a la capa l´ımite. Es claro que Zhdepende de la soluci´on exacta u. Ahora bien, dado que u(x2i) = 1−e2iP eh 1−eP eh, en realidad Zhdepende ´unicamente del n´umero de P´eclet de malla. Esto permite generalizar la definici´on de Zha problemas con velocidad de convecci´on variable, cambiando el n´umero de P´eclet de malla que define cada pendiente a2i−1, por un n´umero de P´eclet de malla promedio en un entorno de x2i−1. Por otro lado observemos que con esta construcci´on Zhes subespacio de H1 0(]0,1[). Adem´as, es inmediato que Xh=Yh+Zh. En efecto, dada una funci´on xh∈Xh, se tiene xh= M X i=1 ˆyiφi+ M X i=1 ˆz2i−1ϕ2i−1,(5) donde los coeficientes ˆyi,ˆz2i−1, est´an dados de forma recursiva, para i=M, M −1, . . . , 1 por, ˆz2i−1=xh(x2i−2)−2xh(x2i−1) + ˆyi −1−ha2i−1 ,ˆyi−1= 2(xh(x2i−1)−ˆz2i−1)−ˆyi. Observemos que las pendientes a2i−1son positivas, como se comprueba sin dificultad, por lo que los denominadores −1−ha2i−1son estrictamente negativos, y por tanto los valores de ˆz2i−1, est´an bien definidos para 1 ≤i≤M. Concluyendo, podemos considerar que el m´etodo presentado con la elecci´on realizada del espacio Zh, puede ser visto como un m´etodo SUPG “a posteriori”, en el que para el caso de coeficientes constantes se recupera la soluci´on exacta sobre la malla de paso 2h. 4 Algunos elementos para un M´etodo de Multiescala Variacional “a posteriori” En este sentido estamos introduciendo un m´etodo Multiescala Variacional “a posteriori”. Nuestras escalas “resueltas”constituyen el espacio Yh, mientras que las “no resueltas”forman el espacio Zh. 3 Ecuaciones de convecci´on-difusi´on no lineales. En esta secci´on aplicamos el procedimiento de autoestabilizaci´on ´optima a un problema de convecci´on difusi´on evolutivo no lineal, concretamente la ecuaci´on del tr´afico. Consideramos el siguiente problema:      u,t(x, t) + F,x(u(x, t)) = ν u,xx(x, t), x ∈[0,1], t ∈[0, T], u(x, 0) = u0(x), x ∈[0,1], (6) donde urepresenta la densidad del tr´afico y la funci´on de flujo Fes dos veces diferenciable, y verifica, F(0) = 0, F(umax) = 0 F(u)>0 y F00(u)<0 en 0 < u < umax .(7) El flujo F(u) depende de la v´ıa en que circula el tr´afico. Aqu´ı hemos considerado el caso modelo F(u) = u(1 −u). Si consideramos un problema de Riemann, u0(x) =      uLif x∈[0,1/2], uRif x∈[1/2,1], (8) con uL< uR,la soluci´on est´a formada por dos estados constantes unidos mediante un choque con difusi´on, que se desplaza hacia la izquierda . Para resolver (6), discretizamos en tiempo mediante el esquema semi-impl´ıcito de Euler, obteniendo una ecuaci´on con la estructura 1 ∆tun+1 +a(x)¡un+1¢0−ν(un+1) 00 =fn+1 +1 ∆tun un+1(0) = uL, un+1(1) = uR,   (9) donde a(x) = F0(un). Discretizamos ahora (9) por el m´etodo Galerkin est´andar, usando para ello el espacio 5 Tom´as Chac´on, Antonio Dom´ınguez de Elementos Finitos Xh. De esta forma, en cada etapa de tiempo resolvemos el problema: Obtener un+1 h∈Xhtal que ∀vh∈Xh, 1 ∆t(un+1 h, vh)+(a¡un+1¢0, vh)−ν((un+1) 0 , v0 h) = =hfn+1, vhi+1 ∆t(un h, vh) un+1 h(0) = uL, un+1 h(1) = uR.                  (10) Al aproximar la soluci´on mediante el m´etodo de Galerkin usando Elementos Finitos P1 Lagrange, la aproximaci´on que se obtiene para altos valores del P´eclet de malla, viene afectada de fuertes oscilaciones en torno al choque, (ver Fig. 3). Es posible aplicar a la soluci´on de (10) el filtrado ´optimo descrito en la secci´on anterior. Con ello se eliminan las oscilaciones esp´ureas. Sin embargo, es mucho mejor estrategia aplicar este filtrado s´olo en el instante de tiempo en que interese calcular la soluci´on con precisi´on: ALGORITMO FILTRADO ´ OPTIMO CASO NO LINEAL: Inicializaci´on. Se determinan N, 4t=T/N, u0 hy4t∗: n´umero cercano a la precisi´on del ordenador (pero mayor). I: Iteraci´on en tiempo sin filtrar. Se resuelve el sistema (10), obteniendo un+1 ha partir de un h, para n= 0,1, . . . , N −1. II: Filtrado. II.1: Se resuelve el sistema (10), obteniendo uN∗ ha partir de uN husando el paso de tiempo 4t∗. II.2: Se efect´ua el filtrado “´optimo”de uN∗ h, usando la base “´optima”de Zhdeterminada mediante los promedios de aN∗ hsobre cada intervalo [x2i−2, x2i]. II.3: Si la soluci´on presenta oscilaciones se hace uN h=uN∗ hy se vuelve a II.1. Observemos de nuevo que el filtrado no se realiza sino hasta el mismo instante de tiempo Ten el que queremos obtener la soluci´on. 4 Tests Num´ericos. Presentamos en primer lugar un test para la ecuaci´on de convecci´on difusi´on lineal (1), con t´ermino fuente lineal. Como se aprecia en la figura, nuestra t´ecnica de postproceso ´optimo permite recuperar la soluci´on exacta sobre la malla de paso 2h. El resultado obtenido se muestra en la Fig.2 6 Algunos elementos para un M´etodo de Multiescala Variacional “a posteriori” Figura 2: Soluci´on exacta (linea marcada por +), de la ecuaci´on de Convecci´on-Difusi´on con t´ermino fuente lineal y condiciones de contorno homog´eneas, aproximaci´on por MEFGalerkin P1(linea marcada por * ) y aproximaci´on filtrada con Zh´optimo (linea marcada por ◦), para Peh= 10. 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Figura 3: Soluci´on Galerkin (l´ınea de puntos), soluci´on filtrada (◦), y soluci´on exacta (solid line), para N= 60, ν = 0.001, t = 0.1. 7 Tom´as Chac´on, Antonio Dom´ınguez En el segundo test se ha aplicado la t´ecnica de postproceso al problema no lineal de tr´afico dado por (6) con las condiciones iniciales (8). Los datos del test son: N= 60, ν = 0.001, T = 0.1,∆t= 1/200, uL= 0.2, uR= 0.9 De aqu´ı Peh= 16.66.El n´umero de filtrados realizados (paso II.3), una vez alcanzado el tiempo final es de 6. Los resultados obtenidos se presentan en la Figura 3. Agradecimientos Este trabajo ha sido parcialmente financiado por el Proyecto del Plan Nacional de I+D MTM2006-1750. Referencias [1] C. Baiocchi, F. Brezzi and L. P. Franca (1993), Virtual Bubbles Galerkin-least-squares type methods (Ga.L.S.) Comput. Methods Appl. Mech. Engrg., 105, pp. 125-141. [2] F. Brezzi, M. O. Bristeau, L. P. Franca, M. Mallet and G. Rog´e (1992), A relationship between stabilized finite element methods the Galerkin method with bubble functions, Comput. Methods Appl. Mech. Engrg., 96, pp. 117-129. [3] F. Brezzi and J. Douglas (1988), Stabilized mixed methods for the Stokes problem, Numer. Math. 53, pp 225-236. [4] A. N. Brooks and T. J. R. Hughes (1982), Streamline/Upwind Petrov-Galerkin formulations for convection dominated flows, with particular emphasis on the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg., 32, pp. 199-259. [5] Tom´as Chac´on Rebollo and Antonio Dom´ınguez Delgado (2000), A Unified Analysis of Mixed and Stabilized Finite Element Solutions of Navier-Stokes Equations, Computer Methods in Applied Mechanics and Engeneering, Vol. 182, pp: 301-331. [6] T. J. R. Hughes, L.P. Franca and M. Balestra (1986), A new Finite Element formulation for CFD : V. Circumventing the Brezzi-Babuˇska condition: A stable Petrov-Galerkin formulation of the Stokes problem accommodating equal-order interpolations. Computer Methods in Applied Mechanics and Engeneering. Vol. 59, pp. 85-99. [7] T. J. R. Hughes (1995), Multiescale phenomena: Grenn’s functions, the Dirichlet-to-Neumann formulation, subgrid-scale models bubbles and the origin of stabilized methods. Computer Methods in Applied Mechanics and Engeneering. Vol. 127, pp. 387-401. [8] Alessandro Russo (1996), Bubble stabilization of finite element methods for the linearized incompressible Navier-Stokes equations., Computer Methods in Applied Mechanics and Engeneering, Vol. 132, pp.335-354. 8