Aspectos matemáticos de algunos métodos numéricos en mecánica de fluidos: el problema incompresible de Navier-Stokes

Revista internacional de métodos numéricos pwa cálculo y diseño en ingenie&. Vol. 3,2,129-152 (1987) ASPECTOS MATEMATICOS DE ALGUNOS METODOS NUMERICOS EN MECANICA DE FLUIDOS: EL PROBLEMA INCOMPRESIBLE DE NAVIEEl-STOKES E. FERNANDEZ CARA Y M. MARIN BELTRAN Dpto. Ecuaciones Funcionales. Facultad de Matemáticas c/ Tarfia, s/n 41 01 2 Sevilla Cátedra de Matemáticas. Facultad de Ciencias c/ San Alberto, s/n 14071 Córdoba RESUMEN El objetivo de este estudio ha sido la obtención de propiedades de convergencia y estabilidad para dos esquemas numkricos que permiten resolver las ecuaciones incompresibles de Navier-Stokes. Dichos esquemas han sido obtenidos modificando ligeramente otros debidos a R. Glowinski, cuya convergencia no había sido estudiada hasta la fecha. En una primera etapa, se usan métodos de direcciones alternadas del tipo de Peaceman-Rachford y de Strang. Esto reduce el problema a la resolución de problemas elípticos lineales del tipo de Stokes y problemas elfpticos quasi-lineales. En la segunda etapa, estos problemas se resuelven numkricamente usando varios mktodos de aproximación en espacio (elementos finitos), (para los problemas no lineales es conveniente introducir una formulación de tipo minimos cuadrados). La convergencia de las soluciones aproximadas hacia la solución del problema inicial se verifica bajo ciertas condiciones especificas de estabilidad. Las propiedades obtenidas vienen a justificar los buenos resultados numkricos conseguidos utilizando los métodos de Glowinski. SUMMARY The goal of this paper is to describe some stability and convergence properties for two numerical schemes which can be used to solve the incompressible time-dependent Navier-Stokes equations. The schemes were derived by modifying slightly others, due to R. Glowinski, for which convergence had not still been proven. At a first stage, alternating direction time-discretization methods of PeacemanRachford and Strang types have been introduced. This reduces the task to the solution of a sequence of (stationary) eliiptic subproblems, some of them linear (quasi-Stokes problems) and some quasilinear. Then, .these are solved using FEMs for the nonlinear sub-problems, it is appropiate to introduce a least-squares reformulation). Under certain specific stability conditions, we establish a convergence result for the computed solutions. This justifies rigorously the fact that Glowinski's methods have provided excellent numerical results. Recibido: Diciembre 1986 QUniversitat Politkcnica de Catalunya (España) ISSN 0213-1315 130 E. FERNANDEZ CARA Y M. MARIN INTRODUCCION Las ecuaciones de Navier-Stokes, objeto de nuestro trabajo, tienen su origen en la Mecánica de Fluidos y juegan un papel muy importante en numerosos campos de la Ciencia y la Ingenieria. Por ello, no sólo nos interesa estudiar diferentes aspectos cualitativos sino tambien cuantitativos, de sus soluciones. Entre sus muchas aplicaciones, citaremos a modo de ejemplo, la determinación del campo de velocidades de un fluido viscoso que se mueve en el interior de un canal en el cual existe un escalón. (Fig. 1 a 3). Figura 1. La geometria del escalón. Las condiciones de contorno son: a) Fluido de Poiseuille sobre F..; b) Condición de no deslizamiento sobre Po; c) Tracción nula (condición natural) sobre F: , Figura 2. La triangulación utilizada (técnica MODULEF; elemento finito de Thomasset, P, -no conforme, base con divergencia nula). Figura 3. Visualización de algunas lfneas de corriente (cálculos realizados utilizando la Biblioteca MODULEF). Valor del número de Reynolds: 190. Vemos como se aprecian dos zonas bien diferenciadas; en la zona situada entre la entrada del canal y el escalón, el flujo transcurre en regimen laminar, debido a la presencia de una vorticidad prácticamente nula (el campo de velocidades posee perfil parabólico); existe una segunda zona en donde el fluido tiende a "rellenar el vacío" originado por el escalón, creándose un torbellino. Las dimensiones de dicho torbellino dependerán de la altura del escalón, la velocidad de llegada del fluido y la viscosidad del mismo. Nuestro objetivo va a ser pues, calcular el campo de velocidades de un fluido en movimiento. Denotemos por u=u(x,t) la velocidad del fluido en 52x(O,T), siendo 52 la región ocupada por el fluido (52 es un abierto acotado de Rn, con n= 2 ó 3 en la práctica) y [O,T] el intervalo de tiempo considerado. Sea p=p(x,t) la presión del fluido. Con estas notaciones, el movimiento de un fluido viscoso (de viscosiclad cinemática u), i.e. cuyas partículas están sometidas a fuerzas internas de rozamiento, incompresible y homogeneo (es decir, con densidad constante) queda modelizado por el problema incompresible de Navier-Stokes: div u = O en 52 x (0,T) (2) u = O sobre 852 x (0,T) (3) La primera de las ecuaciones, llamada ecuación de equilibrio de momentos o ley de conservación de la cantidad de movimiento, es una consecuencia de la segunda Ley de Newton. En ella, el primer y tercer termino provienen de derivar respecto de la variable t a lo largo de las trayectorias; el termino Vp se debe a la existencia de fuerzas superficiales de tensión, mientras que la presencia del termino vAu denota que nuestro fluido es viscoso, (si fuese v=O, se diria que nuestro fluido es ideal); por iiltimo, el segundo miembro, indica la existencia de un campo de fuerzas externas que actúa sobre el fluido. La ecuación (2) constituye la llamada condición de imcompresibilidad y, en nuestro caso, es consecuencia de la ley de conservación de la masa; nos dice que el volumen ocupado por el fluido permanece constante en todo 52 x (0,T). La condición de contorno (3) tiene un significado físico evidente: debido a la presencia de fuerzas viscosas, el fluido se adhiere a las paredes sólidas del dominio; si estas están en reposo, la velocidad sobre la frontera de nuestro dominio debe ser nula (recuérdese el comentario que haciamos sobre el perfil parabólico de velocidades en la entrada del canal). Naturalmente, si alguna de las paredes se mueve con una velocidad v, la condición anterior debe ser sustituida por la igualdad u=v sobre dicha pared. Por último, puesto que nuestro problema es evolutivo, necesitamos de una condici6n inicial de tipo de (4), donde u, es un campo de velocidades inicial conocido. RESULTADOS DE EXISTENCIA Y UNICIDAD Con objeto de enunciar resultados de existencia y unicidad para este problema, vamos a introducir una formulación debil equivalente. Para ello, utilizaremos los siguientes espacios funcionales: EL PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 131 132 E. FERNANDEZ CARA Y M. MARIN b) H1 (a) = 1 u/ v E L (a), Vv E L' (ay 1 , es un espacio de Hilbert para la norma c) H: (a) = la adherencia de 3(a) en H' (a); en H: (a), la seminorma es de hecho una norma equivalente a la norma de H1 (a). e) V = la adherencia de?)' en H: (a)" ; V está dotado con el producto escalar y la norma de H: (a)", que llamaremos respectivamente (1. , . » y 11.11. f) H = la clausura de?)' en L2 (a)"; H está dotado con el producto escalar y la norma de L2 (a)", que llamaremos respectivamente ( . , . ) y l. l. para todo u,v, w E H1 Se verifican las siguientes propiedades para V y H: V = {v/v e ~,'(a)",div v = O en H = )v/v E ~~(a)~,div v = O en a, v.n = O sobreaal V G H H'c, V' con inyecciones compactas y densas gracias a que a es acotado; por esto, el producto de dualidad <. , .>,x,, puede considerarse como una extensi6n , del producto escalar en H. EL PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 133 Por otra parte, usando el Teorema de De Rham, obtenemos que V' puede identificarse con el espacio cociente V' E H-1 (St)" / VL~ (St) (donde H-' (St)" es el dual de H: (a)"), lo que nos permitirá el que, a partir de la formulación débil de nuestro problema, una vez conocida la velocidad, la presión quede determinada salvo una constante aditiva. Con todo, procediendo de la forma usual, es decir, multiplicando la ecuación (1) por una función test, integrando en St y aplicando las fórmulas de integración por partes, la formulación débil del problema de Navier-Stokes queda corno sigue: Dados u, e H y f e L~ (0,T; V') encontrar u e L2 (0,T; V) n L" (0,T; H) tal que La pertenencia de u a L" significa "a grosso modo" que la energía cinetica del fluido va a estar acotada en todo el intervalo (0,T); por otra parte, el que u pertenezca a L2 se traduce físicamente en el hecho de que la pérdida de dicha energía asociada a la presencia de fuerzas viscosas, es finita. Asimismo, conviene señalar que la condición inicial u(0) = u, tiene sentido, dentro de esta formulación débil, ya que, si u E L2 (0,T; V) n LOD (0,T;H) satisface (5), entonces, considerada como una función de [O,T] en V' , u es continua y (6) puede entenderse (e.g.) como una igualdad en V'. Respecto al problema anterior, tenemos un resültado, debido a Hopf, que nos dice que (5)-(6) posee solución, siendo esta única en el caso n= 2. Si, además, f ,u, y S2 son "suficientemente regulares", entonces la función generalizada (vectorial) (x,t) -f u(t) (x) es tal que, para alguna otra función p (escalar), el par (u, p) es una solución clásica de ( 1 )-(4). APROXIMACION NUMERICA Cuando se intenta resolver el problema (1)-(4), nos encontramos con varias dificultades importantes como son: - Se trata de un sistema (no una ecuación) en el que las incógnitas son las n componentes de la velocidad y la presión. - Las incógnitas no son independientes ya que (u,, u,, . . u, ) están ligadas por la condición de incompresibilidad. - El problema es no lineal, debido a la presencia del termino (u . V)u, que aparece en la ecuación de equilibrio de momentos. - Para n= 3 no hay resultados de unicidad. Se trata de un problema abierto desde 1938 en que fue advertido por Leray y puede que, para viscosidades pequeñas, existan soluciones "turbulentas" que presenten grandes oscilaciones. Todas estas dificultades hacen que las soluciones exactas sólo hayan podido obtenerse en algunos casos muy concretos y siempre estacionarios, como P.e. para el movimiento de un fluido viscoso engendrado por la rotación de un disco plano inmerso en el mismo y que gira alrededor de su eje (T. Kármán, 192 l), o el de un fluido que se mueve entre dos paredes planas que forman un ángulo entre ellas (fluido de Poiseuille) (G. Hamel, 19 16). Este hecho justifica plenamente la importancia de buscar algoritmos de resolución numerica, que nos permitan obtener soluciones aproximadas de nuestras ecuaciones. Esta aproximación se lleva a cabo en dos etapas: en una primera etapa se discretiza en tiempo y, posteriormente, los (sub) problemas elipticos que resultan son aproximados en espacio. Existen muchos esquemas clhsicos y que aparecen de una forma natural para llevar a cabo la discretización en tiempo. Expondremos a continuación dos a modo de ejemplo. Consideraremos para ello el intervalo [O,T] dividido en N subintervalos de amplitud k, Asociamos a k y a la función f las fl, . . . , fN, donde En cada esquema definimos recursivamente, para cada k, una familia de pares de funciones (u: ,PO,), . . . , (u! ,pk) donde (u,", p,") es una aproximaci6n de (u,p) en el instante m k, comenzando ambos con 1 ) Esquema totalmente explfcito Para m > O, obtenemos (u,"", pk") de (u:, p,") resolviendo el problema lineal: 134 E. FERNANDEZ CARA Y M. MARIN (*) Supondremos el paso k de discretización en tiempo constante, aunque un análisis similar puede hacerse utilizando paso variable. div u:" = O enS2 Um+l - k - O sobre Este esquema, una vez formulado en forma generalizada, nos permitirfa obtener u:" de u: sin más que despejar, si no fuera porque la condicibn de incompresibilidad lo dificulta ligeramente. Se trata pues de un problema sencillo de resolver. El inconveniente que tiene este esquema es, sin embargo, que hace falta tomar un paso k muy pequeño para obtener soluciones aproximadas estables, además, de una precisión moderada, de primer orden. 2) Esquema de Crank-Nicholson Para m > O, obtenemos (u:", p:'l) de (u: ,p: ) resolviendo el problema: v A[(u:+' +u: ) / 21 + (u: . EL PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 135 u:" = O sobre a!2 Este esquema sigue siendo lineal, ahora la precisión obtenida es de segundo orden, pero de nuevo es necesario tomar un paso k pequeño para que el esquema sea estable. Como alternativa a estos esquemas podemos pensar en otros que conserven la no linealidad del problema original, esperando de esta forma conseguir mejores aproximaciones sin necesidad de tener que utilizar un paso demasiado pequeño. Este es precisamente el objetivo de nuestro trabajo, en el que presentamos dos algoritmos que constituyen una modificación de otros, dados por R. Glowinski et al. en 1980, cuya convergencia no habfa sido estudiada hasta ahora. Ambos esquemas pertenecen a la clase de los metodos de direcciones alternadas, el primero de ellos del tipo de Peaceman-Rachford (de dos pasos intermedios) y el segundo del tipo de Strang (con tres pasos intermedios). Para describirlos, notemos en primer lugar que el problema inicial (1)-(4) equivale a este otro div u = O en ax(0,T) (13) u = O sobre aS2x(O,7') (1 4) 1 Lo que hemos hecho es introducir en la ecuación (12) el termino no linea12(div u)u que en el contexto de nuestro problema, no afecta para nada a su resolución, puesto que u verifica la condición de incompresibilidad. De esta forma, sin embargo, se consigue una mayor simetria en el termino no lineal que será de utilidad en las demostraciones de los resultados de convergencia. Consideremos de nuevo el intervalo de tiempo [O,T] dividido en N subintervalos de amplitud k, y un parámetro 0 perteneciente al intervalo (O, l), que nos servirá para dar peso a la condición de incompresibilidad y a la no linealidad. Para la descripci6n de los algoritmos, con objeto de simplificar la notación, suprimiremos el indice k. ALG 1 Describiremos primero un metodo de direcciones alternadas del tipo.de Peaceman-Rachford. Definimos las funciones u", U%, U', . . . , U~-%, U~ , como sigue: uO = u,; (16) despues, dado m O y um e H: (a)", calculamos la solución (um*',pm+%) de 1 - fm+i1 (um .v)um + - (div um )um 1 en a, 2 um+% = O sobre di2 y calculamos una solución um+ l de 1 + - (div um+l) (( 1-e)um+l + eum+%) = fm+l + vpm+% en a, (18) 2 div um+' = O en S2 Um + 1 = O sobre 136 E. FERNANDBZ CARA Y M. MARIN EL PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 137 En 4 17) y (1 8), f m+'h esta dada por i 2 (m+~)k fm+"(X) =- f(x,t) dt para i=1,2. k lm+z)k 2 El error de discretización local de este algoritmo es de primer orden en k. Desarrollamos un segundo algoritmo, del tipo de Strang, en el cual se introduce un paso intermedio m8s. ALG 2 Definimos ahora las funciones u", u% , U%, U', . . . , uN-%, uN, como sigue: uO = u,; (1 9) después, dados m 2 O y u" e H: (a)", calculamos la solución (um+%, pm+%) de div um+% = O en S2 um+% = O sobre a8 calculamos una soluciBn um+% de um'% = O sobre aS2 .y, de aquf, obtenemos: 3.a etapa: De las estimaciones de las etapas anteriores se deduce la existencia de subN N sucesiones u,. , . ,v, . , , , u,. , . y v , . , . y de una función u E Lm(O,T s2 (a)" ) fl LI (0,~; (52)") n HP(o,T; L2 (a)") tal que dichas subsucesiones convergen hacia u cuando k' y h' + O en los siguientes sentidos: debilmente - * en L"(0,T;L2 (a)") dCbilmente en L2 (0,~; H: (a)") debilmente en H@(O,T ;L2 (a)") Y fuertemente y c.p.d. en L2 (G?x(O,T))" De todas estas convergencias, quizás la menos evidente a partir de las estimaciones anteriores, sea la última, esta se deduce sin más que tener en cuenta las hipótesis de consistencia verificadas por nuestra aproximación espacial. Además, se obtiene que el limite es el mismo para todas ellas, gracias a las relaciones existentes entre cada dos parejas de familias obtenidas en la primera de las etapas. 4.a etapa: El próximo paso consiste en demostrar qbe el limite u es solución. Se deduce fácilmente de las etapas previas que este limite u satisface: (duldt) EL' (0,T; V') con o=4/3 si n=3 y a=2 si n=2. (por densidad, esta ecuación variacional tambien se verifica para todo v de V). Esta igualdad la hemos obtenido tomando limite en las igualdades verificadas por las funciones aproximadas, teniendo en cuenta las propiedades de convergencia de la 3.a etapa. Por último, no es difícil probar que u verifica tambien la condición inicial (6). etapa: Para demostrar la convergencia fuerte tenemos en cuenta que, en el caso n =2, se verifica la llamada igualdad de la energía: 144 E. FERNANDEZ CARA Y M. MARIN EL PROBLEMA INCOMPRESIBLE DE NAVIERSTOKES 145 que nos dice que la energia cinética correspondiente al instante final más la pérdida de dicha energia debido a la presencia de fuerzas de viscosidad, coincide con la energia cinética correspondiente al instante inicial más el trabajo realizado por las fuerzas exteriores. Teniendo en cuenta esta igualdad y, con objeto de probar la convergencia fuerte en L2 (0,T;q (a)"), introducimos la expresión: donde: * wkh : [O,T] -+ Wh es una función constante en cada [mk,(mtl)k) definida por wk, (mk) = BU;+% + (1-o)u;, * u; : [O,T] -+ Wh con u; E L2 (0,T; Wh) y tal que Se verifica que, dada u E L2 (0,T; V), existe siempre una funci6ri u; con las propiedades precedentes (cf. T6mamg). * or, es una sucesión de parámetros positivos que convergen a 1 cuando k y h tienden a cero.. * Ykh es una expresión convenientemente elegida, como veremos m;ls adelante. Si desarrollamos X,, , obtenemos: Xw = xkh + xkh + X; 3 donde x;, = Ilu(T)II:2 + 201, v loT lluh(t)llii dt O Basta, pues, elegir Y,, adecuadamente para que: Con ello conseguiremos que Xkh + O y, de aquf, " .< Teniendo en cuenta (46), concluimos que y, utilizando (40), llegaremos a las mismas propiedades de convergencia para ukh, Vkh 9 Ukh Y Vkh Observación 1. La condición k/h2 < e. es una condición suficiente de estabilidad, usual para problemas parabólicos no lineales. Nos dice, en esencia, que, una vez fijado un pequeño paso de discretización en espacio, el paso de discretización en tiempo no puede ser demasiado grande. Para otros esquemas de discretización en espacio, es decir para otros espacios W,, se obtienen condiciones de estabilidad similares, de la forma k/hq < cte. con q 2 1 Obsewación 2. En el caso tridirnensional no hay resultados de convergencia fuerte, en principio debido a que no se verifica la igualdad de la energía (que ha servido como base en nuestra demostración), Además, la no unicidad $e solución impide un resultado de convergencia para las sucesiones completas. 146 B. FERNANDEZ CARA Y M. MARIN l':L PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 147 Observación 3. Las tecnicas de demostración utilizadas están basadas en metodos de compacidad clásicos, y, en resultados de Témam (1977) para esquemas más simples. Para algunos resultados parciales, veanse los trabajos previos F. Cara - Marín3. 4* '. Observación 4. Se han obtenido resultados análogos para el segundo de los algoritmos con otra constante e, para la correspondiente condici6n de estabilidad; la demostraci6n7 algo más compleja puede encontrarse en Marin? Observación 5. Claramente, ALG 2 es, en la práctica, casi tan econ6mico de usar como ALG 1. Por otra parte, el estudio detallado de algoritmos análogos a ALG 1 y ALG 2 aplicados a un problema modelo sencillo (cf. ~risteau' ) parece indicar que ALG 2 está mucho mejor adaptado a la determinación de soluciones para grandes valores de t y, en particular, al cálculo de soluciones estacionarias. RESOLWCION EFECTIVA DE LOS PROBLEMAS Comentaremos a continuaci6n los metodos numericos de resolución de los problemas lineales y quasi-lineales obtenidos en cada paso de ALG 1 y ALG 2. Por simplicidad, nos limitaremos a trabajar sobre los problemas (23) y (24). El desarrollo que sigue puede adaptarse sin dificultad a los problemas discretizados correspondientes. Resolución de los problemas quasi-lineales El problema (24) no proviene de ningún problema del Cálculo de Variaciones, puesto que no existe ningún funcioisal sobre H: (a)" cuya diferencial coincida con 1 (v.D)v + - (div v)v 2 Sin embargo, usando una formulación conveniente en el sentido de los mfnimos cuadrados, vamos a ser capaces de resolver (24) mediante metodos propios de la Programación no lineal. Así, dada v E H: (a)", definimos y (=y (v)) E H: (a)" como la dnica soluci6n de y = O sobre aC2 Observamos que y se obtiene de v resolviendo n problemas de Dirichlet lineales inde- (pendientes (uno para cada componente de y). 148 E. FERNANDEZ CARA Y M. MARIN Observese que si v es solución del problema no lineal (24), entonces la correspondiente función y, obtenida a traves de (49), es O. Por ello, parece natural introducir la siguiente formulación de tipo minirnos cuadrados del problema (24): "Buscar u E HO (a)" tal que: (50) J(u) < J(v) V v e H: (a)" ". Aqui, J: Ii: (a)" + R está definida por 1 J(v) = - 1 p~~y~2 + hly121 dx, 2 (51) con y definida (a partir de v) resolviendo el problema lineal (49). Los problemas (50) y (24) son equivalentes en el sentido siguiente: Si u es solucidn de (24), también lo es de (SO), con J(u) = O; y, reciprocamente, si u es soluci6n de (50) y J(u) = 0, entonces es tambien soluci6n de (24). Para resolver el problema de minimizaci6n (50) usamos un algoritmo de tipo gradiente conjugado, definido como sigue Paso O Inicialización: u" E H: (a)" dado; a partir de u', calculamos gO e H: (a)" resolviendo <-pago + hgO ,Z> H-l ,H: = G'(u0),z> VZEH~(~)~, (*) w" = go (54) Despues, dados m > O, u", gm y wm , calculamos u"", gm" y wm+l en la forma siguiente : Paso 1 Descenso: (55) (* ) S(.) es el gradiente de J(.). y hacemos EL PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 149 donde X, es soluci6n del problema de minirnización unidimensional Paso 2 Construcción de una nueva dirección de descenso: Calculamos gm'l t. H: (a)", resolviendo En Marin7 puede encontrarse una expresión integral para <J'(um ),z>. Posteriormente, tomamos (variante de Polak-Ribiére): y, finalmente, wm+l = gm+l - ym wm. Cambiamos m por m+ 1 y volvemos a (55). Los dos pasos no triviales del algoritmo (52) - (58) son: (i) La resolución del problema de minirnización en unavariable (56); el correspondiente hm puede calcularse de forma aproximada (por un método de dicotomia o mediante el método de Fibonacci; cf. e.g. ~ao*) o bien de manera exacta, sin más que calcular las raíces de una ecuación algebráica de tercer grado (en efecto, la funci6n J(um hm ) es polinómica de cuarto grado). Por otra parte, como ya dijimos, en cada paso hemos de resolver n problemas de Dirichlet lineales asociados al operador elíptico (Ald.- pA) para, dado v, obtener la correspondiente función y (pi-ob. (49)). (ii) El cálculo de gm a partir de u" requiere la resolución de 2n problemas lineales de Dirichlet del tipo anterior. De ellos, n servirán para calcular <J'(um ),z> y los restantes n para obtener propiamente gm mediante (57). Observese que, para una división uniforme del intervalo [O,T], los parLimetrss A y p , son independientes de m. En consecuencia, la elección de una aproximación particular en las variables espaciales reduce los mencionados problemas de Dirichlet para (AId-pA) a sistemas lineales de ecuaciones donde la matriz de los coeficientes es siempre la misma. Parece, pues, adecuada una factorización de Cholesky para dicha matriz (al principio del prcgrama) y limitar los cálculos a la resolución de sistemas triangulares. Por el contrario, si el metodo de discretización en la variable t es de pasol7ariable, la resolución de los sistemas lineales que provienen de los problemas de Diriciilet para (Ald - pA) resulta costosa por metodos directos, (ahora p depende de m!). Una alternativa interesante consiste en utilizar algoritmos de tipo relajación. Resolución de los problemas de tipo quasi-Stokes 150 E. FERNANDICZ CARA Y M. MARIN Figura 5. Elemento Pl -no conforme, Po (Thomasset) La resolución del problema (23) puede llevarse a cabo de dos maneras. En el caso de que sea posible calcular una base explicita de Vh, basta resolver (29) por un rnetodo directo, ya que se trata de resolver un sistema algebráico cuya matriz es simetrica y l:L PROBLEMA INCOMPRESIBLE DE NAVIER-STOKES 15 1 definida positiva. Esto ocurre, P.e., cuando usamos el elemento finito no conforme de Thomasset (cf. la figura 5). En este caso, se tiene: Wh = \wh /wh r L~ (52)' ; wh 1, E Pl (7')" V TGh ;wh continua en baricentros de caras adyacentes; wh la a = O 1 , Crouzeix2, calculó una base explícita de Vh en el caso n =2. Desafortunadamente, nos encontramos con que no es fácil en general, hallar una base de Vh ; esto sugiere resolver el problema (23) mediante un método iterado, del tipo de Uzawa o Arrow-Hurwicz. Una ventaja de la utilización de un metodo de este tipo es la obtención de una aproximación pr+%en ALG 1, (resp. p,"+x, p,"+l en ALG 2) para1 la presión p en el tiempo (m+l/2)k (resp. (m+1/4)k y (m+l)k). Aqui, nos referiremos exclusivamente al ,algoritmo de Uzawa. Comenzaremos con una función dada Despues, conocido pm , definimos um y pm'l (m 2 0) mediante: X*um - p*Aum = g-vpm en 52, u" = O sobre a52, Se supone que cp > O es un escalar dado. Acerca de la convergencia del algoritmo (5 9) - (6 1 ), tenemos el resultado siguiente: Teorema Supongamos que 152 E. FERNANDEZ CARA Y M. MARIN Entonces, cualquiera que sea p" E L2 (a) se verifica: u" -+ U fuertemente en H: p" -+ p fuertemente en L2 (a), donde 1 u,pl es la única solucibn de (23) que verifica /a p dx = po dx. (Para la dernostrací6r-1, vease Glowinski" p. 303). En cada paso del algoritmo (59) - (61) deben resolverse, de nuevo, n problemas de Dirichlet Independientes para (X* Id - p* A). Si se toma 8 = 112 (resp. 8 = 213) en ALG 1 (resp. AL@ 2), se sigue que (XId -FA) = (A* Id -F* A) = eId -+A> (resp. 4 (Xld - pA) = 2 (A' Id - p* A) = (;I~d -$A), con lo cual muchos de los subprogramas pueden utilizarse tanto para los subproblemas ,lineales como para los no lineales, reduciendo asi memoria y tiempo de programacibn, Los resultados expuestos s0.n generalizables a1 caso de condiciones de contorno de tipo Dirichlet no hornogeneas (paredes deslizantes, entrada de flujo, ets.) y10 condiciones naturales hornogeneas o no (salida de flujo, tracción, etc.). Para detalles, cf. ~lowinski~, ~arin~. REFERENCIAS 1. M. 0~ Bristeau y Y. Periaeix. ""Fite Element methods for the calculatisn of compressible viscous flows using self-adaptive mesh refinements". Rapport de Recherche INRIA, (aparecerá). 2. M. Crsuzeix y P. A. Raviart. "Conforming and non conformirng finite element methods for solving the skationary Stskes equations". RAIRO, Serie Anal. NMm. 3,33-76, (1973). 3. E. Fernandez-Cara y M. Marfn. "The convergence of a numerical method of solution of tke NavierStokes equations". Internacional Syrnposium on Numerical Analysis, Madrid, Sep., (1 985). 4. E. FernBndez-Cara y M. Marfn. "AnBlisis numCrico de varios algoritmss para la resoluci6n de las ecuaciones de N-S". II Simposium sobre Aplicaciones delM.E.F., en Ingeniería. Barcelona, Junio, (1986). 5. E. Fernández-Cara y M. Marin. "The csnvergence of two numerical schemes for the N-S equations9'. Pendlcnte de publicación. 6. R. Glowinski. Nurnerial methods for nonlineal var.Pationa/probems. 2.a Ed. Springer-Verlag, (1 984). 7. M. Marin. Análisis numCrics del problema incompresible de N-S: métodos de tipo direcciones alternadas. 'Tesis, Universidad de Sevilla, Octubre, (1986). 8. SS Ras. Optimization 77zeory and applications. Wiley Eastern Limited, (1979). 9. R. Teman, Theoty aand numerical analysis of the N-S equations. 2: Ed. North-Holland, Amsterdam, (1977).