Repositorio Institucional de Documentos

Métodos numéricos aplicados a la mecánica de fluidos Alejandro López Nieto Trabajo de fin del grado de Matemáticas Universidad de Zaragoza Introducción Este trabajo es una introducción a las ecuaciones que modelan el comportamiento de los fluidos newtonianos incompresibles. A lo largo del mismo se irán introduciendo las herramientas para en última instancia poder resolver satisfactoriamente un problema de mecánica de fluidos y evaluar el resultado obtenido. Sobran motivaciones para estudiar este tema que, habiendo nacido en la física y la ingeniería requiere de las herramientas más potentes de las matemáticas para estudiarse en detalle. Las ecuaciones de Navier-Stokes han demostrado tener la importancia suficiente para que algunos organismos premien el estudio de la existencia y unicidad de soluciones con cuantiosos premios. No es de extrañar teniendo en cuenta que tienen aplicaciones en campos que van desde la aeronáutica hasta medicina, pasando por la producción de energías renovables. La gran complejidad de este campo, unida al interés de la industria han hecho que se creen grupos de estudio dedicado en exclusiva al tema en los que personal de todas las ramas del conocimiento une fuerzas. El primer capítulo del escrito pretende derivar las ecuaciones fundamentales de una manera que resulte asequible al lector. Una vez planteado el sistema de ecuaciones de Navier-Stokes, se introducen algunos conceptos del análisis adimensional para poder llegar a las ecuaciones de Stokes, una linealización no dependiente del tiempo del sistema original que resulta más asequible para el estudio. El segundo capítulo introduce la herramienta fundamental que utilizaremos a la hora de resolver la ecuación de Stokes. Se trata del método de elementos finitos, con una pequeña base de teoría de espacios de Hilbert y distribuciones, estaremos en disposición de hablar de la existencia y unicidad de soluciones de problemas elípticos y podremos al final introducir y hablar de la implementación del método de elementos finitos y de la convergencia del mismo. Por último el capítulo tercero habla del planteamiento de las ecuaciones de Stokes como un problema de punto silla, estudia la existencia y unicidad de soluciones del mismo y la convergencia de los métodos de elementos finitos utilizando herramientas ya presentadas en el tema anterior. Al final se muestran los resultados obtenidos mediante computación para distintos espacios de elementos finitos. Todos los programas utilizados en este último punto serán además adjuntados en forma de anexo. III Índice general Introducción III Summary VII 1. Una introducción al modelo matemático 1 1.1. La derivación de las ecuaciones fundamentales . . . . . . . . . . . . . . . . . . . . 1 1.1.1. La conservación de la masa . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1.2. Los fluidos incompresibles . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.1.3. La ecuación del movimiento . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2. El número de Reynolds y la ecuación de Stokes . . . . . . . . . . . . . . . . . . . . 5 2. El método de los elementos finitos 7 2.1. La necesidad de nuevas soluciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2. Laderivadadébil .................................... 8 2.3. La formulación débil del problema . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.4. Existencia y unicidad de solución . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.5. La discretización del problema variacional . . . . . . . . . . . . . . . . . . . . . . . 13 2.6. Loselementosfinitos.................................. 14 3. Elementos finitos en la ecuación de Stokes 19 3.1. La formulación débil del problema de Stokes . . . . . . . . . . . . . . . . . . . . . 19 3.2. Existencia y unicidad de solución del problema variacional . . . . . . . . . . . . . . 20 3.3. Problema de Galerkin para la ecuación de Stokes . . . . . . . . . . . . . . . . . . . 23 3.4. Unproblemapráctico.................................. 26 3.5. La implementación del método . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 A. Programas utilizados 33 Bibliografía 41 V VI ÍNDICE GENERAL Elementos finitos en la mecánica de fluidos Summary Abstract The main goal of this text is to develop the ideas needed to solve the Stokes equation. To do so, we will first study the derivation of the main equations describing the motion of fluids, these are the well-known Navier Stokes equations. From there we will move on and obtain the Stokes problem that we want to study. Afterwards we will start to introduce the fundamental theory of finite element methods, these include the definition of weak derivatives, the concept of Sobolev spaces, the variational formulation of problems and some different types of domain decompositions. Once all these tools have been explained we will proceed to apply them to the Stokes problem and study its variational formulation, the well-posedness of the problem and the main ideas of mixed finite element methods. In the end we will be prepared to solve a specific real life problem using the MINI-element method, the results will be compared with those obtained using an unstable method. The equations governing the motion of fluids From now on we will consider Ωto be an open, bounded, connected subset of Rn. The derivation of the Navier-Stokes equation system is based in the application of two basic physical laws: 1. The conservation of mass. 2. The conservation of momentum. After some integration on manifolds we prove that the conservation of mass is equivalent to the fact that this equation holds at every point of Ω: 1 ρ Dρ Dt +∇·u=0. Where ρis the density function of the fluid, uis the speed vector and D Dt =∂ ∂t+u·∇denotes the material derivative operator. This equation is normally called equation of mass. From the definition of an incompressible fluid we get that it turns into ∇·u=0 on Ω.(1) Navier-Stokes equations model the behaviour of incompressible newtonian fluids, thus this result holds. We already have one of the equations of our system. The conservation of momentum comes from the expression Zτ Du Dt ρdτ=FV+FS.(2) Where FVare the long-range forces on the fluid and FSare the surface forces. Long-range forces are those like gravity or electromagnetic forces, they can be considered as constants on the fluid since VII VIII Capítulo 0. Summary they do not change a lot with distance. Surface forces appear when molecules interact with each other, they only affect neighbouring particles of fluid and are a much more delicate case. After applying Newton’s hypothesis we derive that for newtonian fluids equation (2) turns into: Du Dt =F−1 ρ∇p+ν∆u.(3) Where Fis the long-range forces acceleration, pis the pressure and νis the cinematic viscosity of the fluid. Equation (3) is usually called equation of motion together with (1) they form the Navier-Stokes equation system:    Du Dt =F−1 ρ∇p+ν∆u, ∇·u=0. (4) Navier-Stokes equations describe the behaviour of some of the most interesting fluids such as water and air. However, they are time dependent and not linear, because of this studying them is extremely complex. We will restrict ourselves to the case of slow flows in small spaces for fluids with a relatively high cinematic viscosity. In such flows viscous forces are much stronger than the inertial ones and the convective term of the material derivative u·∇ucan be neglected. Definition 1. The Reynolds number is the index Re =Inertial forces Viscous forces. Flows with a very small Reynolds number are called Stokes flows. The Navier-Stokes equations for a Stokes flow can be reduced to: −∆u+∇p=fon Ω, ∇·u=0,(5) this system of equations together with appropriate boundary conditions is the one we want to study. The Stokes equation system models the behaviour of several flows such as the movement of microorganisms swimming in a fluid or a lava flow. The weak derivative In this section we will deal with Dirichlet boundary value problems and apply the finite element method to some of them. First of all we have to define the existence of solutions in the weak sense. Classical solutions of partial differential equations problems have strong requirements concerning the softness of the functions. For example, let us consider the following problem −u00(x) = f(x)en (0,1), u(0) = u(1) = 0.(6) A classical solution must fulfil u∈C2(Ω)∩C0Ω. However, the problem stated in (6) models a real life phenomenon. It describes the deformation of an elastic string under a distribution of weights. In case we choose such a distribution to be discrete, physical intuition tells us that the solution does not even have to belong to C1(Ω). Thus the requirements on solutions must be lowered. We introduce the following concepts: Elementos finitos en la mecánica de fluidos IX Definition 2. Given an open set Ω⊂Rn, we define D(Ω)as the space of infinitely differentiable functions from Ωinto Rsuch that their support is compact, i.e. D(Ω) = {f∈C∞(Ω):supp(f)is compact}, where the support of f is defined as supp(f) = {x∈Ω:f(x)6=0}. Definition 3. The dual space of D(Ω)is the space of distributions D(Ω)∗. Any element in D(Ω)∗is called a distribution. Definition 4. Given T ∈D(Ω)∗and f ∈D(Ω), we will write T(f) = hT,fi, and define the derivative of T with respect to xias h∂T ∂xi ,fi=−hT,∂f ∂xii. Iterating this definition, we can define higher order derivatives, given α∈Nna multi-integer . We can define the α-order derivative of T as hDα(T),fi= (−1)|α|hT,Dαfi. For every f ∈D(Ω). We can prove that for every f∈L2(Ω)there exists a unique Tf∈D(Ω)∗such that hTf,gi= RΩfg dΩ. Thus we can introduce the weak derivative of order αof a function in L2as the function Dαfsuch that: hDαf,gi=hf,Dαgifor every g∈D(Ω). Definition 5. A Sobolev space of order k ≥0is a space of functions in L2(Ω)such that their derivatives of up to kth order lie in L2(Ω)as well: Hk(Ω) = {f∈L2(Ω):Dαf∈L2(Ω)for all α:|α|≤k}. Immediately from the definition we get that Hk+1(Ω)⊂Hk(Ω). Lemma 6. Consider the inner product on Hk(Ω)defined by: (u,v)Hk(Ω)=∑ |α|≤kZΩ(Dαu)(Dαv)dΩ, the pair (Hk(Ω),(•,•)Hk(Ω))is a Hilbert space. Definition 7. H1 0(Ω) = D(Ω) = f∈H1(Ω):f=0on Γ. The weak formulation of a general problem The following result is called Green’s formula: Theorem 8. Given u,v∈H1(Ω), it holds: ZΩ ∂u ∂xi vdΩ=ZΓ uvnids−ZΩ u∂v ∂xidΩ, i =1,...,n, where nidenotes the ith component of the unit exterior normal vector on Γ. Autor: Alejandro López Nieto 2Capítulo 1. Una introducción al modelo matemático Como esta relación es cierta para cualquier volumen contenido en Ωpodemos afirmar que en general la ecuación ∂ρ ∂t+∇·(ρu) = 0 (1.2) se cumple para todo punto xy todo instante de tiempo t. Definición 1.1.4. Se llama derivada material de una función cualquiera definida en R3×Ral operador lineal D Dt =∂ ∂t+u·∇. La derivada material puede observarse como el índice de variación de una función en el tiempo, siguiendo la corriente. La ecuación (1.2) se llama ecuación de continuidad o de conservación de la masa y admite la expresión equivalente 1 ρ Dρ Dt +∇·u=0. 1.1.2. Los fluidos incompresibles Las ecuaciones de Navier-Stokes son conocidas por modelar el comportamiento de fluidos incompresibles, sin embargo, aún no sabemos qué quiere decir esto. Una manera intuitiva de ver cuando un fluido es incompresible es la siguiente. Sea τel volumen de un elemento de fluido y consideraremos el índice de expansión local r(x,t) = l´ ım τ→0 1 τ dτ dt=l´ ım τ→0 1 τZδτ ∇·udV≈∇·u Se entenderá entonces que un fluido es no compresible cuando este ratio del volumen instantáneo alrededor de todo punto tienda a 0, es decir, cuando se cumpla ∇·u=0. Esta no es una definición muy rigurosa, en general se dice que un fluido es incompresible cuando la densidad de un elemento de fluido no es dependiente de cambios en la presión. Se puede llegar a demostrar mediante un estudio más profundo que si asumimos que los cambios de temperatura en nuestro material son despreciables (lo cual se puede hacer en muchos casos), entonces la definición formal es equivalente a que Dρ Dt =∂ρ ∂t+u·∇ρ=0 Esto quiere decir que si vamos siguiendo una partícula cualquiera del fluido y estudiamos la densidad del material en un entorno a su alrededor, ésta será constante en el tiempo. En particular obtenemos de nuevo la condición ∇·u=0 sustituyendo en (1.2). Esta condición será la primera de nuestras ecuaciones. 1.1.3. La ecuación del movimiento La ley de conservación del momento (segunda ley de Newton) nos dará la segunda de nuestras ecuaciones, las fuerzas que actúan sobre una porción de nuestro fluido pueden separarse en general en 1. Fuerzas de volumen: se trata de las fuerzas que experimentan variaciones pequeñas con la distancia. Una consecuencia directa es que podemos considerar que una de estas fuerzas afecta por igual a todos los puntos de un elemento de fluido que estemos considerando, de modo que sobre dicho volumen la fuerza será directamente proporcional a la masa contenida por una aceleración constante. Ejemplos serían la gravedad o la fuerza electromagnética. Elementos finitos en la mecánica de fluidos 1.1. La derivación de las ecuaciones fundamentales 3 2. Fuerzas de superficie: se trata de las fuerzas con origen en la interacción molecular que decrecen muy rápidamente cuando aumenta la distancia. Se consideran despreciables a no ser que haya un contacto directo entre dos elementos. Si dos elementos de fluido están en contacto directo, la interacción se producirá sólo sobre las partículas que estén a una distancia corta de la frontera. Considerando elementos con caras planas y considerando la distancia de penetración de las fuerzas pequeña en comparación con la superficie de una cara elemento, la fuerza ejercida sobre dicha cara sera proporcional a la superficie de la misma Σ(n,x,t)δS, donde nes el vector normal unitario exterior de la cara que consideramos. Definición 1.1.5. A la magnitud fuerza por unidad de superficie Σ(n,x,t) = Fuerza ejercida sobre esta cara del elemento Área de dicha cara se la llama esfuerzo. La fuerza aplicada sobre nuestro elemento será de signo contrario a la normal exterior FS= −Σ(n,x,t)δS. Para determinar el esfuerzo sobre un elemento de superficie de un elemento de fluido que corresponde a un punto en un tiempo fijo definimos el tensor de tensiones σ∈Rd×d.σij será la componente idel esfuerzo sobre un elemento de superficie con vector normal en la dirección jdel espacio, así queda Σi(n) = σijnj. Donde njes la componente jdel vector normal unitario exterior con respecto a nuestro elemento de superficie δS. Podemos diferenciar los esfuerzos en función de la dirección de la normal de la superficie en la que actúan. Los esfuerzos en direcciones paralelas al sistema de referencia en el que trabajamos se llaman normales y sus valores se reflejan en la diagonal de la matriz σ. Los esfuerzos asociados con los valores de la matriz fuera de la diagonal se llaman tangenciales. Consideraremos ahora el momento sobre un elemento de fluido, se trata de la cantidad Zτ ρudτ. La fuerza ejercida sobre dicho elemento se define como la variación del momento en el tiempo. Con una serie de deducciones se llega a la expresión d dtZτ ρudτ=Zτ Du Dt ρdτ=FV+FS. Por la conservación del momento se tiene Fuerza ejercida =Fuerzas de volumen+Fuerzas de superficie, y se llega a que Zτ Du Dt ρdτ=FV+FS, con FV=Zτ FρdτyFS=Zτ ∇σdτ, desde el punto de vista de las componentes queda la expresión Zτ Du Dt ρdτ=Zτ Fρdτ+Zτ ∇·σdτ. La igualdad es cierta para cualquier elección del volumen, por tanto tiene que ser que en todo punto de Ωse cumple la ecuación ρDu Dt =ρF+∇·σ.(1.3) Autor: Alejandro López Nieto 4Capítulo 1. Una introducción al modelo matemático Esta ecuación es la que se suele denominar ecuación del movimiento. Sin embargo, para poder determinar el estado de un fluido a través de ella tenemos que entender mejor lo que son Fyσ. Frepresenta la aceleración inducida por las fuerzas que hemos llamado de volumen. En la mayoría de los casos se trata sencillamente de la aceleración gravitatoria, constante en todo el fluido y no hay mayor problema, nos abstendremos de evaluar casos más generales. El estudio del tensor de tensiones σplantea una mayor dificultad. En un fluido estático se tiene σ=−pI,donde Irepresenta el tensor identidad d×dcon dla dimensión del espacio y pes la presión. Sin embargo, cuando el fluido está en movimiento no es cierto que la presión actuará sobre todos los elementos como una constante en todas las direcciones. Podemos dividir el esfuerzo en dos tipos según su origen: 1. Esfuerzo interno: Debido a la presión del fluido, se manifiesta cuando éste se mueve al cambiar la presión respecto a la que había en el estado de reposo. El tensor asociado a este esfuerzo viene dado por σ1=−pI. 2. Esfuerzo viscoso: Se originan como respuesta a la deformación del fluido. Están estrechamente relacionados con la viscosidad µ, una magnitud que mide lo espesos que son los fluidos. La experiencia nos dice que cuesta más mover una masa de fluido más espeso, como puede ser la miel, que desplazar una masa de agua, esto va directamente ligado a la viscosidad. El origen de esta oposición del fluido al movimiento es puramente molecular, surge de las colisiones de moléculas de nuestro material que se producen cuando lo perturbamos. Se dice que un fluido es newtoniano cuando su esfuerzo viscoso satisface la siguiente relación (propuesta en su día por el propio Newton) σ2=µ(∇u+∇uT). Observación 1.1.6. A la forma ∇u+∇uTse la llama también tensor de deformación, así que la característica de los fluidos newtonianos es que el tensor esfuerzo viscoso es una función lineal del tensor deformación. Hay que decir que la familia de los fluidos newtonianos comprende a muchos de los fluidos habituales como el agua o el aire bajo condiciones de temperatura y presión normales. Agrupándolo todo tenemos una expresión del tensor de esfuerzo en un fluido newtoniano σ=σ1+σ2=−pI +µ(∇u+∇uT). Si introducimos todo esto en la ecuación (1.3), obtenemos ρDu Dt =ρF−∇p+µ(∇2u+∇(∇·u)). En forma vectorial y teniendo en cuenta ∇·u=0 se tiene que    Du Dt =F−1 ρ∇p+ν∆u, ∇·u=0. (1.4) Y así hemos llegado por fin al sistema de ecuaciones de Navier-Stokes (1.4) que modela el comportamiento de fluidos newtonianos incompresibles. Sin embargo el objetivo final de este trabajo es una linearización de este sistema que veremos en la sección siguiente. Observación 1.1.7. A la magnitud ν=µ ρse la llama viscosidad cinemática, en contraposición a µ que suele ser la viscosidad dinámica o absoluta. Todas estas variables dependen del estado del fluido, entre otras cosas se ven afectadas por las temperatura pero en muchos casos se pueden considerar como constantes. Elementos finitos en la mecánica de fluidos 1.2. El número de Reynolds y la ecuación de Stokes 5 1.2. El número de Reynolds y la ecuación de Stokes Hay que distinguir que podemos hablar de flujos y fluidos. Un fluido es el material concreto con el que estemos trabajando sea agua, aire o gel. El concepto de flujo es el de un fluido junto con las condiciones en que se encuentra, por ello podemos clasificar los flujos utilizando la siguiente definición. Definición 1.2.1. Se llama número de Reynolds a la relación Re =LU ν.(1.5) A un flujo con un número de Reynolds muy pequeño (cercano a cero) se le llama flujo de Stokes. En (1.5) LyUson una longitud representativa del problema planteado (por ejemplo el diámetro del dominio) y una velocidad representativa (podría ser una velocidad constante en la frontera), con esto podemos permitirnos definir nuevas variables u0=u Ux0=x Lt0=tU Lp0=p−p0 ρU2donde p0es una presión representativa para el problema Con estas nuevas variables podemos escribir (1.4) en forma adimensional. Dos flujos con las mismas condiciones iniciales y de contorno en forma adimensional y mismo número de Reynolds se llamarán dinámicamente similares y admitirán una misma solución en forma adimensional. Esto explica por qué el movimiento de sangre en una arteria principal y en un capilar siguen ecuaciones tan distintas y da una caracterización de los flujos de Stokes. En general los flujos con números de Reynolds bajos son aquéllos tales que la velocidad y longitud características son pequeñas en comparación con las viscosidad cinemática. Además de la definición (1.5), el número de Reynolds admite también la expresión Re =Fuerzas inerciales Fuerzas viscosas .(1.6) Hay que notar las consecuencias directas de (1.6) sobre los flujos de Stokes. En un flujo de este tipo las fuerzas viscosas dominarán a las inerciales. A causa de esto el término convectivo u·∇u que aparece en la derivada material de (1.3) se puede despreciar. Como resultado nuestro sistema de ecuaciones ha pasado a ser lineal y además de esto, por tratarse de un flujo lento, podemos considerar que será constante en el tiempo y conformarnos con resolver el sistema de ecuaciones estacionario (ajustando las constantes como nos convenga) −∆u+∇p=F, ∇·u=0.(1.7) Equipado con las condiciones de contorno que adecuadas. Hemos llegado a nuestra meta, el sistema de ecuaciones de Stokes. Aunque por norma general este sistema es interesante por ser más sencillo de resolver que (1.4), de hecho resulta que sí que modela algunos ejemplos concretos como pueden ser microorganismos nadando, un flujo de lava, algunos polímeros viscosos o el comportamiento de la sangre en los capilares. Autor: Alejandro López Nieto Capítulo 2 El método de los elementos finitos Antes de nada hay que decir que con el objetivo de simplificar todo un poco y aliviar la carga de teoría, en adelante trabajaremos con problemas de ecuaciones en derivadas parciales con condiciones de contorno de tipo Dirichlet. Todo lo que está por desarrollarse a partir de ahora admite una extensión relativamente sencilla a condiciones de contorno de tipo Neumann tan sólo con un poco de trabajo y notación extra. El libro [B] es una excelente guía a este efecto. 2.1. La necesidad de nuevas soluciones Dado un problema de ecuaciones en derivadas parciales, se tiene por norma general la formulación fuerte del mismo. Se conoce como solución clásica a aquélla que es continua en la frontera de un abierto Ω⊂Rny con derivadas continuas del orden de la ecuación en el interior del abierto, por ejemplo: −∆u=f,en Ω, u=0,en Γ=∂Ω.(2.1) Se trata de una ecuación en derivadas parciales de segundo orden, por tanto a la solución clásica u de dicha ecuación se le exigirá u∈C2(Ω)∩C0Ω, se puede ver que existen problemas físicos reales descritos por esta ecuación que admiten soluciones más relajadas. Si consideramos el comportamiento de una cuerda elástica de extremos fijos al ponerle pesos, el fenómeno físico viene descrito por el siguiente problema de valor en el contorno: −u00(x) = f(x)en (0,1), u(0) = u(1) = 0.(2.2) Se trata de una ecuación diferencial de segundo orden y una solución clásica tendrá que ser C2en (0,1)y continua en [0,1]. Si probamos a poner dos pesos en los puntos 0,4y0,6 la intuición física (y la realidad) nos da la solución mostrada en la Figura (2.1). Si esta vez ponemos un único peso en el punto intermedio, la deformación de nuestra cuerda será algo como lo reprensentado en la Figura (2.2). Sin embargo, claramente, estas soluciones no son C2en el interior. Nuestros motivos para buscar nuevas soluciones están ahora más que justificados, y por tanto debemos plantear la formulación débil (variacional) del problema. 7 8Capítulo 2. El método de los elementos finitos Figura 2.1: Deformación de una cuerda elástica por dos masas puntuales. Figura 2.2: Deformación de una cuerda elástica por una única masa puntual. 2.2. La derivada débil Definición 2.2.1. Sea Ω⊂Rnabierto, definimos D(Ω)como el espacio de funciones infinito diferenciables de Ωen Rcon soporte compacto, i.e. D(Ω) = {f∈C∞(Ω):supp(f)es compacto} donde supp(f) = {x∈Ω:f(x)6=0}. Definición 2.2.2. Al espacio dual de D(Ω)se le llama espacio de distribuciones D(Ω)∗. Una vez introducido esto, podemos hablar de derivación en sentido débil. Definición 2.2.3. Sea T ∈D(Ω)∗y f ∈D(Ω), escribimos T(f) = hT,fi, y definimos la derivada de T respecto de la variable xicomo h∂T ∂xi ,fi=−hT,∂f ∂xii. Iterando esta definición, sea α∈Nnun multiíndice1. Podemos definir la derivada de orden αde T hDα(T),fi= (−1)|α|hT,Dαfi. Donde se entiende que si f ∈D(Ω) Dαf=∂|α| ∂xα1 1. . . ∂xαn nf con derivadas en el sentido habitual. 1Cada componente del vector indica el orden de derivación respecto de la variable correspondiente, además |α|= n ∑ i=1 αi . Elementos finitos en la mecánica de fluidos 2.2. La derivada débil 9 Notemos que la derivada débil está bien definida dado que f es infinito derivable y las derivadas parciales conmutan. Está claro que D(Ω)⊂L2(Ω), puesto que si aplicamos las propiedades del soporte compacto: ZΩ|f|2dΩ≤Mfkfk∞<∞∀f∈D(Ω). Existe una relación interesante entre T∈D(Ω)∗y L2(Ω), dado que si g∈D(Ω)yf∈L2(Ω), se tiene que ZΩ|fg|dΩ≤kfk2kgk2<∞ por la desigualdad de Hölder y es lineal por la linealidad de la integral, en otras palabras, la expresión anterior define un funcional lineal en D(Ω)y L2(Ω)puede observarse como un subespacio de D(Ω)∗ como sigue. Para todo f∈L2(Ω)podemos encontrar una distribución Tftal que hTf,gi=ZΩ fg dΩ∀g∈D(Ω). Además la aplicación que lleva la función fal funcional correspondiente Tfes inyectiva. Esto es una consecuencia del siguiente resultado de análisis funcional: Teorema 2.2.4. Sea Ωabierto en Rn, entonces D(Ω)es denso en Lp(Ω)1≤p<∞. Demostración. Se dejará sin demostración, ésta se halla en la página 109 de [Bre]. Entonces si tenemos f,f∗∈L2(Ω)tales que RΩ(f−f∗)gdΩ=0 para todo g∈D(Ω), se sigue f−f∗∈D(Ω)⊥. Considerando L2con el producto interno (f,g)L2(Ω)=RΩfg dΩ.L2,(•,•)L2(Ω)es un espacio de Hilbert y se puede descomponer en un subespacio cerrado y su complemento ortogonal: L2(Ω) = D(Ω)+D(Ω)⊥=D(Ω). Por tanto D(Ω)⊥=0 y se tiene que efectivamente f=f∗. Ahora podemos entonces hablar de derivadas débiles para funciones en el espacio L2(Ω). Por abuso de notación escribimos también Dαfpara la derivada débil de orden α∈Rnde f∈L2(Ω), es decir, la función tal que hf,Dαgi= (−1)|α|hDαf,gipara todo g∈D(Ω). Donde Dαgestá tomada en el sentido habitual. Definición 2.2.5. Se llama espacios de Sobolev de orden k a los espacios de funciones L2(Ω)tales que sus derivadas en sentido débil de hasta orden k se hallan también en L2(Ω): Hk(Ω) = {f∈L2(Ω):Dαf∈L2(Ω)para todo α:|α|≤k}. De la definición se sigue inmediatamente que Hk+1(Ω)⊂Hk(Ω). Podemos generalizar estos espacios para Lp(Ω)con p6=2 pero no será necesario para los propósitos de este trabajo. Observación 2.2.6. Hemos observado ya que D(Ω)⊂L2(Ω). Por otro lado las funciones en D(Ω) son infinito diferenciables en el sentido habitual, en particular las derivadas habituales cumplen la condición de las derivadas en el sentido débil y pertenecen a D(Ω), por tanto D(Ω)⊂H1(Ω). Visto esto la siguiente definición tiene sentido. Definición 2.2.7. 2Definimos la clausura de D(Ω)en H1(Ω) H1 0(Ω) = D(Ω) = f∈H1(Ω):f=0en Γ. 2En el caso unidimensional podemos extender fpor continuidad y definir así el valor en la frontera. Si fes una función vectorial hay que considerar la condición f|Γ=0 en el sentido de trazas como se presenta en la página 325 de [Bre], no entraremos en más detalle. Autor: Alejandro López Nieto 10 Capítulo 2. El método de los elementos finitos 2.3. La formulación débil del problema Sea el problema unidimensional con el que hemos empezado este mismo capítulo: −u00(x) = f(x), u(0) = u(1) = 0.(2.3) Lo pasaremos a formulación variacional. Para ello, multiplicamos la primera ecuación por una función test v∈V, espacio de funciones test que aun no vamos a especificar, e integramos la expresión sobre Ω. Como resultado obtenemos: −Z1 0u00vdx=Z1 0fv dx(2.4) si integramos el término izquierdo por partes −Z1 0u00vdx=−u0v1 0+Z1 0u0v0dx considerando funciones test tales que v(0) = v(1) = 0, la expresión que nos queda para resolver es Z1 0u0v0dx=Z1 0fv dx.(2.5) Para toda función test que valga 0 en los extremos, seguimos con la necesidad de definir cuales van a ser nuestros espacios de funciones test y soluciones. Si impusiésemos la condición de que la solución fuese C1estaríamos siendo muy exigentes puesto que en los ejemplos físicos vistos las derivadas primeras no son continuas. Aquí es donde emplearemos los espacios de Sobolev, por Hölder el término de la izquierda en (2.5) está bien definido para u0,v0funciones L2. Si u,v,fson L2el término de la derecha estará también bien definido. En resumen, la ecuación (2.5) tiene sentido para u,v∈H1 0(0,1). Observación 2.3.1. En particular las soluciones clásicas son soluciones en este nuevo sentido, hemos ampliado el espacio de soluciones y las gráficas (2.1) y (2.2) satisfacen sin problema las nuevas condiciones. Para dimensión mayor que uno, podemos plantear la fórmula de Green, ésta nos da una versión de la integración por partes en espacios de mayor dimensión. Teorema 2.3.2. Sean u,v∈H1(Ω), entonces se tiene: ZΩ ∂u ∂xi vdΩ=ZΓ uvnids−ZΩ u∂v ∂xidΩ, i =1,...,n, donde nies la i-ésima componente del vector normal exterior unitario en Γ. Se trata de una conclusión sencilla del teorema de la divergencia de Gauss. Aplicando esto por ejemplo a la ecuación de Poisson homogénea en dimensión nllegamos a: −∆u=f, u=0, en Γ, (−ZΩ v∆udΩ=ZΩ fv dΩ, u=0, en Γ, −ZΩ v∆udΩ=−ZΓ v∇u·nds+ZΩ ∇u·∇vdΩ=ZΩ ∇u·∇vdΩ. Obtenemos así la formulación débil del problema: Elementos finitos en la mecánica de fluidos 2.4. Existencia y unicidad de solución 11 (Encontrar u∈H1 0(Ω)tal que: ZΩ ∇u·∇vdΩ=ZΩ fv dΩpara todo v∈H1 0(Ω).(2.6) En el caso de que usea una función vectorial y las componentes de este vector sean independientes unas de otras, entonces resolver el problema es tan sencillo como hacerlo componente a componente. Hay que remarcar que esto no es algo que sea siempre posible. 2.4. Existencia y unicidad de solución Definición 2.4.1. Sea una forma bilineal a :V×V→Rcon V un espacio de Hilbert y l :V→Runa forma lineal, entonces se dice: 1. a(•,•)es coerciva o V−elíptica cuando ∃α>0tal que a(v,v)≥αkvk2 V∀v∈V. 2. a(•,•)es continua cuando ∃Ma>0tal que a(u,v)≤MakukVkvkV∀u,v∈V. 3. l(•)es continua cuando ∃K>0tal que l(v)≤KkvkV∀v∈V. A estas funciones nos referimos como funcionales lineales de V a lo largo del trabajo. Lema 2.4.2. Sea V un espacio de Hilbert, a :V×V→Rforma bilineal simétrica, coerciva y continua y l :V→R funcional lineal, entonces son equivalentes: 1. u ∈V es solución de: a(u,v) = l(v)para todo v ∈V. 2. u ∈V minimiza en V la expresión 1 2a(v,v)−l(v). Demostración. Sea u∈Vtal que a(u,v) = l(v)para todo v∈V. Entonces, para todo v∈V 1 2a(u+tv,u+tv)−l(u+tv)−1 2a(u,u)−l(u) = 1 2t2a(v,v)>0 si v,tson distintas de 0. Y por tanto uminimiza la expresión. Veamos el recíproco, sea u∈Vminimizando 1 2a(v,v)−l(v), entonces la función Fv(t) = 1 2a(u+tv,u+tv)−l(u+tv)tiene un mínimo en 0 para todo v∈V, necesariamente su derivada en el origen será 0 por ser éste un punto crítico. Fv(t+h)−Fv(t) h=1 2ha(v,v)+ a(u+tv,v)+l(v). Tomando h→0 y t=0 nos queda F0 v(0) = a(u,v)−l(v) = 0 para todo v∈V Teorema 2.4.3. SeaV un espacio de Hilbert real con producto interno asociado h•,•iV, a :V×V→R forma bilineal coerciva y continua y l :V→R una forma lineal continua. Entonces el problema Encontrar u ∈V tal que: a(u,v) = l(v)∀v∈V,(2.7) tiene solución única. Este resultado se conoce como teorema de Lax-Milgram. Un problema de este tipo se llama problema elíptico. Autor: Alejandro López Nieto 18 Capítulo 2. El método de los elementos finitos Construyendo las matrices Mxx =mxx ij =Zˆ K ∂ˆ φj ∂ˆx ∂ˆ φi ∂ˆx,Mxy =mxy ij =Zˆ K ∂ˆ φj ∂ˆx ∂ˆ φi ∂ˆy,Myy =myy ij =Zˆ K ∂ˆ φj ∂ˆy ∂ˆ φi ∂ˆy, para i,j=1,2,3 y teniendo en cuenta que Mxy = (Myx)TyCes simétrica la expresión superior queda: aij =∑ K|det BK|hC11mxx αβ +C12 mxy αβ +mxy βα+C22myy αβ i. Podemos aplicar por otro lado un procedimiento similar para el término de la derecha en la ecuación de (2.6), de este modo calculamos bi=∑ KZKfφi=∑ K|det BK|Zˆ K ˆ fˆ φα. Donde ˆ f(ˆx,ˆy) = f(FK(ˆx,ˆy)). Lo único que nos faltaría para acabar es introducir las condiciones de tipo Dirichlet en la frontera. Para ello, sea i∈Inodo Dirichlet, sustituiremos ai j por δi,j,j=1. . . , |V| ybipor la condición de contorno en el nodo ique en este caso particular es 0. Con esto se implementa fácilmente el método de los elementos finitos y ya tendríamos un algoritmo para en este caso resolver el problema (2.6) con elementos P1. Para ver la convergencia aplicaríamos nuevamente el lema de Céa (2.5.2) seguido de (2.6.3). Además la velocidad de convergencia será lineal con hen la norma H1. Se pueden definir elementos de muchos otros tipos y distintos espacios de funciones sobre ellos, sin embargo esto que hemos desarrollado servirá para los propósitos finales del trabajo. Con toda esta teoría que hemos venido estudiando hasta ahora hemos podido introducir las ideas fundamentales del método de elementos finitos y estamos ya en condiciones de resolver problemas elípticos de hasta dos dimensiones. Los casos de más dimensión aumentan la dificultad y complejidad pero se dejan intuir a través de las ideas introducidas. Sin embargo, nuestro objetivo es algo aún más complejo que esto. La ecuación de Stokes tiene por incógnitas dos funciones (velocidad del fluido y presión), esto empeora notoriamente la situación. Teorema 2.6.3. Sea el operador de interpolación Πr h:V→Vh, donde Vhes el espacio de elementos finitos Prsobre nuestra división regular6de Ωτh. Entonces, si r ≥1y m =0,1se tiene la acotación: |v−Πr hv|Hm(Ω)≤Chr+1−m|v|Hr+1(Ω)para todo v ∈Hr+1(Ω). Donde |·|Hkes la seminorma |f|Hk=s∑ |α|=kZΩ(Dαf)2dΩ, equivalente en Hk 0a la norma habitual. Esto asegura la convergencia del método para la mayoría de las funciones que podamos imaginar. Además da una estimación de la velocidad de convergencia. Tendremos convergencia con orden hr+1−m. Así pues cuanto más “grados de libertad” haya sobre el espacio de polinomios, mejor será nuestra convergencia. Demostración. Puede verse en la página 93 de [Q] 6No vamos a entrar en detalle respecto a lo que es un mallado regular, asumiremos que las divisiones del dominio que utilizamos son regulares. Para más información sobre este tema consultar la página 96 de [Q] Elementos finitos en la mecánica de fluidos Capítulo 3 Elementos finitos en la ecuación de Stokes 3.1. La formulación débil del problema de Stokes Planteamos el problema:    −∆u+∇p=fen Ω, ∇·u=0, u=u0en Γ, (3.1) donde esta vez Ωes un abierto de Rn,urepresenta la velocidad del fluido en un punto y pes la presión del fluido en un punto. Notemos que uestá escrito en negrita puesto que ahora la velocidad es una función vectorial. Resulta sencillo ver que a una solución clásica del problema se le exigirá (u,p)∈C2(Ω)∩C0Ωn×C1(Ω). Por lo tanto sobran los motivos para querer hallar una formulación débil del problema. Para ello empezaremos trabajando con un problema con condiciones de contorno homogéneas:    −∆u+∇p=fen Ω, ∇·u=0, u=0en Γ. (3.2) Consideraremos los espacios de funciones test V=H1 0(Ω)n,(3.3) M=q∈L2(Ω):ZΩ q=0=L2,0(Ω).(3.4) Aplicamos el procedimiento habitual a (3.1) multiplicando las correspondientes ecuaciones por funciones test v∈V,q∈Me integramos las dos ecuaciones. ZΩ v·∆udΩ+ZΩ v·∇pdΩ=ZΩ v·fdΩ,(3.5) ZΩ q∇·udΩ=0.(3.6) Aplicando la fórmula de Green obtenemos: ZΩ v∆udΩ=ZΩ ∇u:∇vdΩ−Z∂Ω(∇u·n)vds=ZΩ ∇u:∇vdΩ,(3.7) ZΩ v·∇pdΩ=−ZΩ p∇·vdΩ+Z∂Ω(v·n)pds=−ZΩ p∇·vdΩ.(3.8) Donde nes el vector normal unitario exterior y ∇u:∇v=∑ ij ∂ui ∂xj ∂vj ∂xi . 19 20 Capítulo 3. Elementos finitos en la ecuación de Stokes Una de las primeras cosas que notamos es que si (u,p)es una solución del problema, entonces (u,p+C)es también solución ∀C∈R. Para evitar este pequeño contratiempo normalizaremos pcon RΩpdΩ=0. Observación 3.1.1. Además, para que el problema esté bien planteado, de la ecuación de incompresibilidad (3.6) se deriva que para que la condición en la frontera Dirichlet sea compatible debe cumplirse: ZΩ ∇·udΩ=ZΓ u·nds =ZΓ u0·nds =0. Esto quiere decir que se preserva la cantidad de fluido dentro de la región estudiada. El balance entre fluido que entra y sale de la región debe ser nulo. 3.2. Existencia y unicidad de solución del problema variacional El problema variacional de Stokes tiene la forma:    Hallar (u,p)∈V×Mtales que: a(u,v)+b(v,p) = (f,v)L2(Ω)=l(v), b(u,q) = 0,para todo v∈Vyq∈M.(3.9) Un problema del tipo de (3.9) se llama problema de punto silla. Para esta familia de problemas se tiene el siguiente resultado: Lema 3.2.1. Sea el problema (3.9) con a :V×V→Rforma bilineal simétrica, coerciva y continua con V un espacio de Hilbert, l :V→Rlineal y continua y b :V×M→Rbilineal y continua. Sea F(v) = 1 2a(v,v)−l(v), entonces son equivalentes: 1. El problema del mínimo F(u) = m´ ın v∈VF(v)bajo la restricción b(v,q) = 0para todo q ∈M. 2. ues solución de (3.9). Demostración. Tomamos L(v,q) = F(v)+b(v,q). Solucionar el problema de minimización planteado es lo mismo que calcular la solución de m´ ın v∈Vm´ ax q∈ML(v,q).(3.10) Si denotamos Lv,q(t,h) = L(u+tv,p+hq)donde (u,p)es una solución de (3.10). Entonces (0,0) es un punto crítico de Lv,q(t,h)para todo v∈Vyq∈M. Y por tanto las derivadas parciales de esta función respecto a cada componente en (0,0)son 0. Lv,q ∂t(0,0) = l´ ım t→0 Lv,q(t,0)−Lv,q(0,0) t=l´ ım t→0 1 2(a(u+tv,a(u+tv)−a(u,u))−tl(v)+tb(v,p) t =a(u,v)−l(v)+b(v,p) = 0. (3.11) Lv,q ∂h(0,0) = l´ ım h→0 Lv,q(0,h)−Lv,q(0,0) h=l´ ım h→0 b(u,p+hq)−b(u,p) h=b(u,q) = 0.(3.12) Para todo v∈Vyq∈M. Esto que hemos obtenido no son otra cosa que las ecuaciones que aparecen en el problema (3.9). Para ver el recíproco, sea (u,v)una solución de (3.9). Entonces si tenemos que v∈K={v∈V:b(v,q) = 0 para todo q∈M}. Elementos finitos en la mecánica de fluidos 3.2. Existencia y unicidad de solución del problema variacional 21 se cumple que F(u+v)−F(u) = 1 2(a(u+v,u+v)−a(u,u))+ l(u)−l(u+v) =1 2a(v,v)≥0, y por tanto uminimiza a Fen K. Viendo el problema desde este punto de vista podemos aprovecharnos de 2.4.3 para demostrar la existencia y unicidad de solución de (3.9). Cabe remarcar que probar este resultado para la función u no tiene una gran dificultad, es en la otra variable de nuestro problema donde surgen los problemas como veremos ahora mismo. Observación 3.2.2. K={v∈V:b(v,q) = 0para todo q ∈M}es un espacio de Hilbert con el producto interior heredado (se ve que es cerrado, contiene a sus puntos de acumulación). De esta observación se desprende por (2.4.3) que a(u,v)−l(v) = 0 para todo v∈Ktiene solución única en K. Queda entonces ver bajo que condiciones existe una única p∈Mtal que: b(v,p) = l(v)−a(u,v)para todo v∈V.(3.13) En otras palabras, ver en qué casos existe una función p∈Mtal que el funcional lineal definido por l(•)−a(u,•)∈V∗puede ser expresado en la forma b(•,p) = l(•)−a(u,•). La existencia y unicidad de soluciones de (3.9) queda entonces ligada a la de este nuevo problema. No se trata de un tema para nada evidente, pero se puede deducir la existencia de una condición que caracterizará los casos en que (3.13) tiene solución única. Definición 3.2.3. Se dice que una forma bilineal b :V×M→Rsatisface la condición LBB cuando: ∃β>0tal que se cumple: ´ ınf q∈Msup v∈V b(v,q) kvkkqk≥β.(3.14) Se puede contemplar esta propiedad como una condición de compatibilidad necesaria entre los espacios de soluciones. La condición LBB impone automáticamente unicidad para la solución (si existe) del problema (3.13). En efecto, sea 0 6=p∈Mtal que b(v,p) = 0 para todo v∈V, entonces si se cumple la condición LBB debe ser p=0. De otra manera se tendría ´ ınf q∈Msup v∈V b(v,q) kvkkqk≤0. Por la linealidad en la segunda variable de la forma bsi hay solución será única. Falta por ver entonces que, efectivamente, la solución existe. Sea K0={f∈V∗tales que f(v) = 0 para todo v∈K}. Tenemos entonces el siguiente teorema. Teorema 3.2.4. Definimos el operador B:M→K0⊂V∗ p→b(•,p). Entonces son equivalentes: 1. b cumple la condición LBB. Autor: Alejandro López Nieto 22 Capítulo 3. Elementos finitos en la ecuación de Stokes 2. B es un isomorfismo lineal y se tiene kBpk≥βkpkpara todo p ∈M. Demostración. Asumiendo que bcumple la condición LBB, definimos el inverso de Ben B(M) B−1:B(M)→M. El operador Bes lineal y continuo por serlo también b, utilizando la condición LBB vemos que el inverso también es lineal continuo: kB(p)k=kb(•,p)k=sup v6=0kb(v,p)k kvk≥βkpk. Como B−1es continuo y M=B−1(B(M)) es cerrado, debe ser B(M)cerrado. Definimos el operador dual de Bcomo: B∗:V∗∗ →M∗ f→f(B(•)). Veremos ahora que (Ker B∗)0={g∈V∗:fg =0 para todo f∈Ker B∗}=K0. Aplicando el teorema de representación de Riesz dos veces, sea f∈Ker B∗yg∈V∗, entonces fg =uf,ugpara dos elementos uf,ug∈V, se tiene que f(B(p)) = uf,uB(p)=0 para todo p∈M y por tanto uf∈K. De hecho sea u∈K, el funcional definido por hu,ugipertenece claramente a Ker B∗. Sea g∈(Ker B∗)0, entonces fg =uf,ug=0 para todo f∈(Ker B∗) = Kpor lo que se cumple ug∈K⊥yg(v) = hv,ugi=0 si v∈K, es decir, g∈K0. Sea g∈V0, entonces ug∈K⊥y como uf∈Kpara todo fde (Ker B∗)0, está claro que fg = uf,ug=0. Queda por demostrar ahora que B(M) = B(M) = (Ker B∗)0=K0. Para ello tenemos el siguiente resultado: Lema 3.2.5. Sean E y F espacios de Banach y B ∈L(E,F)son equivalentes: 1. B(E)es cerrado. 2. B(E) = (Ker B∗)0. Demostración. Sea B(E)cerrado veremos que B(E)=(Ker B∗)0. Tomamos f∈Ker B∗,g∈B(E) entonces g=B(p)para algún p∈Eyfg =0 por tanto g∈(Ker B∗)0. Para ver (Ker B∗)0⊂B(E) aplicamos Hahn-Banach, en particular se deriva la existencia de un funcional lineal l∈F∗tal que l(x) = ´ ınfy∈B(E)kx−yk. Se ve trivialmente que este funcional pertenece a Ker B∗, así que si v∈ (Ker B∗)0, necesariamente l(v) = d(v,B(E)) = 0 y así v∈B(E) = B(E). El recíproco es trivial porque (Ker B∗)0es cerrado. Como MyV∗son espacios de Banach, se tiene el resultado. Aplicando el resultado anterior tenemos una dirección del teorema, el recíproco se deriva automáticamente de la continuidad de B−1. Con todo esto podemos por fin enunciar un resultado de existencia y unicidad similar (y con motivos sobrados) al de Lax-Milgram en la sección anterior. En la demostración anterior hemos utilizado la condición LBB para demostrar la continuidad de un operador inverso, la propiedad LBB puede entonces ser observada como una especie de coercividad mixta. Teorema 3.2.6. Sea el problema (3.9) con las condiciones: 1. La forma a(•,•)es coerciva y continua. Elementos finitos en la mecánica de fluidos 3.3. Problema de Galerkin para la ecuación de Stokes 23 2. b(•,•)es continua y satisface la condición LBB. Entonces existe una única solución (u,p)para (3.9). Tomando las formas: a(u,v) = ZΩ ∇u:∇vdΩ,(3.15) b(u,q) = ZΩ q∇·udΩ.(3.16) Se puede escribir el problema de Stokes como un problema de punto silla. Así pues la existencia y unicidad de soluciones depende íntegramente de que se cumplan las condiciones de (3.2.6). La coercividad de a(•,•)se comprueba de manera directa con la norma definida para H1(Ω)en (2.4.4). La demostración1de que se cumple la condición LBB es bastante más compleja y por desgracia no habrá tiempo de realizarla. Daremos por hecho que se cumple. 3.3. Problema de Galerkin para la ecuación de Stokes Tomamos subespacios de dimensión finita, Vh×Mh. Debemos buscar en particular familias de subespacios que cumplan (3.2.6). Si no es así veremos que el método será inestable y no convergerá. Queremos ahora solucionar el problema:    Encontrar (uh,ph)∈Vh×Mhtal que: a(uh,vh)+b(vh,ph) = (f,vh)L2(Ω)para todo vh∈Vh, b(uh,qh) = 0,para todo qh∈Mh.(3.17) La coercividad de a(•,•)se cumple suponiendo que el problema variacional la cumpla en general. Si las familias de subespacios {Vh×Mh}hcumplen la condición LBB se tiene existencia de una única solución. Para ver la convergencia del método tenemos el siguiente resultado: Teorema 3.3.1. Sea el problema (3.17) con las condiciones: 1. a(•,•)es coerciva. 2. Se satisface la condición LBB. Entonces se tiene que para (u,p)solución del problema variacional de Stokes: ku−uhkV+kp−phkM≤c´ ınf vh∈Vhku−vhkV+´ ınf qh∈Mhkp−qhkM.(3.18) Se trata de una condición muy similar a la del lema de Céa. Demostración. Observemos que siendo (u,p)la solución de nuestro problema en V×My(uh,ph)la solución de la versión discreta del problema variacional, cualquier par de funciones (uI,pI)∈Vh×Mh será solución de: a(uh−uI,vh)+b(vh,ph−pI) = a(u−uI,vh)+b(vh,p−pI), b(uh−uI,qh) = b(u−uI,qh).para todo vh∈Vh,qh∈Mh(3.19) 1Una idea de la misma se puede encontrar en las páginas 159-160 de [B]. Autor: Alejandro López Nieto 24 Capítulo 3. Elementos finitos en la ecuación de Stokes Como trabajamos en un espacio de dimensión finita esto es equivalente a siendo Vh=span{Φi}con dimensión d1yMh=span{Ψj}con dimensión d2, si consideramos A= (aij) = ZΩ ∇Φj·∇ΦidΩmatriz d1×d1, B= (bij) = ZΩ Φj·∇ΨidΩmatriz d1×d2, F= (Fi) = a(u−uI,Φi)+b(Φi,p−pI)matriz d1×1, G= (Gi) = b(u−uI,Ψi)matriz d2×1. Llamaremos fyga los elementos de V∗yM∗con coordenadas en la base dual las componentes de FyGrespectivamente. Si denotamos como UyPa los vectores cuyas componentes son respectivamente las coordenadas de uh−uIyph−pIrespecto de las bases tomadas, entonces se cumple el sistema AU +BP =F, BTU=G.(3.20) Podemos escribir U=Uf+UgyP=Pf+Pg. Elegimos estos vectores de tal modo que AUf+BPf=F, BTUf=0,(3.21) AUg+BPg=0, BTUg=G.(3.22) Y llamamos uf,ug,pfypga los elementos de VhyMhtales que sus coordenadas respecto de las bases dadas son Uf,Ug,PfyPgrespectivamente. Si multiplicamos (3.21) por UT f, teniendo en cuenta que UT fBPf=PT fBTUf=0 y del hecho de que a(•,•)es coerciva nos queda αkufk2 V≤uT fAuf=UT fF≤kufkVkfkV∗. Y por tanto αkufkV≤kfkV∗.(3.23) Sean MayMblas constantes de continuidad de las formas aybrespectivamente, tenemos que ka(•,uf)kV∗≤Ma αkfkV∗(3.24) La condición LBB restringida al espacio Vh×Mhpuede escribirse como ∃β>0 independiente de h tal que ´ ınf qh∈Mh sup vh∈Vh vTBq kvhkVkqhkM≥β. Considerando la dualidad entre el espacio de matrices d1×1 y V∗ h, entonces está claro que kb(•,qh)k=sup 06=vh∈Vh b(vh,qh) kvhkV=sup 06=vh∈Vh vTBq kvhkV=kBqkV∗≥βkqhkM. Es una definición equivalente de la condición LBB. Además αkvk2 V≤a(v,v)≤Makvk2 V, por tanto α≤May así obtenemos que kpfkM≤1 βkb(•,pf)kV∗=1 βkf−a(•,uf)kV∗≤1 β1+Ma αkfkM∗≤2Ma αβ kfkM∗.(3.25) Elementos finitos en la mecánica de fluidos 3.3. Problema de Galerkin para la ecuación de Stokes 25 Ahora nos toca estudiar que sucede en la parte gdel sistema, multiplicando los términos por uT g se tiene la siguiente desigualdad αkugk2 V≤uT gAug=−uT gBpg=−pT gBTug=−pT gg≤kpgkMkgkM∗.(3.26) Por otro lado, aplicando la condición LBB como antes kpgkM≤1 βkb(•,pg)kV∗=1 βka(•,ug)kV∗≤Ma βkugkV.(3.27) Si combinamos (3.26) y (3.27), llegamos a kugkV≤Ma αβ kgkM∗.(3.28) Y una combinación de (3.27) y (3.28) nos da kpgkM≤M2 a αβ2kgkM∗.(3.29) Sumando por un lado (3.23) y (3.28) y por otro (3.25) y (3.29) y aplicando la desigualdad triangular se tiene el resultado:        kuh−uIkV≤1 αkfkV∗+Ma αβ kgkM∗, kph−pIkM≤2Ma αβ kfkV∗+M2 a αβ2kgkM∗. (3.30) Ahora bien, por la forma en que hemos construido fyg, podemos afirmar que kfkV∗≤ kfkV∗ h≤Maku−uIkV+Mbkp−pIkM, kgkM∗≤ kgkM∗ h≤Mbku−uIkV. Aplicando la desigualdad triangular queda que ku−uhkV≤C1´ ınf vh∈Vhku−vhkV+C2´ ınf qh∈Mhkp−qhkM,(3.31) kp−phkM≤C3´ ınf vh∈Vhku−vhkV+C4´ ınf qh∈Mhkp−qhkM.(3.32) Sumando las expresiones (3.31) y (3.32) y tomando cconstante suficientemente grande tenemos nuestro resultado. Con toda esta teoría queda visto que todo método numérico que cumpla las condiciones de compatibilidad convergerá, nuevamente tendremos que encontrar las proyecciones ortogonales de la solución sobre los espacios de dimensión finita. Nuestro par de espacios de elementos finitos debe aspirar a crecer hacia el espacio total en el sentido de ´ ınfvh∈Vhku−vhk→0,´ ınfqh∈Mhkp−qhk→0. Sabiendo todo esto ya podemos permitirnos plantear el problema que vamos a resolver, se tratará de un rectángulo con condiciones de contorno de tipo Dirichlet. Así la velocidad en el contorno será cero en todos los lados de dicho cuadrado salvo en uno en el que será constante en una dirección tangente a la frontera. Autor: Alejandro López Nieto 26 Capítulo 3. Elementos finitos en la ecuación de Stokes 3.4. Un problema práctico Vamos a plantear el problema que queremos resolver:    Encontrar (u,p)∈V×Mtal que: a(u,v)+b(v,p) = (f,v)L2(Ω), b(u,q) = 0,for all (v,q)∈V0×M.(3.33) Donde se tiene: Ω= (0,1)×(0,1), V=v∈H1(Ω)×H1(Ω):v(x,1) = (1,0),x∈[0,1]yv=0 en el resto de Γ, V0=v∈H1(Ω):v|Γ=0, M=q∈L2(Ω):ZΩ qdΩ=0. Nos surge la duda de si este problema tendrá realmente solución y si en tal caso será única, por lo visto anteriormente tenemos claro que el caso Dirichlet homogéneo si que tiene solución única, tomamos así el par de problemas siguiente:    Encontrar (u,p)∈V0×Mtal que: a(u,v)+b(v,p) = (f,v)L2(Ω), b(u,q) = 0,para todo (v,q)∈V0.×M(3.34)    Encontrar (u,p)∈V×Mtal que: a(u,v) = 0, b(u,q) = 0,para todo (v,q)∈V0×M.(3.35) Sean (ˆ u,p)solución de (3.34) que sabemos tiene solución única y (u0,0)solución de (3.35) que admite la solución u(x,1) = (1,0)si x∈[0,1], u(x,y) = (0,0)en el resto de Ω. Se tiene entonces que (u=ˆ u+u0,p)cumple: a(u,v)+b(v,p) = a(ˆ u,v)+b(v,p) +a(u0,v) = (f,v)L2(Ω), b(u,q) = b(ˆ u,q)+b(u0,q) = 0, u(x,1) = u0= (1,0)si x∈[0,1], u(x,y) = (0,0)en el resto de Ω. Por tanto ues solución de (3.33) y se tiene además que es la única solución. En efecto, sea ˆ u0otra solución de (3.35), entonces u0−ˆ u0es solución de:    Encontrar (u,p)∈V0×Mtal que: a(u,v) = 0, b(u,q) = 0,para todo (v,q)∈V0×M.(3.36) Sabemos que (0,p)es solución de este problema y además por Lax-Milgram es única, así que llegamos a lo que queríamos, u0=ˆ u0. Elementos finitos en la mecánica de fluidos 3.4. Un problema práctico 27 La convergencia de cualquier método será en este caso igual que en el homogéneo dado que la diferencia de soluciones del problema a resolver y del problema de Galerkin no homogéneo será una solución del problema (3.34), con todo esto podemos plantear el método de elementos finitos y ver la solución. Probaremos en primer lugar el método P1−P1. Este método resultará ser inestable puesto que, como veremos en (3.4.1), no satisface la condición LBB. Sin embargo una pequeña ampliación de los subespacios de soluciones de las velocidades nos dará un buen método convergente llamado método del MINI-elemento. Se trata de seleccionar las componentes de la velocidad en los espacios de funciones P1⊕B3donde B3denota las funciones burbuja sobre cada uno de los elementos. En nuestro elemento triangular de funciones nodales en los vértices φ1,φ2,φ3la función burbuja es el producto de estas tres funciones Φ=φ1φ2φ3, se trata de una función cúbica y considerando el elemento de referencia habitual para elementos triangulares tenemos la función burbuja de referencia ˆ Φ(ˆx,ˆy) = ˆxˆy−ˆx2ˆy−ˆxˆy2. Para una triangulación dada la dimensión del espacio B3tiene dimensión |F|=Número de caras (triángulos) de la triangulación. Los nodos para la velocidad y la presión sobre el elemento de referencia se reflejan en la figura (3.1). De esta forma tendremos 11 grados de libertad por elemento, en lugar de los 9 que tiene el método P1−P1. Figura 3.1: Nodos de presión (izquierda) y velocidades (derecha) sobre el elemento de referencia para el método MINI. Proposición 3.4.1. El método de elementos finitos P1−P1no cumple la condición LBB. Demostración. Consiste en ver que realmente existen elementos ph∈Mhtales que ZKph∇·vhdK=0 para todo vh∈Vh,(3.37) para todo K∈τhtriangulación con la que estamos trabajando. En particular, en el algoritmo que utilizaremos para la simulación más adelante trabajaremos con una triangulación del tipo de la figura (3.2). Si ph=∑|V| i=1piφiy tomamos valores de pi=1 para iun nodo cuadrado,pi=−1 para iun nodo círculo ypi=0 en los nodos triángulo, entonces se trata de un simple ejercicio de comprobación ver que, efectivamente, 0 6=ph∈Mhy cumple (3.37). Estas familias de soluciones falsas carecen de significado físico. Son ellas las que perturban la solución real del problema y la hacen comportarse de forma errática, este hecho hace que se las llame spurious modes (modos falsos). Autor: Alejandro López Nieto 34 Capítulo A. Programas utilizados rk ( 1 : 3 , : ) = d ete r ∗(C(2 ,1) ∗mx+C(2 ,2) ∗my) ; g2 ( nodos , nodos )=g2 ( nodos , nodos )+rk ( 1 : 3 , : ) ; rk ( 4 , : ) = d e t e r ∗(C(1 ,1) ∗[−1/120 ,1/120 ,0;]+C(1 , 2) ∗[−1/120 ,0 ,1/120;]) ; g3 ( k , nodos ) =g3 ( k , nodos ) +rk ( 4 , : ) ; rk ( 4 , : ) = d e t e r ∗(C(2 ,1) ∗[−1/120 ,1/120 ,0;]+C(2 , 2) ∗[−1/120 ,0 ,1/120;]) ; g4 ( k , nodos ) =g4 ( k , nodos ) +rk ( 4 , : ) ; t ( nodos )= t ( nodos ) +1.0/6∗d e t e r ; end G1=g1 . ’ ; G2=g2 . ’ ; G3=g3 . ’ ; G4=g4 . ’ ; % Construimos la mat riz por bloques A(1: nt , 1 : nt )=a ; A( nt +1:2∗nt , nt +1:2∗nt )=a ; A(2∗nt +1:2∗nt +nel ,2∗nt +1:2∗nt +nel )=m; A(2∗nt+ nel +1:2∗nt +2∗nel ,2∗nt+ nel +1:2∗nt +2∗nel )=m; A(2∗nt +1:2∗nt +nel ,2∗nt +1+2∗nel :3∗nt +2∗nel ) =g3 ; A(2∗nt+ nel +1:2∗nt +2∗nel ,2∗nt +1+2∗nel :3∗nt +2∗nel )=g4 ; A( 1 : nt ,2∗nt +1+2∗nel :3∗nt +2∗nel )=g1 ; A( nt +1:2∗nt ,2∗nt +1+2∗nel :3∗nt +2∗nel )=g2 ; A(2∗nt +1+2∗nel :3∗nt +2∗nel , 1 : nt ) =G1; A(2∗nt +1+2∗nel :3∗nt +2∗nel , nt +1:2∗nt )=G2; A(2∗nt +1+2∗nel :3∗nt +2∗nel ,2∗nt +1:2∗nt +nel )=G3; A(2∗nt +1+2∗nel :3∗nt +2∗nel ,2∗nt+nel +1:2∗nt +2∗nel )=G4; A(3∗nt +2∗nel +1 ,2∗nt +2∗nel +1:3∗nt +2∗nel )= t ; % Introdu cimos condición D i r i c h l e t en la d ir ec ci ó n x % 0 en [0 ,1] x {0} for k =1:( ns1 +1) for i =1:3∗nt +2∗nel A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 0 en {0} x [0 ,1] for k=1: ns1 +1: ns2∗ns1+ns2+1 for i =1:3∗nt +2∗nel A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 0 en {1} x [0 ,1] for k=ns1 +1: ns1 +1: nt for i =1:3∗nt +2∗nel A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 1 en [0 ,1] x {1} for k=ns2∗ns1+ns2 +1: nt for i =1:3∗nt +2∗nel A( k , i ) =0; Elementos finitos en la mecánica de fluidos 35 end b ( k ) =1; A( k , k ) =1; end % Introdu cimos condición D i r i c h l e t en la d ir ec ci ón y % 0 en [0 ,1] x {0} for k =1:( ns1 +1) for i =1:3∗nt +2∗nel A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % 0 en {0} x [0 ,1] for k=1: ns1 +1: ns2∗ns1+ns2+1 for i =1:3∗nt +2∗nel A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % 0 en {1} x [0 ,1] for k=ns1 +1: ns1 +1: nt for i =1:3∗nt +2∗nel A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % 1 en [0 ,1] x {1} for k=ns2∗ns1+ns2 +1: nt for i =1:3∗nt +2∗nel A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % Resolvemos s istem a so l =A\ b ; % Dibujamos la s olu ció n X=reshape( coord ( : , 1 ) , ns1 +1 , ns2 +1) ; vx=reshape( s o l ( 1 : n t ) , ns1 +1 , ns2 +1) ; Y=reshape( coord ( : , 2 ) , ns1 +1 , ns2 +1) ; vy=reshape( so l ( nt +1:2∗nt ) , ns1 +1 , ns2 +1) ; % Normalización de l as v elo ci dad es for i =1: ns1+1 for j =1: ns2+1 k=norm ( [ vx ( i , j ) vy ( i , j ) ] ) ; i f k<0.00000001 vy ( i , j ) =0; vx ( i , j )=vy ( i , j ) ; else vy ( i , j )=vy ( i , j ) / k ; vx ( i , j )=vx ( i , j ) / k ; Autor: Alejandro López Nieto 36 Capítulo A. Programas utilizados end end end Z=reshape( s ol (2∗nt +2∗nel +1:3∗nt +2∗nel ) , ns1 +1 , ns2 +1) ; % Gráfica de la pr esión surf (X,Y, Z) pause % Gráfica de ve lo cid ad es h=quiver (X,Y, vx , vy ) ; axis square shading interp title( ’ Solución ’ ) Programa principal para el método P1−P1. [mxx , myy , mxy]= matriz2 ; ns1 = 30; % S u b d i v i s i o n e s en l a d i r e c c i o n x ns2 = 30; % S u b d i v i s i o n e s en l a d i r e c c i o n y . d = 3; % Grados de l i b e r t a d por elemento . nt = ( ns1 +1) ∗( ns2 +1) ; % Numero t o t a l de v e r t i c e s . [ globales , x , y , nel ]= gen2 ( ns1 , ns2 ) ; mx=zeros( d , d ) ; my=zeros( d , d ) ; mx( 1 , 1 : 3 ) = −1/6; mx( 2 , : )=−mx ( 1 , : ) ; my( 1 , : ) =mx( 1 , : ) ; my( 3 , : )=−my ( 1 , : ) ; mx=−mx; my=−my; coord=zeros( nt , 2 ) ; a=zeros( nt , nt ) ; g1=zeros( nt , nt ) ; g2=zeros( nt , nt ) ; m=zeros( nel , nel ) ; G1=zeros( nt , nt ) ; G2=zeros( nt , nt ) ; A=zeros(3∗nt +1 ,3∗nt ) ; b=zeros(3∗nt +1 ,1) ; nodos=zeros(1 ,3) ; c=zeros(2 ,2) ; mk=zeros( d , d ) ; lk=zeros( d , 1 ) ; t =zeros(1 , nt ) ; so l =zeros(3∗nt , 1 ) ; coord = [ x ; y ] ’ ; % Ensamblado for k=1: nel nodos= glo ba les ( k , : ) ; [ c , C, d e t e r ]= InvDet ( coord ( nodos , : ) ) ; mk= d e t e r ∗( c ( 1 , 1) ∗mxx+c (1 , 2) ∗(mxy+mxy ’ ) +c ( 2 , 2) ∗myy) ; a ( nodos , nodos )=a ( nodos , nodos )+mk; t ( nodos )= t ( nodos ) +1.0/6∗d e t e r ; end for k=1: nel nodos= glo ba les ( k , : ) ; [ c , C, d e t e r ]= InvDet ( coord ( nodos , : ) ) ; Elementos finitos en la mecánica de fluidos 37 mk= d e t e r ∗(C( 1 , 1) ∗mx+C(1 ,2 ) ∗my) ; g1 ( nodos , nodos )=g1 ( nodos , nodos )+mk; mk= d e t e r ∗(C( 2 , 1) ∗mx+C(2 ,2 ) ∗my) ; g2 ( nodos , nodos )=g2 ( nodos , nodos )+mk; end G1=g1 . ’ ; G2=g2 . ’ ; A(1: nt , 1 : nt )=a ; A( nt +1:2∗nt , nt +1:2∗nt )=a ; A( 1 : nt ,2∗nt +1:3∗nt )=g1 ; A( nt +1:2∗nt ,2∗nt +1:3∗nt )=g2 ; A(2∗nt +1:3∗nt , 1 : nt ) =G1; A(2∗nt +1:3∗nt , nt +1:2∗nt )=G2; A(3∗nt +1 ,2∗nt +1:3∗nt ) = t ; % Introdu cimos condición D i r i c h l e t % Condiciones de contorno sobre la componente x de la vel ocida d % 0 en [0 ,1] x {0} for k =1:( ns1 +1) for i =1:3∗nt A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 0 en {1} x [0 ,1] for k=ns1 +1: ns1 +1: nt for i =1:3∗nt A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 0 en {0} x [0 ,1] for k=1: ns1 +1: ns2∗ns1+ns2+1 for i =1:3∗nt A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 0 en en [0 ,1] x {0} for k =1:( ns1 +1) for i =1:3∗nt A( k , i ) =0; end b ( k ) =0; A( k , k ) =1; end % 1 en [0 ,1] x {1} for k=ns2∗ns1+ns2 +1: nt for i =1:3∗nt A( k , i ) =0; end b ( k ) =1; Autor: Alejandro López Nieto 38 Capítulo A. Programas utilizados A( k , k ) =1; end % Condiciones de contorno sobre la componente y de la vel ocida d % 0 en [0 ,1] x {0} for k =1:( ns1 +1) for i =1:3∗nt A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % 0 en {1} x [0 ,1] for k=ns1 +1: ns1 +1: nt for i =1:3∗nt A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % 0 en en [0 ,1] x {0} for k=1: ns1 +1: ns2∗ns1+ns2+1 for i =1:3∗nt A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % 0 en [0 ,1] x {1} for k=ns2∗ns1+ns2 +1: nt for i =1:3∗nt A( nt+k , i ) =0; end b ( nt+k ) =0; A( nt+k , nt+k ) =1; end % Resolvemos s istem a so l =A\ b ; X=reshape( coord ( : , 1 ) , ns1 +1 , ns2 +1) ; vx=reshape( s o l ( 1 : n t ) , ns1 +1 , ns2 +1) ; Y=reshape( coord ( : , 2 ) , ns1 +1 , ns2 +1) ; vy=reshape( so l ( nt +1:2∗nt ) , ns1 +1 , ns2 +1) ; Z=reshape( s ol (2∗nt +1:3∗nt , : ) , ns1 +1 , ns2 +1) ; % Normalización de l as v elo ci dad es for i =1: ns1+1 for j =1: ns2+1 k=norm ( [ vx ( i , j ) vy ( i , j ) ] ) ; i f k<0.00000001 vy ( i , j ) =0; vx ( i , j )=vy ( i , j ) ; else vy ( i , j )=vy ( i , j ) / k ; vx ( i , j )=vx ( i , j ) / k ; end Elementos finitos en la mecánica de fluidos 39 end end % Gráfica de pres ión surf (X,Y, Z) shading interp pause % Gráfica de ve lo cid ad es h=quiver (X,Y, vx , vy ) ; axis ([ −0.2 1.2 −0.2 1 . 2 ] ) title( ’ Solución ’ ) Programa auxiliar gen2 que genera la malla. % Programa para generar mallas 2D. function [ globales , x , y , nel ]= gen2 ( ns1 , ns2 ) nel = ns1∗ns2 ∗2; % Numero t o t a l de elementos . nt = ( ns1 +1) ∗( ns2 +1) ; aux = zeros( nt , 1 ) ; % En el s i g u i e n t e bucle creamos la matr iz de c on ec ti vi da d . for j = 1: ns2 ind = 1+( j −1) ∗( ns1 +1) ; for i = 1:2:2∗ns1−1 elem = i + ( j −1)∗2∗ns1 ; gl ob ales ( elem , 1 : 3 ) = [ ind , ind +1 , ind+ns1 +1]; ind = ind +1; end ind = ns1 +( j −1) ∗( ns1 +1) +3; for i = 2:2:2∗ns1 elem = i + ( j −1)∗2∗ns1 ; gl ob ales ( elem , 1 : 3 ) = [ ind , ind −1,ind−ns1 −1]; ind = ind +1; end end % En el s i g u i e n t e bucle de finim os l as coordenadas de l os nodos de la malla . x1 = 0; x2 = 1; y1 = 0; y2 = 1; h1 = ( x2−x1 ) / ns1 ; h2 = ( y2−y1 ) / ns2 ; for j = 1: ns2+1 for i = 1: ns1+1 nodo = ( j −1)∗( ns1 +1) + i ; x ( nodo ) = h1 ∗( i −1) + x1 ; y ( nodo ) = h2 ∗( j −1) + y1 ; end end % Dibujamos la malla de t r i a n g u l o s . for i =1:2: nel nodos = g lo bal es ( i , : ) ; hold on f i l l ( x ( nodos ) , y ( nodos ) , ’w’ ) Autor: Alejandro López Nieto 40 Capítulo A. Programas utilizados x1 = int2str( i ) ; x2 = int2str( nodos (1 ) ) ; text( x ( nodos ( 1) )+ 0.35∗h1 , y ( nodos (1 ) ) +0.25∗h2 , x1 , ’ Color ’ ,[1 0 0 ]) ; text( x ( nodos ( 1) ) +0.1∗h1 , y ( nodos (1 ) ) +0.1∗h2 , x2 , ’ Color ’ ,[0 0 0 ]) ; end for i =2:2: nel nodos = g lo bal es ( i , : ) ; hold on f i l l ( x ( nodos ) , y ( nodos ) , ’w’ ) x1 = int2str( i ) ; text( x ( nodos ( 1) )−0.25∗h1 , y ( nodos (1 ) ) −0.25∗h2 , x1 , ’ Color ’ ,[1 0 0 ]) ; end pause hold of f return Programa auxiliar matriz2 que contiene matrices de rigidez en el elemento de referencia. function [mxx , myy , mxy] = matriz2 mxx=[0.5 −0.5 0 . ; −0.5 0.5 0 . ; 0. 0. 0 . ; ] ; myy=[0.5 0. −0.5; 0. 0. 0 . ; −0.5 0. 0 . 5 ; ] ; mxy=[0.5 0. −0.5; −0.5 0 . 0 . 5 ; 0. 0. 0 . ; ] ; return Programa auxiliar InvDet que calcula la matriz de transformación afín, su determinante, su inversa traspuesta y el producto de su inversa por su inversa traspuesta. function [ c , C, d e t e r ]= InvDet ( v ) % Calcula c=B_k ^{ −1}( B_k ^{ −1}) ^T y d et er=abs ( detB_k ) M=[v (2 , 1)−v (1 ,1) , v (3 ,1)−v (1 ,1) ; v ( 2 ,2 )−v (1 , 2) ,v ( 3 , 2)−v (1 , 2) ] ; % Matriz de la % transformación afin d e t e r = det(M) ; M=inv (M) ; C=M’ ; d e t e r =abs( d e t e r ) ; c=M∗M’ ; return Programa auxiliar pf que calcula el promedio de la función fsobre el elemento de referencia. function prom=pf ( v ) % Calcula el promedio de la fun cion f ( f =0) prom = 0; return Elementos finitos en la mecánica de fluidos Bibliografía [B] D. Braess, Finite elements: theory, fast solvers, and applications in eslasticity theory, Cambridge University Press, Cambridge, 2007. [BA] G.K. Batchelor, An introduction to fluid dynamics, Cambridge University Press, Cambridge, 1967. [BBF] D. Boffi, F. Brezzi, M. Fortin, Mixed finite element methods and applications, Springer, Berlin, 2013. [Bre] H. Brezis, Functional analysis, Sobolev spaces and partial differential equations, Springer, 2010. [ESW] H. Elman, D. Silvester, A. Wathen, Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics, Oxford University Press, Norfolk, 2006. [LB] M.G. Larson, F. Bengzon, The finite element method: theory, implementation, and applications, Springer, Berlin, 2013. [Q] A. Quarteroni, Numerical models for differential problems, Springer, Milano, 2009. 41