Métodos de elementos finitos estabilizados para flujos de fluidos incompresibles

Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 TEMAS MÉTODOS DE ELEMENTOS FINITOS ESTABILIZADOS PARA FLUJOS DE FLUIDOS INCOMPRESIBLES) TOMÁS CHACÓN REBOLLO DEPARTAMENTO DE ECUACIONES DIFERENCIALES yANÁLISIS NUMÉRICO UNIVERSIDAD DE SEVILLA el TARF1A, S/N. 41080 SEVll.LA e·mail: chaconQnumer.us.es 1 Introducción En este t.rabajo abordaremosciertosaspectosdel desarrollo de técnicas eficientes de simulación numérica de flujos de fluidos incompresibles mediante el Método de los Elementos Finitos. También nos ocuparemos del análisis numérico de las mismas. El estudio teórico de las ecuaciones relevantes que gobiernan los flujos de fluidos incompresibles, esencialmente las deStokes yNavier-Stokes, fue realizado y completado en los años 60 yprimeros 70. Como es bien conocido, se han demostrado resultados de existencia para las Ecuaciones de Navier-Stokes en dos y tres dimensiones de espacio. Sin embargo, únicamente se conocen resultados de unicidad en dos dimensiones, debido a la falta de regularidad de la solución en el caso tridimensional. Sirvan como referencia los trabajos, hoy clásicos, de Lions (1969) yTemam (1977). La resolución numérica de flujos de fluidos incompresibles presenta dos res-- tricciones de estabilidad de relevancia, bien conocidas por la comunidad de investigadores yusuarios en Mecánica de Fluidos computacional. Por una parte, la condición de incompresibilidad impone condiciones de compatibilidad bastante estrictas entre los espacios de discretización de velocidad ypresión empleados. Ilnvestigación fina.ncia.da parcialmente por los proyectos HCM ERBCHB lCI' 94.18'23 de la Unión Europea, yPB93-1196 de la DCICYT 5 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 Por otra, el carácter fuertemente convectivo de la. maYOl"la de los flujos de interés práctico produce la aparición de fenómenos de capa límite, zonas cerca de las paredes sólidas en que la velocidad presenta muy fuertes gradientes. La simulación de tales flujos mediante el uso directo del método de discretización Galerkin -Elementos Finitos produce en general soluciones numéricas presentando fuertes oscilaciones en las zonas de capa límite. Tales oscilaciones de~ sapare<::en únicamente si la talla de la malla es suficientemente pequeii.a. Ello origina discretizaciones con gran número de grados de libertad, resultando inabordables con los medios actuales en muchos casos de interés. Yello, incluso si los espacios de discretización de velocidad ypresión son compatibles. La primera técnica general para el tratamiento numérico de la condición de incompresibilidad fue introducida en Babuska (1973) yBrezzi (1974), para resolver la ecuación de Stokes. Se trata de introducir espacios de discretización diferentes para velocidad ypresión, de forma que satisfagan una condición de compatibilidad específica, conocida como de Bre-t.zj-BabuSka. Posteriormente, hubo ciertos esfuerzos tendentes a la utilización de bases del espacio de velocidades con divergencia nula, con objeto de eliminar la presión (ef. Griffits (1979), Hecht (1981)). Sin embargo, la construcción de tales bases resulta. una tarea ardua, ylos sistemas lineales resultantes están relativamente mal condicionados. Por ello, esta técnica es poco utilizada. Algo más recientemente, se introdujo en Brooks yHughes (1982) el método SUPG (Stl'eamline Upwind /Petrov-Galerkin). Se trata de una modificación del método estándar de Galerkin, en que las funciones test se añaden términos constantes por elemento, consiguiendo estabilizar no sólo la discretización de la presión, sino también el operador de convección "corriente arriba". Este método puede también ser visto como una técnica de Lagrangiano aumentado, consistente en "aumentar" la formulación Galerkin estándar con un término que trata de minimizar una norma del residuo. Adaptaciones de tal formulación fueron utilizada en Hughes et al. (1986), Douglas yWang (1989), Franca yFrey (1992) yotros para construir Uknicas estables, utilizando virtualmente cualquier palo de espacios de Elementos Finitos para velocidad ypresión. Al mismo tiempo, cuando se aplican ala resolución de las ecuaciones de Navier-Stokes, estas técnicas incorporan la difusión corriente arriba del método SUPG, con lo que consiguen estabilizar también el cálculo de flujos con convección dominante. Curiosamente, algunos estudios recientes muestran que existe una reladón muy estrecha entre los métodos de estabilización yalgunos de los métodos clásicos arriba mencionados -hasta el punto de ser equivalentes en ciertos casos. 6 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 En este trabajo nos ocuparemos de la descripción y el análisis numérico de los métodos de estabilización más utilizados. También de la relación entre éstos yla utilización de pares de espacios que satisfagan 1& condición de Brezú- Babu!ka. En la Sección 2repasaremos la teoría clásica sobre aproximación de las ecuaciones de 5tokes mediante el M.E.F.. La Sección 3 estará dedicada. al análisis de algunos métodos de estabilización, aplicados alas ecuaciones de 5t.okes. En la Sección 4estudiaremos la equivalencia entre el método SUPC yun método clásico usando funciones "burbuja". Por último, en la. Sección 5veremos cómo esta. Cormulación puede ser utilizada para realizar el análisis numérico del método SUPC en el caso de las ecuaciones de Navier-Stokes. Como referencias generales sobre métodos clásicos, citemos los libros de Temam (1977), Girault -Raviart (1986) yBrezzi yFortin (1991). Existen pocas referencias generales sobre métodos de estabilizo.ciÓn. Citemos, por ejem· plo, el artículo de Franca, Aughes yStenberg (1993), en que puede encontrarse un enfoque bastante general yabundantes referencias. 2 Teoría clásica de aproximación de las ecuaciones de Stokes por el M. E. F. En esta Sección harem06 un breve repaso sobre la teoría referente aexistencia yunicidad de soluciones de las ecuaciones de Stokes, asr como sobre su apro.'dmación mediante el M.E.f. En particula.r, recordaremos el papel de la condición inC· sup discreta en la convergencia de estas aproximacic;lOes. Consideremos un dominio acotado n e Rd (d =2 o 3), con frontera r. Nos damos un campo de velocidades "vector" u : n _ Rd. que supondremos con divergencia. nula.. Nuestro propósito es resolver numéricamente el siguiente problema de contorno para las ecuaciones de Stokes estacionarias, incluyendo un término convectivo : Obtener y: n _ Rd, u . 'Vy - vó.y +'Vp = f, p,O-R "·y~O y=O tales que } en n, sobre r (1 ) Aquí, v>Oes el coeficiente de viscosidad cinemática., y fE [H-I(O)]d es lIn término fuente dado. Consideraremos el problema estacionario con condiciones de Dirichlet homogéneas, con objeto de no introducir dificultades no esenciales en nuestro desarrollo. Estas ecuaciones constituyen una !inealización de las ecuaciones de Navier· 7 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 Stokes, que gobiernan alos flujos de Ruidos incompresibles en general ; Obtener u: O ---+ Rd, p; 51 __ Rtales que } u· \7u -v.6.u + \lp = f, \7. u=Oen O, u=O sobrer Analizaremos la aproximación de estas ecuaciones en la Sección 5_ Definamos la forma biJineal sobre [HJ(O>]d x{HJ(n)]cf, a(v, w) = (u· \7v, w) + v(\7v, \7w), Vv, w€ [HJ(O)]cf, (2) (3) (5) donde denotamos por (-,.) el producto escalar L2, para funciones sea escalares, vectoriales otensoriales. Supondremosque uE[V(O)]d con \1 -u = O. Entonces a(-,·) está bien definida, es oontinua y[HJ(n)]d-ooercitiva; osea, verifica a(v, w):5 eIvlllwh, a(v, v) 2: vlvl~ Vv, w E [HJ(O)]d (4) Aquí, e es una constante positiva, y 1- h denota la seminorma sobre [H1(n)]d, Ivl, = [In l'7vl' dX]'" Podemos dar al problema (1) una formulación variacional mixta (esto es, incluyendo la presión además de la velocidad), como sigue; donde Obtener (y,p) E Y x M tal que B(y,p;v,q) =< r,v >, V(v,q) Eyx M; }(6) y= [HJ(O)J', M ~ L'(O), yBes la forma bilineal sobre el espacio producto (Y x M) x (Y x M) definida PO' B(y,p; v,q) = a(y, v) + b(v,p) +b(u,q). (7) A su vez, b(·,·) es la forma bilineal sobre YxM dada por b{v, q) = -(\7 -v, q). Por \iltimo, (.,.) representa la dualidad Y'- Y. La existencia yunicidad de soluciones del problema es consecuencia del siguiente resultado general, que enunciamos con evidente abuso de notación: Consideremos dos espacios de Hilbert Y y M. Sea b(·,·) una forma bilineal continua sobre Yx M, ydefinamos los operadores R:Y --t M', y su adjunto R· :M --t Y' por (Rv,q)~(R'q,v)=b(v,q), VVEY,VqEM. (8) Sea además a( -, -) una forma. bilineal continua sobre Yx Y, que supondremos coercitiva. Se tiene entonces 8 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 Teorema 2.1 Si 1m Res armdo, entonces el problema (6) posu una solución. Además, esta solución es únial sobre Yx M, donde id es el espacio cociente M{KeT R'. La prueba de este teorema está basada esencialmente en el Teorema del Rango Cerrado de Banach. (Ver Girault -RAviart (1986». Notemos que en el caso de las ecuaciones de Stokes generalizadas (1), Ker R' está. formado por las funciones constantes, de modo que A su vez, ImR es cerrado, como consecuencia del resultado que enunciamos a continuación : Teorema 2.2 Las siguientes proposiciones son equivalentes: i) 1m R es cerrado en M'. ii) El par de espacios (Y, M) sawface la condición inf· 81lp continua: Existe una constante {J >Otal que . b(v,q) mf sup ~ {J. ,EM VEY II v llyllqllA1 (9) (10) iü) La forma bilineal Btkfinida por (7) sobre YxMsatisface la siguiente condici6n inf -sup :Existe una constante ..., >Otal que . r B(U,T¡V,q) In _ sup > 'Y. {U,r)EYxM (V.q)EYXÜ lI(u, r)lIyx},f lI(v,q)lh'xM La condición inf-sup (9) es satisfecha por el par de espacios ([HJ(O)]d ,L2 (0)), como es bien conocido (Ver Girault yRaviart (1986)). Por otra parte, la condición (10) asegura la coercitividad de la forma Bsobre el espacio producto yx M. Esta forma B esta ligada ala formulación como punto de silla del problema de Stokes: Si la forma aes simétrica, la solución (u,p) es también la. solución del problema de optimización 'nr sup (J(v,q)), J(v,q) ~ a(v,v) + b(v,q) -(r, v) (U) VEV qEM Las ecuaciones (6) son justamente las ecuaciones de Euler·Lagrange asociadas aest.e problema de punto silla.. 9 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 La aproximación clásicadel problema (6) porel métodode Galerkinse realiza como sigue: Sean Yh, Mh dos subespacios. de dimensión finita de YyM, respectivamente. Estaremos interesados en los casos en que ambos son espacios de Elementos Finitos. El parámetro hrepresentará el mayor diámetro de los elementos de la malla. Buscamos un par (Yh,Ph) EYh x Mh tal que Obtener (Yh,Ph) EYhxMh tal que B(Yh,Ph;vh,qh) =< f,vh >, V(Vh,qh) E Yh XMh; }(12) Sea Rhel operador de Yh a. valores sobre Mhdefinido por (8). Rhes UD operador de rango cerrado, al ser lineal entre espacios de dimensión finita. Según el Teorema 2.2, este problema tiene solución única sobre el espado cociente rVh = Mh/Ker R'h. Sin embargo, de (12) se deduce fácilmente que debe tenerse Yh E Ker Rh' Esta condición puede ser excesivamente restrictiva. si el espacio Yh no es suficientemente "rico" en grados de libertad respecto al espacio Mh. En Pironneau (1991) se muestran \'8Xios ejemplos en que Ker.Rh se reduce al cero. Es claro, pues, que los espacios Yhy Mh deben ser compatibles de alguna forma para. poder garantizar la convergencia. de la discretizaciÓn. Una condición suficiente para. ello es la llamada "inf-sup" discreta, ode Brezzi - Babuska : Lema 2.1 Supongamos que lcu con.stantesl3h >Odefinidas por " . r b(Vh,'h} IJh = In sup q"UI" V"eY .. Il v hllyllqhlL .. ! están acotadas inferiormente por una constante {j' independiente de h. tDnces, existe una constante e>odependiendo de {j', lIallyllbll tal que IIY-Yhlly+IIP-PhIlM"SC{ 'nf lIy-vhll+ 'nf IIP-Qhll,,}, V"EY¡, qhEM" (13) En- (14) De este modo, el error en la aproximación considerada vendrá dado por el enor de interpolación sobre los espacios Yh yMh• La condición inf - sup (13) puede interpretarse en el sentido de que el operador nh admite una prolongación a todo Y, uniformemente continua en h. Para significar su importancia., escribiremos el problema (6) en forma matricial. Sean N'U y N" las dimensiones de YhyMh,respectivamente (Observemos que en el caso de las ecuaciones de Stokes, Ñfh =Mh () LW1)). Denotemos por {4>.Ii = 1"" ,N,,} y {,pi IJ = 1,'" ,N p} una base de Yh y otra de Mh, resp. 10 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 Definamos las matrices AyRde dimensiones NvxNvyNpxNv,resp., yel vector columna rde dimensión Nv,por (15) N. N., Supongamos que buscamos Yh y Ph en la forma Yh = LYi ¡Pi> PI,. = LPi tPj· 1_1 J"1 Entonces, el problema (6) equivale a! sistema linea! Ay + R'j; d, RY = O; (16) \ donde yypson los vectores de componentes y.¡ yPi' Puesto que la forma aes coercitiva, la matriz Aes definida positiva. Por ello, la primera. ecuación en (16) puede ser resuelta en y, de modo que (12) equivale al sistema lineal RA-1 Rtp =A-if. El número de condición de la matriz que aparece en este problema resulta ser Cond(RA-' R') = I~~I. De este modo, la. condición inf -sup discreta garantiza no solamente que este problema admite una única solución sobre MI,., sino que además el número de condición de la matriz está uniformemente acotado en h. La obtención de pares de espacios que satisfagan la condición inf - sup discreta ha sido objeto de estudio de bastantes investigadores desde 1975. Estalabor ha permitido descubrir una amplia gama de pares satisfaciendo esta condición. En el caso de presiones continuas, el elemento más "económico" (en el sentido de que el espacio de velocidades resulta de añadir un mínimo de grados de libertad al espacio de presiones) es el p¡- Burbuja - PI' también llamado minielemento. Con objeto de recordar la definición de este elemento, denotemos por TI,. auna triangulación de n. Suponemos que 7h está formado por triángulos si d=2 Ypor tetraedros si d -= 3. Denotemos por Ai K ), .. " A~~~ las coordenadas baricéntricas de un elemento I< E TI,.. Denotemos por bKel polinomio de grado d+1dado por bK= ).~K) . ).~K) .' .. >.~~l. Notemos que b1< E [HJ(]<)( Suponemos bKextendido por cero fuera de K. Consideraremos el siguiente espacio de elementos finitos "burbuja" sobre Th: Bh(O) = { ÍI :n -+ Rtal que til .... ESpan{bK}, VI< E TI,. }. Definamos por otra parte Vh(m) = {r EL2(O) tal que rl .... E P m, VI< ETh }, } Yh(m) ~ Wh(mll' n[HJ(n)]', Mh(m) =Vh(m). (17) 11 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 (18) Aquí, Pmdenota el espacio de polinomios sobre Rdde grado menor oigual que m. Entonces, el Mini-elemento es (YIl{l) EIl [B/¡(O)]d ,1111>(1»). Si se pretende trabajar con presiones discontinuas, el espacio de velocidad debe ser más rico en grados de libertad. En este caso, el par de espacios más económico conocido que satisface lacondición inf-sup es (YI>(d), MI> (O) n CO(ñ»). 3 Métodos de Estabilización Según lo dicho en la Sección anterior, la clave para conseguir estabilizar Elementos Finitos "incompresibles" ~s debilitar la condición de divergencia nula. ElJo es realizado en los métodos clásicos enriqueciendo el espacio de velocidades en grados de libertad frente al espacio de presiones. Otra forma de lograrlo es relajar explícitamente la condición \l .u= 0, cambiándola en 'V ' Uh =gil, siendo 9h una función "pequeña" conveniente. Un método de este tipo fue introducido en Brezzi yPitkfu:anta (1984), donde se considera la condición relajada 'V .UIl = a: h'2 t::.PI>' Con ello se consigue cierto control sobre el gradiente de presión (y no sólo sobre la presión), impidiéndose la aparición de modos de presión espúreos. Como contrapartida, al tratarse de un método de penalización, un error de consistencia (de orden hpara la norma L2de la pt'esión) es inevitable. Esta dificulta.d es salvada por los métodos de estabilización, basa.- dos esencialmente en formulaciones de tipo Lagrangiano aumentado. Cuando la forma a(·,·) es simétrica, toda solución del problema de punto silla (11) lo es también de 'nf sup{J(v,q)-QIIAu+l7q-fll_¡), vev qEM donde Av = -v~v ya es un coeficiente adeterminar. En Hughes, Franca y Balestra. (1986) se propone discretizar este último problema reemplazando la norma H- ldel operador por una aproximación: Se trata del conocido método "Galerkin -Least Squares" .. El coeficiente h1< actúa como un factor de escala que relaciona la norma L2 con la norma 11-1, elemento aelemento. Las ecuaciones de EuJer -Lagrange asociadas a(19) son Obtener (y",pl» E Yh xM" tal que BOLS(Yh,P,,; vh,ql» =FOLS(VI¡,q,,), \i(v",q},) E YI> XMh; 12 }(20) Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 f 2k = -o: Lhidf, 'V1j;¡JK> {(ET" siendo BGLS(y,p;v,q) = B(y,p;v,q)-a 2.: hid-Ay+\lp,A*v+\lq)J(; J(ETr, FGLS(V,q) = (f,v/-a 2.: hk(f,A*v+\lq)K' KET" Aquí, A* es el operador adjunto de A. Cuando el op€rador no es simétrico (caso u ,¡. O), el método GLS es también (20), con Av = u' 'Vv -vf:..v, (y A~v = -u· \Iv -v.6.v). Sin embargo, ahora ia solución (u,p) no es punto de silla de ningún funcionol. En el caso de elementos finitos afines atrozos, el método GLS se conoce como "Streamline ~ Upwind /Petrov -Galerkin" o, abreviadamente, SUPG. De hecho, este método fue desarrollado en primer lugar en Brooks yHughes (1982), ydio lugar posteriormente alos métodos GLS. Observemos que en el método SUPG se tiene .6. vh=Osobre cada demento K, con la consiguiente simplificación de la estructura de BCLS y FGLS Una interesante propiedad del método GLS es que es consistente, en el sen~ tida de que la solución exacta verifica el problema. discreto (20), al igual que el método de Galerkin. Sin embargo, es estable únicamente para volores del parámetro Qpequefios, por debajo de un valor crítico específico para cada Elemento Finito concreto, yque resulta difícil de calcular. Desde el punto de vista computacional, el método GLS presenta la ventaja de proporcionar un sistema lineol con matriz simétrica si el operador A lo es. En efecto, supo11garnos u= O. Definamos las matrices A, R, Ode dimensiones NvxNv,NpxNvyNpxNp,respectivamente, así como los vectores columna f\ y(2, de dimensiones NvyN rn respectivamente, por R..j = -(\7. (Pi, ?}Jj) -a v Lhk('V1/Jj, fl4J;)K, KeT,. Okj = -a 2.: fllCvt/!j, 'V1/Jdl\. J(ET" fH= {f,fú -QV L hJ«(f,flf¡)l\, KET" Entonces, el problema (20) es equivalente al sistema lineal (21) 13 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 Teorema 5.2 Supongamcs que 11~I-dv < 'l. Supcmgamc3 que u[H2(O)r':, P E [lIt (O)]. Entonces existe una constante dependiente de uypy tm ho >O tales que si O<h< 110, entonces Este teorema utiliza de forma esencial dos elementos: la regularidad de uyp, y ciertas estimaciones de normas L', 2 ~ q ~ 6, de funciones burbuja: IJVhllo,p ~ ehQI"hlt, VVh E(B,,(O)Jd, donde overifica Ot/2 + (1 - Ot)/6 = I/q. Queda aún la labor de determinar los coeficientes de estabilización. Se trata de un problema mayor en el desarrollo de estas técnicas, que áun está por resolver plenamente. La calidad de uno. solución numérica concreta asociada a una malla depende fuertemente de los parámetros de estabilización que hayan sido usados. En la práctica, se definen los 'Th,K de modo Que se tenga en cuenta la talla relativa local de los términ06 de convección y de difusión. 'Th,K es una función del número de Peclet local , PeK __ UKhK con UK=lIul!o,p;K' v Habitualmente, Th,K es de la forma si PeK ~ P, (37) si PeK > P; donde A es una constante, yP es un "escalón" predeterminado para el número de Péclet. Esto permite por una parte introducir una estabilización conveniente en zonas de grandes gradientes de velocidad debidos aconvección dominante (grandes PeK)' Por otra parte, se introducen niveles bajos de difusión numérica (de orden h'k) en regiones donde la difusión es dominante (bajos PeK)' Observemos que en la discretizaci6n (37) los coeficientes Th,K son funciones no lineales de la solución u". Si embargo, el análisis anterior se puede extender a este caso, ylos Teoremas 5.1 y5.2 siguen siendo válidos. Digamos, por último, que el análisis anterior puede e,.'Xtenderse aelementos finitos de órdenes superiores de interpolación, obien amétodos de tipo Douglas .Wang yotros, con ciertas adaptaciones técnicas. 20 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 En definitiva., la formulación (26) de los métodos de estabilización en térmi· nos de operadores "de condensación estática" tiene en cuenta de forma natural la regularidad H- 1del operador aestabilizar. Esta fonnuJación proporciona un marco matemático para el análisis numérico de estos métodos --entendiendo como tal el análisis de estabilidad, convergencia yestimaciones de error-, completamente análogo al que se tiene en el caso de los métodos clásicos. Ello permite incluso realizar tal análisis numérico en el caso de las ecuadones de Navier-Stokes, mejorando el análisis realizado hasta ahora con las formulaciones habituales. Referencias [11 I. BABUSKA (1973), The Finite Element Method with Lagrange multipliers. Numer. Math. 20, pp. 179-192. [21 F. BREZZI (1974), On tM eristence, uniqueness and approximation 01 saddle-point problems arising from Lagronge Multipliers. R.A.I.n.O. Anal. Numer., R2, pp. 129-151. [3] F. BREZZI, M. FORTIN (1991), Mi:ted and Hybrid Finite Elements. Springer·Verlag. [4J F. BREZZI, J. PITKARANTA (1984), On the stabilization 01 Finite Elements approximations of!he Stokes problem, Springer.Verlag, Berlín. In Efficient Solutions oC ElIiptic Systems, Notes on Numerical Fluid Mecltanics, Vol. 10 (W. Hnckb",h e<k), pp. 11-19. 15J A. N. BnooKS, T. J. R. RUGHES (1982), Streamline/Uptuind PetrotJ- Galerkin lormulations for convection dominated flow.s, with particular emphasi" on fll.e incompressib~ Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg., 32, pp. 199-259. [61 T. CHACÓN REBOLLO (1996), Aterm by tcrm stabilizati011 algorithm for Finite Element solution of incompressible fiow problems. Propuesto para publicación en Numer. Math. [7J J. DOUGLAS, J. WANG (1989), An absolutely stabilized Finite Elcnent method for the Stokes problem, Math. Comp. 52, pp. 459-508. [81 L.P. FRANCA, S. L. FnEY (1992), Stabüized FiniU Elements: l/. The in· compressible Navier·Stokes equotions, Comput. Methods Appl. Mech. En~ grg., 99, pp. 209-233. 21 Bol. Soc. Esp. Mat. Apl. nº9 (1996) 5-22 191 L.P. FRANCA, T. J. R. HUGHES, R. STENBERC (1993), Stabilized Finite E/ement Meth0d3. In "Incompressible Computational Fluid Dynamics" (iVI. D. Gunzburger and R. A. Nicolaides Eds.), Cambridge Univ. Press. [10] V. GIRAULT, P. A. RAVIART (1986), Finite Element MetJwd.s far NatJi.er· Stakes equanons. Springer-Verlag. (11} D. GRIFFITS (1979), The construction of approximately divergence-free fi· nite elements. In "Mathematics of Finite Elements and Applications" (R. Whiteman, ed.) Academk Press. [12} F. HECHT (1981), Construction d'une base de fonctions PJ non-conformes áditJef'9ence nuIIe dans Jil. R.A.I.R.O. Anal. Numer. 15, pp. 119-150. [13J T. J. R. BUCHES, L.? FRANCA, M. BALESTRA (1986), A new Finile Etement formulabon for CFD: V. Circumventing the Brez~BabtiSka rondition: A stabk Petrov-Gaierkin formulation of lhe Stakes problem accom· modating equal-arder interpolations. Comput. Methods Appl Mech. Engrg. 59, pp. 85-99. (14J T. J. R. HUCHES, L.? FRANCA, G. M. HUBERT (1989), A new Finite Element lormtdatilm lor CPD : VIll. The Galc-kin-least-squares method lor advedive·diffi¡sive equations. Comput. Methods Appl Mech. Engrg. 73, pp. 173-189. [151 J. L. LIONS (1969), Quelques Méthodes de Résolution des Pl'Oblémes a,llX Limites non Linéaires. Gauthier-VilIars, París. /161 PIRONNEAU, O. (1991), Finite EIement Methodslar Fluids Wiley -Masson. 1171 R. TEMAM (1977), Theory and Numerirol AnaJysis 01 tk NatJier-Stokes Equations. North-Holland, Amsterdam. 22