Formulación de tipo Petrov-Galerkin de algunos métodos distributivos: Aplicación a las ecuaciones de Navier-Stokes

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Formulaci´on de tipo Petrov-Galerkin de algunos m´etodos distributivos: Aplicaci´on a las ecuaciones de Navier-Stokes. T. Chac´ on Rebollo1, M. G´ omez M´ armol1, G. Narbona Reina2 1Dpto. E.D.A.N., Universidad de Sevilla, Aptdo. 1160, 41012 Sevilla. E-mails: [email protected], [email protected]. 2Dpto. Matem´atica Aplicada I, Universidad Sevilla, 41012 Sevilla. E-mail: [email protected]. Palabras clave: M´etodos distributivos, Petrov-Galerkin, Navier-Stokes Resumen En este trabajo estudiamos la resoluci´on de las Ecuaciones de Navier-Stokes estacionarias mediante m´etodos distributivos no lineales. Formulamos estos m´etodos como m´etodos de tipo Petrov-Galerkin, en un contexto de discretizaci´on por el m´etodo de los elementos finitos. Utilizamos funciones tests descentradas “corriente arriba”para el tratamiento del t´ermino de convecci´on. Esta formulaci´on nos permite realizar el an´alisis de los m´etodos distributivos que consideramos como una extensi´on del an´alisis est´andar. Presentamos resultados de existencia de soluci´on del problema discreto, convergencia y estimaciones de error. Por ´ultimo, presentamos algunos test num´ericos resueltos mediante un esquema de tipo distributivo no lineal, el PSI. Estos tests muestran un comportamiento resistente a la generaci´on de oscilaciones par´asitas, y una mayor exactitud que un m´etodo de las caracter´ısticas de primer orden. 1. Ecuaciones de Navier-Stokes En esta secci´on desarrollamos la aproximaci´on del problema estacionario de NavierStokes. Para ello consideramos estas ecuaciones de tipo convecci´on-difusi´on, con el objetivo de realizar la discretizaci´on del t´ermino de transporte utilizando el m´etodo PSI. Sea Ω un dominio acotado en Rd(d= 2 o 3) con frontera Lipschitz Γ = ∂Ω. Consideramos el siguiente problema:    Hallar u: Ω −→ Rd, p : Ω −→ Rtal que u· ∇u−ν∆u+∇p=f, ∇ · u= 0 en Ω; u= 0 sobre Γ, (1) 1 T. Chac´on, M. G´omez, G. Narbona siendo uel campo de velocidad, pla presi´on y νel coeficiente de difusi´on. Y denotamos por f: Ω →Rd, el t´ermino fuente. Consideramos los espacios V= (H1 0(Ω))dyL2 0(Ω) donde vamos a buscar la soluci´on u yp, respectivamente. Suponemos fen (H−1(Ω))d. Por ´ultimo, definimos la forma trilineal basociada al primer t´ermino de la ecuaci´on: b(u, v, w) = RΩ(u· ∇v)w dx, con u, v, w ∈V. Y denotamos su norma por M. Consideramos la formulaci´on variacional est´andar del problema (1):    Hallar u∈V, p ∈L2 0(Ω) tal que b(u, u, v) + ν(∇u, ∇v)−(p, ∇ · v) =< f, v > ∀v∈V; (∇ · u, q) = 0 ∀q∈L2 0(Ω). (2) Se puede probar que el problema (2) admite soluci´on (Cf. [6], [8]) y que adem´as es ´unica bajo la condici´on Mkfk−1< ν2. A partir de ahora supondremos que Ω es un dominio poligonal. Sea Thuna triangulaci´on de Ω. Definimos los siguientes espacios discretos: V∗ h={vh∈ C0(¯ Ω)/ vh|T∈P1,∀T∈ Th}, Vh={vh∈V∗ h/ vh= 0 sobre Γ}.(3) Nuestra aproximaci´on interna {Vh}h>0para el espacio de velocidades Ves Vh= (Vh)d. Consideramos tambi´en una aproximaci´on interna por tipo elementos finitos {Mh}h>0del espacio de presiones L2 0(Ω) tal que se satisface una condici´on inf-sup: ∃ˆ β > 0 tal que ˆ βkrhk ≤ sup vh∈Vh (rh,∇ · vh) k∇vhk0 ,∀rh∈Mh.(4) 1.1. Esquema PSI En esta secci´on describimos la idea general para la construcci´on de los m´etodos distributivos, y en particular para el esquema PSI. Por simplicidad en la notaci´on lo haremos en dimensi´on 2. Los m´etodos distributivos se centran en el tratamiento del t´ermino de convecci´on. Vamos a considerar el flujo convectivo de la magnitud ρtransportada por la velocidad u, definido por qconv =uρ. El balance total del flujo convectivo que atraviesa la frontera de un elemento T∈ Thviene dado por ΦT=Z∂T qconv ·n=ZT ∇ · (uρ) = ZT u· ∇ρ, donde hemos utilizado que ∇ · u= 0. La idea b´asica de los m´etodos distributivos es repartir el flujo ΦTentre los v´ertices de los elementos vecinos a T. Esta distribuci´on se hace mediante unos coeficientes que denotamos por {βT j}N j=1, de forma que el flujo enviado al nodo bjes Φj=X T∈Th βT jΦT.(5) 2 Formulaci´on Petrov-Galerkin de algunos m´etodos distributivos Normalmente, el flujo generado en el elemento Ts´olo se distribuye entre los v´ertices de este elemento, as´ı que podemos escribir: βT j= 0 si bjno es un v´ertice de T. As´ı, si denotamos por Ejel conjunto de todos los elementos de los cuales bjes v´ertice, podemos escribir: Φj=X T∈Ej βT jZT u· ∇ρ. (6) Para obtener un flujo conservativo y un m´etodo estable L∞, se consideran las siguientes propiedades sobre los coeficientes de distribuci´on (Cf. [3]): d+1 X j=1 βjT T= 1, βT jT≥0,para todo nodo bj,para todo tri´angulo T∈ Th.(7) El esquema PSI se construye a partir de un m´etodo distributivo lineal llamado N-esquema, el cual viene determinado por la definici´on del flujo ΦT i(sh) = βT iΦT(sh),para sh∈Vh. Para ello vamos a introducir el vector normal interior a la frontera de T, opuesto al nodo bi, nT i=d|T| ∇ϕiT,(8) y los valores KT i=1 d¯uT·nT i,con ¯uT=1 |T|ZT u. De manera que para el N-esquema tenemos, ΦT i(sh) = 3 X i=1 cT ij (sT i−sT j), cij = (KT i)+MT(KT j)−,(9) donde (KT i)+= m´ax{KT i,0},(KT j)−= m´ın{KT j,0}, MT= 3 X j=1 (KT j)−, ysT 1,sT 2,sT 3son los valores de shen los v´ertices de T. Observar que el signo de KT iindica si el v´ertice biest´a “corriente abajo” con respecto au(KT i≤0). Para que el m´etodo sea estable, los nodos que est´an “corriente arriba” en el tri´angulo Tno reciben ninguna aportaci´on del flujo de este elemento. En algunos casos, el flujo ΦT(sh) puede anularse, mientras que el flujo ΦT i(sh) toma un valor finito, y por tanto, el coeficiente βT ino estar´ıa definido. Es por esta raz´on por la que el N-esquema no puede escribirse bajo una formulaci´on Petrov-Galerkin (ver [5]). El esquema PSI se construye, como hemos dicho antes, a partir del N-esquema, buscando unas nuevas funciones de flujo Φ∗ itales que Φ∗ i= (1 −µi) Φi,para unos coeficientes 0 ≤µi≤1, i = 1,· · · , d + 1 (10) 3 T. Chac´on, M. G´omez, G. Narbona d+1 X i=1 Φ∗ i= Φ.(11) Los coeficientes β∗ i=Φ∗ i Φest´an acotados independientemente de la malla.(12) (Hemos suprimido el super´ındice Ty la dependencia de shpor claridad en la notaci´on). Se prueba que µison funciones continuas de Φi, y que β∗ icumplen la propiedad (7). Observar adem´as que tanto los coeficientes de distribuci´on como las funciones de flujo Φ∗ y Φ dependen de sh. La resoluci´on del problema (10)-(12) y por tanto la definici´on expl´ıcita del esquema PSI, puede encontrarse en [3], no la detallamos aqu´ı, ya que no es relevante en nuestro an´alisis. 1.2. Formulaci´on Petrov-Galerkin Podemos dar una formulaci´on abstracta de tipo Petrov-Galerkin del m´etodo PSI y de otros m´etodos distributivos no lineales. Al igual que en [5], definimos un nuevo espacio discreto cuyas funciones de base nodales denotamos por λ1, λ2,· · · , λM(dependiendo de un elemento σhde Vh): Wh(σh) = span{λ1(σh), λ2(σh), . . . , λM(σh)}.(13) En particular para el m´etodo PSI, dado un nodo bide la malla, definiremos la funci´on de base asociada como sigue: λi(σh) = ½β∗ i T(σh) si bi∈T; 0 en otro caso. Observemos que la dependencia respecto de la funci´on σhse debe al car´acter no lineal del m´etodo PSI. Consideramos tambi´en el operador de interpolaci´on asociado, que toma valores en Wh, y que denotamos por Πσh: Πσh:C0(Ω) −→ Wh z−→ Πσhz= M X i=1 z(bi)λi,(14) donde {bi}M i=1 denota los nodos de la malla situados en el interior del dominio. Llamaremos a Πσhel operaci´on de interpolaci´on distribuida generado por la funci´on σh. La extensi´on de esta definici´on al espacio vectorial Vhse realiza por componentes, de la siguiente forma. Definimos Wh(wh) = Wh(w1h)×Wh(w2h) y denotamos por Πwhel operador vectorial de interpolaci´on distribuida generado por wh∈Vh. Utilizaremos las funciones de Wh(wh) como funciones test para el t´ermino de convecci´on. Definimos ahora la forma discreta asociada al t´ermino de convecci´on bh(rh) : Vh×Vh×Vh7→ Rcomo bh(rh;uh, vh, wh) = ZΩ (uh·∇vh)Πrhwh, para una funci´on 4 Formulaci´on Petrov-Galerkin de algunos m´etodos distributivos dada rh∈Vh. Escribimos entonces la aproximaci´on variacional del problema (2):    Hallar uh∈Vh, ph∈Mhtal que bh(uh;uh, uh, vh) + ν(∇uh,∇vh)−(ph,∇ · vh) = hf, Πuhvhi ∀ vh∈Vh; (∇ · uh, qh) = 0 ∀qh∈Mh. (15) Hemos discretizado tambi´en de forma descentrada el t´ermino fuente con el fin de obtener un esquema bien equilibrado al segundo orden en r´egimen de convecci´on dominante, ver [5]. Para que este t´ermino hf, Πuhvhitenga sentido basta considerar fen un espacio m´as regular que V0, por ejemplo f∈L2(Ω)d. 2. Existencia, convergencia y estimaciones de error Nuestro an´alisis de existencia de soluciones del problema discreto (2), convergencia y estimaciones de error est´a basado en ciertas propiedades satisfechas por nuestra formulaci´on : Propiedad 1 Para todo elemento T∈ Thy para todo σh∈Vh,λT i(σh)≥0, i= 1,· · · , d+1 y d+1 X i=1 λT iT(σh) = 1,donde por iTdenotamos el ´ındice global correspondiente al ´ındice local i, del elemento T, y siendo λT iT(σh)la restricci´on de λiT(σh)al elemento T. ¤ Definamos la matriz de convecci´on asociada a una velocidad vh∈Vhy a un elemento σh∈Vhdenotada por C(vh;σh)∈MM×M(R) (el espacio de las matrices reales de dimensi´on M×M) como sigue: Cij(vh;σh) = ZΩ (vh· ∇ϕj) Πσhϕi=ZΩ (vh· ∇ϕj)λi(σh), Entonces Propiedad 2 La matriz C(vh;σh)es semi-definida positiva para cualesquiera vh∈Vh yσh∈Vh.¤ Por ´ultimo se tiene una hip´otesis t´ecnica sobre el comportamiento de la matriz de convecci´on respecto a sus argumentos vhyσh: Propiedad 3 La matriz de convecci´on C(vh;σh)es una matriz continua de Vh×Vh al espacio MM×M(R).¤ Bajo las Propiedades 1 y 2, la forma bes continua sobre H1 0(Ω)3ysemidefinida positiva. Observaci´on 1 Una t´ecnica habitual para la aproximaci´on de las ecuaciones de NavierStokes por m´etodos mixtos en elementos finitos, consiste en reemplazar la forma bpor una forma antisim´etrica e bque satisface e b(uh;vh, vh) = 0. En nuestro caso, la forma bhno satisface esta condici´on, pero es semi-positiva. Esta propiedad nos permite realizar el an´alisis del problema discreto usando el an´alisis est´andar con algunas modificaciones adecuadas. En particular, el an´alisis de existencia de soluciones del problema discreto se sigue del Teorema del Punto Fijo de Brouwer: 5 T. Chac´on, M. G´omez, G. Narbona Teorema 2.1 Supongamos f∈Lr(Ω)dpara r > rmin, con rmin = 1 para d= 2 y rmin = 6/5para d= 3. Entonces el problema (15) admite al menos una soluci´on que satisface la estimaci´on: k∇uhk0≤ν−1Krkfkr kphk0≤ˆ β−1(Nν−2K2 rkfkr+ 2Kr)kfkr (16) siendo ˆ βla constante de la condici´on inf-sup (4), Kruna constante que depende de r, y N= sup 0<h≤h0 sup u,v,w∈Vh bh(r;u, v, w) k∇uk0k∇vk0k∇wk0 . ¤ El an´alisis de convergencia y estimaciones de error se siguen de una modificaci´on bastante elaborada del an´alisis est´andar de aproximaci´on de las Ecuaciones de NavierStokes mediante m´etodos mixtos. Tenemos los siguientes resultados. Teorema 2.2 Bajo las hip´otesis del Teorema 2.1, existe una subsucesi´on de la sucesi´on {(uh, ph)}h≥0, soluci´on de (15) que converge fuerte en V×L2 0(Ω) a una soluci´on del problema (2). Si el problema continuo tiene una ´unica soluci´on, entonces toda la sucesi´on converge a ella. ¤ Teorema 2.3 Supongamos ciertas las hip´otesis del Teorema 2.2 y adem´as suponemos que Nkfk−1< ν2yMkfk−1< ν2. Entonces existe una constante positiva Ctal que k∇(u−uh)k0+kp−phk0≤Chd1(u, Vh) + d0(p, Mh) + h1−d ˆsi,∀ˆs∈(d, rmax] con rmax =½∞d= 2; 6d= 3. donde d1(u, Vh) = ´ınf vh∈Vh k∇(u−vh)k0yd0(p, Mh) = ´ınf qh∈Mh kp−qhk0.¤ 3. Resultados num´ericos En esta secci´on presentamos algunos test num´ericos cl´asicos para el problema de Navier-Stokes. Se han obtenido utilizando el m´etodo distributivo PSI, que comparamos con la soluci´on obtenida por el m´etodo de las caracter´ısticas. Observamos una mayor precisi´on y una escasa formaci´on de oscilaciones par´asitas del m´etodo PSI, especialmente en zonas de flujo de fuerte gradiente. Para la resoluci´on num´erica hemos utilizado del software Freefem++ (www.freefem.org). 3.1. Test 1: Problema de la cavidad Consideramos como dominio el cuadrado unidad, Ω = [0,1]2. Y resolvemos el problema (15) para ν= 0,001 y f= 0. Como condici´on inicial tomamos u= 0 en todo el dominio y como condici´on de contorno: u=½(1,0) sobre y= 1; (0,0) en otro caso. 6 Formulaci´on Petrov-Galerkin de algunos m´etodos distributivos Figura 1: Test 1: M´etodo de las Caracter´ısticas (izquierda), Esquema PSI (derecha). En la Figura 1, representamos la soluci´on obtenida para tiempo final T= 2,4. Observamos el buen comportamiento del m´etodo PSI y una mejor aproximaci´on en la zona superior izquierda, donde los gradientes son fuertes. Esta diferencia es a´un m´as clara en el Test 2, como veremos a continuaci´on. 3.2. Test 2: Problema del escal´on Se trata de resolver mediante las ecuaciones de Navier-Stokes incompresibles la evoluci´on de un fluido que se mueve en un canal estrecho con un escal´on en la parte inferior. El dominio est´a representado en la Figura 2, donde hemos dividido la frontera en cuatro partes. Consideramos en este caso el coeficiente de difusi´on ν= 0,0025 y t´ermino fuente f= 0. Las condiciones de contorno vienen dadas por: u=   (4y(1 −y),0) sobre Γ1; (0,0) sobre Γ2∪Γ4; ν∂u ∂n +p n = 0 sobre Γ3. Como condici´on inicial tomamos: u0=½(4y(1 −y),0) sobre Γ1; (0,0) en otro caso. En este caso podemos apreciar una mayor diferencia entre las soluciones, notando que la soluci´on obtenida por el m´etodo PSI es m´as regular y captura mejor el torbellino formado tras el escal´on. En la Figura 3 mostramos la soluci´on para ambos m´etodos, para un tiempo final de T= 20 segundos. Agradecimientos Investigaci´on parcialmente financiada por Proyecto de Investigaci´on MEC MTM200601275. 7 T. Chac´on, M. G´omez, G. Narbona 1 3 2 4 Γ Γ Γ Γ Figura 2: Test 2: Dominio. Figura 3: Test 2: M´etodo de las Caracter´ısticas (arriba), Esquema PSI (abajo). Referencias [1] F. Clarke, Optimization and Nonsmooth Analysis, John Wiley & Sons, Inc., New York, 1983. [2] J. Simon. Compact sets in Lp(0, T ;B). Ann. Mat. Pura Appl., s´er. IV, CXLVI (1987), 65–96. [3] R. Abgrall, M. Mezine, Construction of second-order accurate monotone and stable residual distribution schemes for steady problems, Journal of Computational Physics 195, pp. 474-507 (2004). [4] S. Brenner, L. Ridgway Scott, The Mathematical Theory of Finite Element Methods, Springer (2002). [5] T. Chac´on Rebollo, M. G´omez M´armol, G. Narbona Reina, Numerical analysis of the PSI solution of the convection-diffusion problem through a Petrov-Galerkin formulation, Mathematical models and methods in applied sciences. En prensa. [6] V. Girault, P.A. Raviart, Finite Element Methods for Navier-Stokes Equations, Springer-Verlag, Berlin (1986). [7] O. Pironneau, On the transport-diffusion algorithm and its applications to the Navier-Stokes equations. Numer. Math. 38, no. 3, 309–332 (1981/82). [8] R. Temam, Navier-Stokes equations, North-Holland Publishing Company, (1977). 8