Eficiencia computacional de los elementos de Scott-Vogelius

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Ingeniería Aeroespacial Eficiencia computacional de los elementos de Scott-Vogelius Autor: Carlota Gutiérrez Amor Tutor: Juan Bosco García Archilla Dpto. Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Grado Ingeniería Aeroespacial Eficiencia computacional de los elementos de Scott-Vogelius Autor: Carlota Gutiérrez Amor Tutor: Juan Bosco García Archilla Catedrático de Universidad Dpto. Matemática Aplicada II Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2022 Trabajo Fin de Grado: Eficiencia computacional de los elementos de Scott-Vogelius Autor: Carlota Gutiérrez Amor Tutor: Juan Bosco García Archilla El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: Agradecimientos En primer lugar me gustaría agradecer a mi tutor, Bosco García Archilla, toda su ayuda y dedicación a lo largo de este trabajo. A mi familia, especialmente a mis padres y a mi hermana Alejandra. Carlota Gutiérrez Amor Sevilla, 2022 I Resumen Los elementos de Scott-Vogelius son unos elementos finitos diseñados para la resolución de las ecuaciones de Navier-Stokes. El interés de estos elementos es que poseen divergencia nula para las velocidades satisfaciendo automáticamente la ley de conservación de masa. En el presente trabajo se analiza la eficiencia computacional de los elementos de Scott-Vogelius frente a otros elementos de divergencia no nula acompañados de métodos de estabilización (elementos de Taylor-Hood con estabilización grad-div). III XÍndice A.3.2 Taylor-Hood con estabilización grad-div 33 A.3.3 Scott-Vogelius 34 Apéndice B Programas de MATLAB®35 B.1 Problema de Stokes 35 B.2 Problema de Oseen 39 B.3 Ecuaciones de Navier-Stokes estacionarias 44 Índice de Figuras 59 Bibliografía 61 Notación AMatriz de coeficientes lineales. bCampo vectorial de la convección en el problema de Oseen. c.c.En caso contrario. dDimensión del dominio de definición Ω. ´ ınf Ínfimo. int Interior. f Campo vectorial de las aceleraciones debidas a las fuerzas externas que actúan sobre el fluido. gen Espacio vectorial generado por. hMáximo diámetro del mallado. L2Espacio de funciones cuadrado integrables. L2 0Espacio de funciones cuadrado integrables de media nula. nNúmero de subdivisiones. nhDimensión asociada a una componente del espacio Vh. nnNúmero de nodos. nTNúmero de triángulos. NConjunto de los números naturales. mhDimensión del espacio Qh. OOrden de convergencia. pCampo escalar de la presión. PkEspacio de elementos finitos polinomiales de grado k. PlCoeficiente de la l-ésima función base del espacio Qh. QEspacio de campos escalares definidos sobre Ω. Qh Subespacio de Q de dimensión finita al que se restringe la solución numérica aproximada de la presión. RConjunto de los números reales. sup Supremo. supp Soporte. tVariable temporal. T(x,y) Producto tensorial del tensor de tercer orden T:=Ti,j,k con los vectores x,ydefinido como: ∑ j,k Ti,j,kxjyk. Tref Triángulo de referencia. TsElemento triangular. {Ts}Triangulación, conjunto de triángulos. uCampo vectorial de la velocidad. uComponente de la velocidad uen la dirección OX. Uj Coeficiente de la componente de la velocidad en la dirección OX en la j-ésima función base del espacio Vh. vComponente de la velocidad uen la dirección OY. XI XII Notación VEspacio de campos vectoriales definidos sobre Ω. Vh Subespacio de V de dimensión finita al que se restringe la solución numérica aproximada de la velocidad. Vj Coeficiente de la componente de la velocidad en la dirección OY en la j-ésima función base del espacio Vh. xVector de incógnitas. αConstante del término de reacción adicional. γParámetro de la estabilización grad-div. ∆Operador laplaciano. ηFunción test de la presión. λParámetro de escala de la presión. νViscosidad cinemática del fluido. ρDensidad, masa por unidad de volumen del fluido. φ φ φFunción test de la velocidad. φComponente en la dirección OX de φ φ φ. ψComponente en la dirección OY de φ φ φ. ΩDominio de definición de los campos. ΩClausura de Ω. ∂ΩFrontera de Ω. ∩Intersección. ∀Para todo. ∈Perteneciente a. ∇Operador gradiente. ∇·Operador divergencia. ∂ ∂tDerivada parcial respecto de t. ∂ ∂x,·xDerivada parcial respecto de x. ∂ ∂y,·yDerivada parcial respecto de y. ∅Conjunto vacío. (·,·)Producto escalar en L2, integral del producto. |·| Determinante. ||·||L2 Norma del espacio L2 o norma do, ||u||L2=ZΩ u·udΩ1/2 . :Tal que. 1 Introducción Los elementos de Scott-Vogelius son unos elementos finitos diseñados para la resolución de las ecuaciones de Navier-Stokes que rigen la dinámica de un fluido. El interés de estos elementos es que poseen divergencia nula para las velocidades satisfaciendo automáticamente la ley de conservación de masa. En este trabajo se analizan la eficiencia computacional y la convergencia del error de forma experimental de los elementos finitos de Scott-Vogelius frente a otros elementos de divergencia no nula: los elementos de Taylor-Hood con y sin estabilización grad-div. Este trabajo consta de las siguientes secciones: 1. En el capítulo 2 se presentan las ecuaciones de Navier-Stokes y los tres problemas de aplicación derivados de estas ecuaciones: el problema de Stokes, el problema de Oseen y las ecuaciones estacionarias de Navier-Stokes. 2. En el capítulo 3 se describe el método de elementos finitos, empleado para obtener una solución aproximada de un sistema de ecuaciones diferenciales. 3. En el capítulo 4 se definen los elementos finitos de Scott-Vogelius y Taylor-Hood, así como las condiciones suficientes para garantizar la existencia y unicidad de solución. 4. Los capítulos 5, 6, 7 tratan los problemas ordenados de menor a mayor complejidad (problema de Stokes, problema de Oseen y ecuaciones estacionarias de Navier-Stokes), exponiendo y analizando los resultados obtenidos. 5. Al final de este trabajo se encuentran dos apéndices. El Apéndice A recoge los resultados numéricos obtenidos de los errores cometidos en el cálculo de la velocidad, la divergencia de la velocidad y la presión y del tiempo de ejecución de los algoritmos. En el Apéndice B se proporcionan los códigos de MATLAB ® empleados. 1 2 Ecuaciones de Navier-Stokes Las ecuaciones de Navier-Stokes consisten en un sistema de ecuaciones en derivadas parciales que rigen el movimiento de los fluidos. En el caso de fluidos newtonianos (con viscosidad constante independiente de la velocidad del fluido), en flujos isotérmicos e incompresibles (a temperatura y densidad del fluido constantes en todo el espacio fluido) estas ecuaciones son ∇·u=0(2.1) ∂u ∂t+(u·∇)u=−1 ρ∇p+ν∆u+f(2.2) donde u es la velocidad y p la presión del fluido definidas en cada punto del espacio y cada instante de tiempo, ρ es la densidad, ν es la viscosidad cinemática (ambas constantes) y f son las aceleraciones debidas a las fuerzas externas que actúan sobre el fluido en cada punto e instante de tiempo. Si se considera además un flujo en régimen estacionario entonces ∂u ∂t=0 y las ecuaciones de Navier-Stokes pueden simplificarse ∇·u=0(2.3) (u·∇)u=−1 ρ∇p+ν∆u+f(2.4) Si además la velocidad del fluido es pequeña, el término no lineal de fuerzas de inercia es despreciable frente al término de las fuerzas viscosas y el sistema se reduce al llamado problema de Stokes: ∇·u=0(2.5) 0=−1 ρ∇p+ν∆u+f(2.6) Si no puede hacerse esta aproximación porque el término de fuerzas de inercia no sea despreciable, para simplificar el problema numéricamente se puede considerar el problema de Oseen: ∇·u=0(2.7) (b·∇)u=−1 ρ∇p+ν∆u+f(2.8) donde bes una función de divergencia nula ∇·b=0que es una aproximación de la solución para la u. 3 3 Método de los Elementos Finitos Las ecuaciones de Navier-Stokes constituyen un problema cuya solución está definida en una región Ω⊂Rn con int(Ω)=∅, de un espacio de funciones de dimensión infinita. Para poder obtener una solución numérica aproximada se usan conjuntamente el método de elementos finitos para discretizar el espacio de definición y el método de Galerkin para reducir la dimensión del espacio de funciones a uno de dimensión finita que permita resolver el sistema de ecuaciones diferenciales de forma aproximada. El método de elementos finitos requiere dividir la clausura del dominio de definición Ω en elementos, usualmente símplices1o prismas. El conjunto de elementos M(Ω)es tal que int(C)∩int(D) = ∅∀C,D∈M(Ω):C=D. En el presente documento, se van a considerar solo dominios bidimensionales Ω⊆R2 por lo que los recubrimientos también serán bidimensionales 2 . El mallado seleccionado es una malla uniforme triangular con diagonales en dirección suroeste a noreste (SWNE). Dado un número natural n∈N , el mallado viene dado por la triangulación {Ts}con s∈ {1,..,nT}donde nTes el número de triángulos. Figura 3.1 Malla triangular uniforme SWNE para n=4. En este tipo de mallado el número de triángulos viene dado por nT=2n2 y el número de nodos o vértices por nn= (n+1)2. El diámetro de un símplice se define como la máxima distancia entre dos puntos pertenecientes a dicho símplice. Los resultados teóricos de convergencia del error de la solución se suelen expresar respecto del 1 Sea un subconjunto abierto Ω∈Rd , un símplice en dicho espacio se define como la envolvente convexa de d+1 puntos de Ω no habiendo tres de ellos alineados. En el caso R2un símplice es un triángulo y en R3un tetraedro. [7] 2Se define el espacio de definición Ωcomo un cuadrado de lado unidad Ω= [0,1]2. 5 6Capítulo 3. Método de los Elementos Finitos máximo diámetro del mallado, denotado por h . Este parámetro es inversamente proporcional al número de subdivisiones n. 3.1 Descripción del método Para aplicar el método de elementos finitos al problema de interés es necesario realizar los siguientes dos pasos: 1. Obtener la formulación débil del problema. 2. Restringir el espacio solución a un subespacio de dimensión finita de funciones polinómicas a trozos a fin de obtener la solución aproximada. La obtención de la forma débil requiere primero establecer la integral del producto de la función test por el residuo. Para las ecuaciones de Navier-Stokes estas integrales se pueden escribir en coordenadas cartesianas como: ZΩ φuux+vuy−νuxx +px−fxdΩ(3.1) ZΩ ψuvx+vvy−νvyy +py−fydΩ(3.2) ZΩ ηux+vydΩ(3.3) donde φ y ψ son las funciones test para la velocidad y η la función test para la presión, u= (u,v) es la velocidad, p la presión y f= (fx,fy) las fuerzas externas y los subíndices denotan las derivadas parciales respecto a las coordenadas correspondientes. Después es necesario relajar la derivada, a fin de eliminar los componentes con derivadas segundas. Para ello, en el ejemplo anterior, teniendo en cuenta que uxx =∆u=∇·(∇u)y la relación ∇·(φ∇u) = ∇φ·∇u+φ∇·(∇u), se puede desarrollar la integral de φ∆ucomo ZΩ φ∆u dΩ=ZΩ ∇·(φ∇u)−∇φ·∇u dΩ por lo que la expresión (3.1) resulta ZΩ φuux+vuy+px−fxdΩ+νZΩφxux+φyuydΩ−νZΩ(∇·(φ∇u))dΩ.(3.4) Aplicando el teorema de la divergencia a la tercera integral se obtiene ZΩ φuux+vuy+px−fxdΩ+νZΩφxux+φyuydΩ−νZ∂Ω φ∇u·nds.(3.5) Teniendo en cuenta que u=0 en ∂Ω , la función test verifica que φ,ψ=0 en ∂Ω , por lo que la última integral es nula quedando ZΩ φuux+vuy+px−fxdΩ+νZΩφxux+φyuydΩ.(3.6) Además, teniendo en cuenta que φpx+ψpy=φ ψ·∇p=∇·pφ ψ−p∇·φ ψ , integrando en todo el espacio y aplicando el teorema de la divergencia ZΩ φpx+ψpydΩ=ZΩφ ψ·∇pdΩ=Z∂Ω pφ ψ·nds −ZΩ p∇·φ ψdΩ.(3.7) 3.1 Descripción del método 7 Como φ,ψ=0en ∂Ω, la primera integral del segundo término se anula ZΩ φpx+ψpydΩ=−ZΩ p∇·φ ψdΩ=ZΩ φxp+ψyp dΩ.(3.8) Además, esta igualdad debe ser cierta para φ,ψ cualesquiera, en particular para ψ≡0 , quedando la ecuación 3.6 como: ZΩ φuux+vuydΩ+νZΩφxux+φyuydΩ−ZΩ φxp dΩ−ZΩ φfxdΩ.(3.9) Operando de forma análoga para la expresión (3.2) y teniendo en cuenta que la expresión (3.3) no tiene derivadas de orden dos, el problema consiste en encontrar u∈V,p∈Q con V=H1 0(Ω)2,Q=L2 0(Ω)3 tales que verifiquen: ZΩ φuux+vuydΩ+νZΩφxux+φyuydΩ−ZΩ φxp dΩ=ZΩ φfxdΩ ZΩ ψuvx+vvydΩ+νZΩψxvx+ψyvydΩ−ZΩ ψyp dΩ=ZΩ ψfydΩ ZΩ ηux+vydΩ=0 u=0en ∂Ω (3.10) para todo φ φ φ=φ ψ∈V,η∈Qdefinidas en Ωcon φ,ψ=0en ∂Ω. A este problema se le conoce como forma débil de las ecuaciones de Navier-Stokes. Tal y como se ha comentado anteriormente al inicio de la sección, el segundo paso consiste en establecer un subespacio de funciones de dimensión finita Vh⊂V,Qh⊂Q del cual obtener la solución aproximada uh∈Vh,ph∈Qh. Estos espacios Vh,Qh tienen por elementos funciones continuas, polinómicas a trozos y de grado fijo en cada elemento. En este documento se considera un espacio de funciones cuadráticas para la velocidad (para Vh) y lineales para la presión (para Qh). Como Vh,Qh son espacios de funciones de dimensión finita, se pueden generar a partir de una base finita. Un conjunto de funciones φ φ φh=φ1 0,φ2 0,...,φnh 0,0 ψ1,0 ψ2,...,0 ψnh:gen(φ φ φh) = Vh se denomina como conjunto de funciones base de la velocidad. Análogamente, un conjunto de funciones ηh={η1,η2,...,ηmh}:gen(ηh) = Qhse denomina como conjunto de funciones base de la presión. Por ejemplo, para un espacio de funciones cuadráticas en dos dimensiones una base es el conjunto {x2,y2,xy,x,y,1} , pues cualquier función cuadrática puede escribirse como combinación lineal de las seis funciones de la base, h(x,y) = ax2+by2+cxy +dx +ey +fcon a,b,c,d,e,f∈R. 3Sea el espacio de funciones cuadrado integrables L2(Ω) = nf:Ω→R:RΩf2dΩ<∞o. Se define L2 0(Ω) = nq∈L2(Ω):(q,1) = 0o. Se define H1 0(Ω)2=nv∈{f∈L2(Ω):D f ∈L2(Ω)}2:supp(v)es compactoo como el conjunto de funciones en dos dimensiones cuyas componentes son cuadrado integrables con primera derivada cuadrado integrables que se anulan en ∂Ω. 14 Capítulo 4. Elementos de Scott-Vogelius y Taylor-Hood Se ha demostrado en [ 7 ] que en ocasiones es necesario introducir un término de reacción adicional a las ecuaciones 3.10 siendo las ecuaciones resultantes: ZΩ φuux+vuydΩ+νZΩφxux+φyuydΩ+αZΩ φu dΩ−ZΩ φxp dΩ=ZΩ φfxdΩ ZΩ ψuvx+vvydΩ+νZΩψxvx+ψyvydΩ+αZΩ ψv dΩ−ZΩ ψyp dΩ=ZΩ ψfydΩ ZΩ ηux+vydΩ=0 u=0en ∂Ω (4.6) Este término, ponderado por el parámetro α≥0 , permite garantizar la existencia y unicidad de solución en determinados problemas numéricos o mejorar la eficiencia de resolución. 4.4 Cambio de formulación de elementos de Taylor-Hood a elementos de Scott- Vogelius Para reformular el problema escrito en el espacio de elementos finitos de Taylor-Hood como uno sobre el espacio de elementos finitos de Scott-Vogelius es necesario realizar los siguientes cambios: 1. Sustituir la base de funciones del espacio de elementos de Taylor-Hood por la de Scott-Vogelius. 2. Expandir el sistema de ecuaciones obtenido a partir de las funciones base de las presiones del número de nodos (para una malla refinada en el baricentro (n+1)2+2n2 ) a tres veces el número de triángulos (para una malla refinada en el baricentro 18n2). 5 Problema de Stokes En este capítulo y los dos capítulos siguientes (6, 7) se realizan una serie de experimentos para analizar la convergencia y realizar una comparativa entre los distintos métodos empleados: elementos finitos de Taylor-Hood con y sin estabilización grad-div frente a elementos finitos de Scott-Vogelius. A partir de los códigos programados con MATLAB ® para la resolución de los problemas de Stokes y Oseen con elementos finitos de Taylor-Hood se ha programado la resolución de ambos problemas con los elementos de Scott-Vogelius y la resolución de las ecuaciones de Navier-Stokes estacionarias con ambos elementos. En esta sección se analiza el resultado obtenido del problema de Stokes para mallados con diferente número de subdivisiones ( n=6,12,24,48,96 ) para el caso de un fluido de baja viscosidad ν=10−6 . Las condiciones de contorno de Dirichlet y el término independiente fse han seleccionado para obtener como solución exacta: u=8sen2(πx)y(1−y)(1−2y) (1.2+0.8cos(4t)) −8π(y(1−y))2sen(πx)cos(πx) (1.2+0.8cos(4t)) p=λsen(πx)cos(πy) (0.5(1.2+0.8cos(4t))) donde λ es un parámetro utilizado para escalar las presiones y t=0 . El parámetro γ de la estabilización grad-div se ha fijado en 0.05 [3]. Figura 5.1 Convergencia del error en la velocidad ufrente al número de subdivisiones, n. La figura 5.1 muestra que las respectivas soluciones uh1 convergen a la solución exacta del problema para los tres espacios de elementos finitos ya que el error decrece al aumentar la dimensión de dichos espacios. Además, el orden de convergencia dado por la pendiente de la gráfica es O(1/n3) . La figura 5.2 muestra que 1Las conclusiones que se obtienen para la componente vhson equivalentes a los obtenidos para uh. 15 16 Capítulo 5. Problema de Stokes las soluciones para la presión ph convergen de forma cuadrática, O(1/n2) . Ambos resultados coinciden con el resultado teórico (véase [7]). Figura 5.2 Convergencia del error en la presión pfrente al número de subdivisiones, n. Se han realizado los experimentos para valores de λ=0 , λ=10 y λ=100 . Las gráficas 5.1 y 5.2 se consideran representativas de todos los casos de velocidades y presiones obtenidos en cuanto al orden de convergencia. Figura 5.3 Evolución del máximo error (arriba, error en velocidad u ; abajo, error en presión p ) y el tiempo de computación para distintos valores del parámetro λ de las presiones para diferentes mallados (n=6,12,24,48,96). 17 La figura 5.3 recoge los resultados obtenidos para el problema de Stokes para diferentes valores de λ (λ=0,10,100). De esta figura se extrapolan los siguientes resultados: 1. Para gradientes de presiones medios y bajos, el uso de los elementos de Taylor-Hood es la opción más eficiente desde el punto de vista computacional. El principal inconveniente de estos elementos es el elevado error del término de la divergencia de las velocidades (véase figura 5.4). Sin embargo, se observa que el empleo de un método de estabilización (en este caso estabilización grad-div) corrige este defecto, mejorando los valores obtenidos para gradientes de presión no nulos. Figura 5.4 Error de la divergencia en norma dos, ||∇·(u−uh)||L2, para distintos valores del parámetro λ. 2. Los elementos de Scott-Vogelius no aproximan los valores de las presiones tan bien como los elementos de Taylor-Hood debido a la discontinuidad de estos elementos. Sin embargo, para gradientes elevados de presiones, los elementos de Scott-Vogelius son los que mejor aproximan la solución exacta de las velocidades por lo que en estos casos resulta muy recomendable su uso. 6 Problema de Oseen En esta sección se analiza el resultado obtenido en el problema de Oseen para mallados con diferente número de subdivisiones ( n=6,12,24,48 ) para el caso ν=10−6 . Las condiciones de contorno de Dirichlet y el término independiente fse han seleccionado para obtener como solución exacta: u=8sen2(πx)y(1−y)(1−2y) (1.2+0.8cos(4t)) −8π(y(1−y))2sen(πx)cos(πx) (1.2+0.8cos(4t)) p=λsen(πx)cos(πy) (0.5(1.2+0.8cos(4t))) donde λes un parámetro utilizado para escalar las presiones y t=0. Asimismo el coeficiente bes b=8sen2(πx)y(1−y)(1−2y) (1.2+0.8cos(4t)) −8π(y(1−y))2sen(πx)cos(πx) (1.2+0.8cos(4t)). Nótese que el coeficiente bcoincide con u, lo que implica que la solución exacta elegida para el problema de Oseen es también una solución exacta para el problema de Navier-Stokes estacionario. En el problema de Oseen las soluciones también convergen a la solución exacta y lo hacen con orden O(1/n2) para las velocidades y para las presiones (véase [7]). Figura 6.1 Convergencia del error en la velocidad ufrente al número de subdivisiones, n. En las figuras 6.1 y 6.2 se refleja el orden de convergencia de los errores. Para los elementos de Taylor-Hood y de Scott-Vogelius se verifica el resultado teórico de convergencia cuadrática. El hecho de que los elementos de Taylor-Hood con estabilización grad-div converjan con un orden ligeramente inferior a dos puede deberse a que el número de subdivisiones no es suficientemente alto para observar la tendencia teórica. 19 20 Capítulo 6. Problema de Oseen Figura 6.2 Convergencia del error en la presión pfrente al número de subdivisiones, n. Se han realizado los experimentos para valores de λ=0 , λ=10 y λ=100 . Las gráficas 5.1 y 5.2 se consideran representativas de todos los casos de velocidades y presiones obtenidos en cuanto al orden de convergencia. La figura 6.3 recoge los resultados obtenidos para el problema de Stokes para diferentes valores de λ (λ=0,10,100). Figura 6.3 Evolución del máximo error (arriba, error en velocidad u ; abajo, error en presión p ) y el tiempo de computación para distintos valores del parámetro λ de las presiones para diferentes mallados (n=6,12,24,48). 21 De la figura anterior se extraen conclusiones similares a las observadas en el problema de Stokes: 1. Para gradientes de presiones nulos, el uso de los elementos de Taylor-Hood es la opción más eficiente desde el punto de vista computacional. 2. Para gradientes de presiones medios y bajos, el uso de los elementos de Taylor-Hood resulta poco eficiente debido al elevado error de la divergencia de las velocidades (véase figura 6.4). Sin embargo, se observa que el empleo de un método de estabilización (en este caso, estabilización grad-div) corrige este defecto, mejorando los valores obtenidos para gradientes de presión no nulos y siendo el método más eficiente en estos casos. Figura 6.4 Error de la divergencia en norma dos, ||∇·(u−uh)||L2, para distintos valores del parámetro λ. 3. Los elementos de Scott-Vogelius no aproximan los valores de las presiones tan bien como los elementos de Taylor-Hood debido a la discontinuidad de estos elementos. Sin embargo, para gradientes elevados de presiones, los elementos de Scott-Vogelius son los que mejor aproximan la solución exacta de las velocidades por lo que en estos casos resulta muy recomendable su uso. 7 Problema de Navier-Stokes estacionario En este capítulo se trata el problema de Navier-Stokes estacionario. Recordando el problema dado por las ecuaciones (3.11), con ψh=φh y sustituyendo uh= nh ∑ j=1Ujφj Vjφj∈ Vh×Vhyph= mh ∑ l=1 Plηl∈Qhse obtiene nh ∑ j,r=1 UrZΩ φiφjUjφrx +VjφrydΩ+ν nh ∑ j=1 UjZΩφixφjx +φiyφjydΩ− mh ∑ l=1 PlZΩ φixηldΩ=ZΩ φifxdΩ nh ∑ j,r=1 VrZΩ φiφjUjφrx +VjφrydΩ+ν nh ∑ j=1 VjZΩφixφjx +φiyφjydΩ− mh ∑ l=1 PlZΩ φiyηldΩ=ZΩ φifydΩ nh ∑ j=1 UjZΩ ηkφjxdΩ+VjZΩ ηkφjydΩ=0 nh ∑ j=1 Ujφj= nh ∑ j=1 Vjφj=0en ∂Ω(7.1) para todo i=1,...,nhy para todo k=1,...,mh. Este problema no es lineal sino cuadrático y por tanto, a diferencia de los problemas anteriores, no se puede resolver mediante eliminación gaussiana. Es necesario otro método, en este caso, el método de Newton. 7.1 Método de Newton El método de Newton es un método iterativo para el cálculo de una solución de una ecuación no lineal. Este método se basa en la obtención de un cero de una función, que representa la ecuación a resolver, y de su diferencial. La función Fconsiderada para la resolución mediante el método de Newton es F(x) = L(x,x)+Ax−f(7.2) donde A , x y f son la matriz de coeficientes lineales, el vector de incógnitas y el término independiente definidos en (3.13) (con D=0 ya que este término es sustituido por el término no lineal) y L es un tensor de tercer orden definido como 23 30 Capítulo A. Resultados numéricos Tabla A.5 Taylor-Hood, n=96 γ=0. λt max error L2error 0 7.67e+00 5.93e-06 4.88e-06 1.05e-09 1.10e-06 1.19e-06 1.67e-03 1.07e-10 10 7.81e+00 0.77785 0.77735 0.0013197 2.46e-01 2.46e-01 2.78e+02 5.70e-04 100 7.15e+00 7.7785 7.7735 0.013197 2.46e+00 2.46e+00 2.78e+03 5.70e-03 A.1.2 Taylor-Hood con estabilización grad-div Tabla A.6 Taylor-Hood, n=6γ=0.05. λt max error L2error 0 1.94e-01 0.0696 0.0628 1.58e-06 2.21e-02 2.29e-02 5.80e-05 3.90e-07 10 1.54e-01 0.4629 0.5065 0.323 1.27e-01 1.30e-01 1.79e+00 1.40e-01 100 1.33e-01 4.6062 5.0654 3.230 1.25e+00 1.28e+00 1.79e+01 1.40e+00 Tabla A.7 Taylor-Hood, n=12 γ=0.05. λt max error L2error 0 2.02e-01 0.0105 0.0099 3.36e-07 2.77e-03 2.87e-03 2.09e-05 7.36e-08 10 1.98e-01 0.0178 0.0191 0.08305 5.64e-03 5.62e-03 4.39e-01 3.60e-02 100 2.02e-01 0.1368 0.1517 0.8305 4.92e-02 4.84e-02 4.39e+00 3.60e-01 Tabla A.8 Taylor-Hood, n=24 γ=0.05. λt max error L2error 0 4.86e-01 0.00145 0.00142 6.96e-08 3.17e-04 3.31e-04 6.32e-06 1.19e-08 10 2.84e-01 0.00218 0.00233 0.0209 6.16e-04 6.22e-04 1.09e-01 9.07e-03 100 3.23e-01 0.01625 0.01658 0.2090 6.16e-03 5.28e-03 1.09e+00 9.07e-02 Tabla A.9 Taylor-Hood, n=48 γ=0.05. λt max error L2error 0 1.52e+00 1.89e-04 1.93e-04 1.86e-08 3.72e-05 3.92e-05 1.97e-06 2.04e-09 10 1.44e+00 0.00218 0.00029 0.00523 7.53e-05 7.63e-05 2.72e-02 2.27e-03 100 1.41e+00 0.01625 0.00205 0.05233 6.56e-04 6.56e-04 2.72e-01 2.27e-02 Tabla A.10 Taylor-Hood, n=96 γ=0.05. λt max error L2error 0 1.59e+01 2.44e-05 2.51e-05 4.73e-09 4.54e-06 4.80e-06 3.73e-06 4.19e-10 10 1.60e+01 3.38e-05 3.70e-05 0.00131 9.35e-06 9.48e-06 6.78e-03 5.68e-04 100 1.51e+01 0.000253 0.000256 0.01308 8.19e-05 8.19e-05 6.78e-02 5.68e-03 A.1.3 Scott-Vogelius Tabla A.11 Scott-Vogelius, n=6. λt max error L2error 0 8.18e-02 0.0697 0.0628 1.79e-05 2.21e-02 2.29e-02 3.17e-06 10 8.40e-02 0.5106 0.5607 0.4733 1.17e-01 1.23e-01 1.63e-01 100 9.20e-02 5.084 5.6072 4.7329 1.15e+00 1.21e+00 1.63e+00 A.2 Problema de Oseen 31 Tabla A.12 Scott-Vogelius, n=12. λt max error L2error 0 1.55e-01 0.0011 0.0099 6.65e-06 2.77e-03 2.87e-03 1.09e-06 10 1.55e-01 0.0125 0.0136 0.12096 3.47e-03 3.58e-03 4.14e-02 100 1.48e-01 0.0906 0.0947 1.2096 2.11e-02 2.16e-02 4.14e-01 Tabla A.13 Scott-Vogelius, n=24. λt max error L2error 0 6.14e-01 0.0015 0.0015 1.90e-06 3.17e-04 3.32e-04 3.22e-07 10 6.52e-01 0.0015 0.0015 0.0305 3.19e-04 3.33e-04 1.04e-02 100 6.39e-01 0.0019 0.0020 0.3041 4.72e-04 4.85e-04 1.04e-01 Tabla A.14 Scott-Vogelius, n=48. λt max error L2error 0 6.35e+00 0.00019 0.00019 5.22e-07 3.72e-05 3.92e-05 8.58e-08 10 6.41e+00 0.00019 0.00019 0.00761 3.72e-05 3.92e-05 2.60e-03 100 6.63e+00 0.00020 0.00020 0.07612 3.76e-05 3.96e-05 2.60e-02 Tabla A.15 Scott-Vogelius, n=96. λt max error L2error 0 1.23e+02 2.44e-05 2.51e-05 1.36e-07 4.54e-06 4.80e-06 2.19e-08 10 1.24e+02 2.44e-05 2.51e-05 0.0019 4.54e-06 4.80e-06 6.50e-04 100 1.25e+02 2.45e-05 2.52e-05 0.0190 4.54e-06 4.80e-06 6.50e-03 A.2 Problema de Oseen A.2.1 Taylor-Hood sin estabilización grad-div Tabla A.16 Taylor-Hood, n=6γ=0. λt max error L2error 0 2.47e-02 6.0192 4.7195 1.2757 9.40e-01 8.84e-01 4.06e+01 4.98e-01 10 2.62e-02 182.82 147.87 6.8071 1.26e+01 1.56e+01 1.62e+03 1.1942 100 2.45e-02 1774 1434.5 66.222 1.26e+02 1.56e+02 1.62e+04 1.09e+01 Tabla A.17 Taylor-Hood, n=12 γ=0. λt max error L2error 0 1.71e-01 0.374 0.3185 0.0497 4.31e-02 3.93e-02 4.47e+00 1.65e-02 10 1.68e-01 39.414 43.237 3.5073 3.7027 4.9738 6.17e+02 3.43e-01 100 1.69e-01 390.92 429.87 34.933 3.72e+01 4.97e+01 6.17e+03 3.4306 Tabla A.18 Taylor-Hood, n=24 γ=0. λt max error L2error 0 2.1588 0.0471 0.0453 0.0197 1.19e-02 1.19e-02 1.4031 7.14e-03 10 2.1787 6.8543 7.6112 0.3862 4.91e-01 6.58e-01 1.68e+02 3.33e-02 100 2.1676 68.276 75.704 3.7409 4.91e+00 6.58e+00 1.68e+03 3.25e-01 32 Capítulo A. Resultados numéricos Tabla A.19 Taylor-Hood, n=48 γ=0. λt max error L2error 0 5.52e+01 0.0078 0.0078 0.00803 2.68e-03 2.70e-03 4.30e-01 2.51e-03 10 5.58e+01 0.8225 1.1162 0.04143 6.52e-02 8.02e-02 4.65e+01 4.95e-03 100 5.25e+01 8.2201 11.147 0.41381 6.51e-01 8.01e-01 4.65e+02 4.27e-02 A.2.2 Taylor-Hood con estabilización grad-div Tabla A.20 Taylor-Hood, n=6γ=0.05. λt max error L2error 0 3.57e-02 0.2581 0.2786 0.245 1.12e-01 1.12e-01 2.16e-01 2.30e-02 10 3.19e-02 0.5832 0.6009 0.364 1.26e-01 1.26e-01 1.72e+00 1.70e-01 100 1.04e-01 4.6852 4.8182 2.999 5.87e-01 5.87e-01 1.70e+01 1.43e+00 Tabla A.21 Taylor-Hood, n=12 γ=0.05. λt max error L2error 0 1.85e-01 0.1653 0.159 0.19808 6.82e-02 6.81e-02 5.14e-02 6.67e-02 10 1.81e-01 0.2227 0.288 0.20772 7.00e-02 7.00e-02 4.19e-01 7.58e-02 100 1.83e-01 1.9586 2.270 0.83115 1.72e-01 1.77e-01 4.16e+00 3.67e-01 Tabla A.22 Taylor-Hood, n=24 γ=0.05. λt max error L2error 0 2.6657 0.05572 0.05502 0.0468 2.13e-02 2.17e-02 1.35e-02 1.94e-02 10 2.7428 0.06652 0.08147 0.0542 2.17e-02 2.21e-02 1.04e-01 2.14e-02 100 2.7275 0.58377 0.52259 0.2144 4.65e-02 4.74e-02 1.03e+00 9.26e-02 Tabla A.23 Taylor-Hood, n=48 γ=0.05. λt max error L2error 0 3.96e+01 0.02054 0.02067 0.03163 6.39e-03 6.45e-03 3.55e-03 2.68e-03 10 3.67e+01 0.02255 0.02125 0.03163 6.47e-03 6.52e-03 2.58e-02 6.12e-03 100 3.72e+01 0.01173 0.1248 0.05405 1.19e-02 1.16e-02 2.56e-01 2.34e-02 A.2.3 Scott-Vogelius Tabla A.24 Scott-Vogelius, n=6. λt max error L2error 0 4.55e-02 8.571 8.1928 11.2 3.7374 3.7806 4.3668 10 3.24e-02 8.5752 8.1972 11.373 3.7374 3.7806 4.3698 100 3.28e-02 8.6126 8.2263 12.922 3.7376 3.7808 4.6625 Tabla A.25 Scott-Vogelius, n=12. λt max error L2error 0 3.20e-01 1.3314 1.3512 0.76542 2.44e-01 2.53e-01 2.37e-01 10 3.61e-01 1.3325 1.3523 0.80839 2.44e-01 2.23e-01 2.41e-01 100 3.47e-01 1.3424 1.3622 1.509 2.44e-01 2.53e-01 4.77e-01 A.3 Ecuaciones de Navier-Stokes estacionarias 33 Tabla A.26 Scott-Vogelius, n=24. λt max error L2error 0 5.4353 0.6476 0.7306 0.8277 2.20e-01 2.19e-01 2.90e-01 10 5.4396 0.6476 0.7306 0.8294 2.20e-01 2.19e-01 2.90e-01 100 5.5667 0.6466 0.6306 0.8879 2.20e-01 2.19e-01 3.08e-01 Tabla A.27 Scott-Vogelius, n=48. λt max error L2error 0 4.49e+01 0.09613 0.09984 0.07368 2.15e-02 2.23e-02 1.57e-02 10 5.34e+01 0.09613 0.09984 0.07403 2.15e-02 2.23e-02 1.59e-02 100 4.28e+01 0.09613 0.09984 0.09781 2.15e-02 2.23e-02 3.04e-02 A.3 Ecuaciones de Navier-Stokes estacionarias A.3.1 Taylor-Hood sin estabilización grad-div Tabla A.28 Taylor-Hood, n=3γ=0. λt max error L2error 0 2.12e-02 0.1141 0.1753 0.0664 3.09e-02 3.98e-02 2.395 1.98e-02 10 2.91e-02 2.2493 2.9767 4.3496 5.04e-01 5.27e-01 2.44e+01 7.27e-01 100 2.93e-01 231.95 192.25 3327.4 4.67e+01 4.59e+01 1.92e+03 9.44e+02 Tabla A.29 Taylor-Hood, n=6γ=0. λt max error L2error 0 4.06e-01 0.0401 0.0483 0.0160 7.66e-03 7.94e-03 9.38e-01 3.57e-03 10 6.6522 85.796 82.083 488.4 1.45e+01 1.42e+01 1.30e+03 1.14e+02 100 6.7661 300.28 195.21 879.75 3.40e+01 3.06e+01 3.07e+03 2.22e+02 Tabla A.30 Taylor-Hood, n=12 γ=0. λt max error L2error 0 8.5636 0.0013 0.0012 0.0022 1.79e-03 1.74e-03 3.71e-01 4.31e-04 10 1.54e+02 48677 29439 2.11e+05 2.83e+03 2.77e+03 5.04e+05 2.29e+04 100 1.33e+02 329.22 282.9 6926.4 3.87e+01 3.78e+01 6.43e+03 1.12e+03 A.3.2 Taylor-Hood con estabilización grad-div Tabla A.31 Taylor-Hood, n=3γ=0.05. λt max error L2error 0 2.62e-02 0.16981 0.20904 0.16189 4.18e-02 4.53e-02 1.02e+00 4.12e-02 10 2.86e-02 0.61076 0.54925 1.3119 1.73e-01 1.53e-01 4.80e+00 4.51e-01 100 5.23e-02 5.9775 4.8919 26.094 1.59e+00 1.37e+00 4.44e+01 6.02e+00 34 Capítulo A. Resultados numéricos Tabla A.32 Taylor-Hood, n=6γ=0.05. λt max error L2error 0 3.41e-01 0.05527 0.04871 0.03046 1.17e-02 1.03e-02 2.37e-01 6.39e-03 10 4.27e-01 0.12569 0.10432 0.28895 3.09e-02 2.90e-02 1.53e+00 1.36e-01 100 3.18e-01 1.208 0.85456 2.8681 2.87e-01 2.76e-01 1.51e+01 1.32e+00 Tabla A.33 Taylor-Hood, n=12 γ=0.05. λt max error L2error 0 7.63e+00 0.01717 0.01149 0.00488 2.87e-03 2.41e-03 5.53e-02 6.58e-04 10 7.97e+00 0.03651 0.03214 0.07979 5.98e-03 5.51e-03 4.05e-01 3.59e-02 100 5.39e+00 0.26359 0.23821 0.80983 5.27e-02 4.99e-02 4.01e+00 3.59e-01 Tabla A.34 Taylor-Hood, n=24 γ=0.05. λt max error L2error 0 2.41e+02 0.00496 0.00369 6.78e-04 6.86e-04 5.88e-04 1.33e-02 7.18e-05 10 2.38e+02 0.01152 0.01032 0.02065 1.37e-03 1.20e-03 1.02e-01 9.05e-03 100 1.79e+02 0.07274 0.07129 0.20675 1.19e-02 1.05e-02 1.02e+00 9.05e-02 A.3.3 Scott-Vogelius Tabla A.35 Scott-Vogelius, n=3. λt max error L2error 0 2.86e-02 0.44653 0.39499 1.148 1.13e-01 1.07e-01 3.12e-01 10 3.26e-02 0.44761 0.39607 2.6464 1.13e-01 1.07e-01 6.92e-01 100 2.36e-02 0.45732 0.40579 18.44 1.13e-01 1.08e-01 6.18e+00 Tabla A.36 Scott-Vogelius, n=6. λt max error L2error 0 5.60e-01 0.12345 0.09873 0.27149 2.64e-02 2.27e-02 4.99e-02 10 5.45e-01 0.12347 0.09878 0.66403 2.64e-02 2.27e-02 1.71e-01 100 4.28e-01 0.12368 0.09925 4.8741 2.64e-02 2.27e-02 1.63e+00 Tabla A.37 Scott-Vogelius, n=12. λt max error L2error 0 1.65e+01 0.04277 0.03332 0.04907 5.88e-03 5.38e-03 9.15e-03 10 1.39e+01 0.04277 0.03333 0.15582 5.88e-03 5.38e-03 4.24e-02 100 1.06e+01 0.04279 0.03335 1.233 5.88e-03 5.38e-03 4.14e-01 Tabla A.38 Scott-Vogelius, n=24. λt max error L2error 0 6.54e+02 0.01195 0.00899 0.01088 1.37e-03 1.31e-03 2.01e-03 10 6.29e+02 0.01195 0.00899 0.03899 1.37e-03 1.31e-03 1.06e-02 100 4.97e+02 0.01199 0.00892 0.30862 1.38e-03 1.31e-03 1.04e-01 Apéndice B Programas de MATLAB® B.1 Problema de Stokes Código B.1 Cálculo de las matrices de los coeficientes, caso Scott-Vogelius. makronvq.m. function [A,M,Dx,Dy,Mp,Sp,J,Se,Me]=makronvq(T,z,epsilon) % INPUT % T ne x nt matriz con los vértices de cada triángulo T(i,j), % i=1,2,3, j=1,...,nt es el i-ésimo vértice del % triángulo número j. nt=2*N*N % z nn x 2 matriz con las coordenadas xy de los vértices de % la triangulación. z(n,1) y z(n,2), son las % coordenadas x e y del n-ésimo nodo. nn=(N+1)*(N+1) % epsilon parámetro de difusión % OUTPUT % A nn x nn Matriz de la forma bilineal. epsilon*(grad(phi_i),grad(phi_j)) % M nn x nn Matriz de masa % Dx,Dy n1 x nn matrices de las formas bilineales: % (\psi_i, partial_x\phi_j), (\psi_i, partial_y\phi_j) % con i=1,2,...,n1, j=1,2,...,nn, % (divergencia de la velocidad) % J nt x 1 vector con las jacobianas de los elementos. % Se ne*nt x ne matrices elementales de a(u,v)= (grad(u), grad(v)) % Me ne*nt x ne matrices elementales de masa % Mp n1xn1 matriz de masa para las presiones lineales. % Sp n1xn1 matriz de rigidez para las presiones lineales. % (grad(phi_i),grad(phi_j)) % nodos de cuadratura c1=[1 1]/3; c2=(6+sqrt(15))*[1 1]/21; c3=[(9-2*sqrt(15)),(6+sqrt(15))]/21; c4=[c3(2) c3(1)]; c5=(6-sqrt(15))*[1 1]/21; c6=[(9+2*sqrt(15)),(6-sqrt(15))]/21; c7=[ c6(2) c6(1)]; c=[c1;c2;c3;c4;c5;c6;c7]’; cx=(c(1,:)); 35 36 Capítulo B. Programas de MATLAB® cy=(c(2,:)); % pesos en los nodos de cuadratura w1=0.1125; w2=(155+sqrt(15))/2400; w3=w2; w4=w2; w5=(155-sqrt(15))/2400; w6=w5; w7=w5; w=[w1 w2 w3 w4 w5 w6 w7]’; % funciones base en los nodos de cuadratura Nc=[(1-cx-cy).*(1-2*cx-2*cy); cx.*(2*cx-1); cy.*(2*cy-1); ... 4*cx.*(1-cx-cy); 4*cx.*cy; 4*cy.*(1-cx-cy)]; Ncxi=[-(1-2*cx-2*cy)-2*(1-cx-cy); 4*cx-1; zeros(size(cx));... 4*(1-cx-cy)-4*cx; 4*cy; -4*cy]; Ncet=[-(1-2*cx-2*cy)-2*(1-cx-cy); zeros(size(cx)); 4*cy-1;... -4*cx; 4*cx; 4*(1-cx-cy)-4*cy]; ne=6; ne0=3; %velocidades cuadráticas y presiones lineales nt=size(T,2); nn=size(z,1); x1=z(:,1);y1=z(:,2);x=x1(T’);y=y1(T’); clear x1 y1;% b=reshape(b,2,1); % Derivadas parciales de x e y con respecto a xi y esta xxi=x(:,2)-x(:,1); xet=x(:,3)-x(:,1); yxi=y(:,2)-y(:,1); yet=y(:,3)-y(:,1); % Jacobiano J=(x(:,2)-x(:,1)).*(y(:,3)-y(:,1))-(x(:,3)-x(:,1)).*(y(:,2)-y(:,1)); % Derivadas parciales de xi y eta con respecto a x e y xix=(y(:,3)-y(:,1))./J; xiy=(x(:,1)-x(:,3))./J; etx=(y(:,1)-y(:,2))./J; ety=(x(:,2)-x(:,1))./J; % MATRICES ELEMENTALES % Derivadas parciales con respecto a xi y eta de las funciones base u=[-1 1 0]’; v=[-1 0 1]’; e=ones(ne,1); S1=[ 310-400 130-400 000000 -4 -4 0 8 0 0 00008-8 0 0 0 0 -8 8]/6; %(\partial_xi phi_i,\partial_xi phi_j) S2 =[ 6 1 1 -4 0 -4 10-1-440 1 -1 0 0 4 -4 -4 -4 0 8 -8 8 0 4 4 -8 8 -8 -4 0 -4 8 -8 8]/6; % S12+S12’ S3 =[ B.1 Problema de Stokes 37 30100-4 000000 10300-4 0008-80 000-880 -4 0 -4 0 0 8]/6; %(\partial_eta phi_i,\partial_eta phi_j) S12 =[ 30100-4 10-1-440 000000 -4 0 0 4 -4 4 0 0 4 -4 4 -4 0 0 -4 4 -4 4]/6; %(\partial_xi phi_i,\partial_eta phi_j) S4=[ ... 6 -1 -1 0 -4 0; ... -1 6 -1 0 0 -4; ... -1 -1 6 -4 0 0; ... 0 0 -4 32 16 16; ... -4 0 0 16 32 16; ... 0 -4 0 16 16 32]/360; % (phi_i,phi_j) Se=(kron(abs(J).*(xix.^2+xiy.^2),S1)+... kron(abs(J).*(etx.^2+ety.^2),S3)+... kron(abs(J).*(xix.*etx+xiy.*ety),S2)); Ae=epsilon*Se; % Matriz elemental de masa Me=kron(abs(J),S4);Me=reshape(Me,(ne^2)*nt,1); % índices I1=reshape(T(1:ne,:),ne*nt,1); I=I1*e’; J1=kron(T(1:ne,:)’,e); I=reshape(I,ne*ne*nt,1); J1=reshape(J1,ne*ne*nt,1); % MATRICES FEM A=sparse(I,J1,Ae,nn,nn); M=sparse(I,J1,Me,nn,nn); Me=reshape(Me,ne*nt,ne); % Divergencia ne0=3; N0c=[ones(1,size(cx,2)); 2*cx+cy-1; 2*cy+cx-1]; e0=ones(ne0,1); nn0=ne0*nt; Dxi=[-1 1 0 0 4 -4 1 1 0 -2 0 0 1 0 0 -1 1 -1]/6; Det=[-1 0 1 -4 4 0 1 0 0 -1 1 -1 1 0 1 0 0 -2]/6; II01= (1:(ne0*nt))’; %reshape(T(1:ne0,:),nt*ne0,1); I01=reshape(kron(II01,e’),nt*ne0*ne,1); J01=reshape(kron(T(1:ne,:)’,e0),ne0*ne*nt,1); Dxe=kron(abs(J).*xix,Dxi)+kron(abs(J).*etx,Det); Dye=kron(abs(J).*xiy,Dxi)+kron(abs(J).*ety,Det); Dx=sparse(I01,J01,reshape(Dxe,numel(Dxe),1),nn0,nn); Dy=sparse(I01,J01,reshape(Dye,numel(Dye),1),nn0,nn); 38 Capítulo B. Programas de MATLAB® % PRESIÓN % Matriz de masa de la presión S04=(1/24)*[12 0 0;0 2 1; 0 1 2]; M0e=kron(abs(J),S04); I10=(1:(ne0*nt))’; J0=kron(reshape(I10,ne0,nt)’,e0); IIq=reshape(I10*e0’,numel(I10)*numel(e0),1); JJq=reshape(J0,numel(J0),1); Mp=sparse(IIq,JJq,reshape(M0e,numel(M0e),1),nn0,nn0); % Matriz de rigidez para la presión ex0=[0 2 1]’; ey0=[0 1 2]’; S01=0.5*ex0*ex0’; S03=0.5*ey0*ey0’; S02=0.5*ex0*ey0’ + 0.5*ey0*ex0’; S0e=(kron(abs(J).*(xix.^2+xiy.^2),S01)+... kron(abs(J).*(etx.^2+ety.^2),S03)+... kron(abs(J).*(xix.*etx+xiy.*ety),S02)); Sp=sparse(IIq,JJq,reshape(S0e,numel(S0e),1),nn0,nn0); Código B.2 Resolución Problema de Stokes, caso Scott-Vogelius. pruesto2d.m. N=96; % número de subdivisiones n uex=@bwinb; % solución exacta force=@forceb; % término independiente epsilon=1e-6; % parámetro de difusión alpha=0; % constante del término de reacción adicional lambda=100; % constante para escalar las presiones tf=0; global LANU global ELALPHA global ELLAMBDA LANU=epsilon; ELALPHA=alpha; ELLAMBDA=lambda; % MALLADO [T1,z1,is1,ldir]=swne(N); % mallado regular con diagonales swne [T1b,z1b]=red_sv(T1,z1); % mallado refinado en el baricentro [T,z,is]=datosq(T1b,z1b’,is1’,is1’); z=z’; is=is’; x=z(:,1); y=z(:,2); X=x(is(1:ldir,:)); Y=y(is(1:ldir,:)); uis=zeros(size(is(1:ldir,:))); vis=uis; [u00,v00]=uex(x,y,0); u0=u00; v0=v00; ne0=3; nt=size(T,2); % Índices necesarios para representación presiones en Scott-Vogelius Ip=reshape(T(1:ne0,1:nt),nt*ne0,1); z1p=z1b(Ip,:); B.2 Problema de Oseen 39 % Resolución tic; [u,v,p,rho,Mn,Mp,Se,Me]=sto2d_steady(T,z,is,ldir,epsilon,force,uis,vis,alpha); t_gal=toc TT=[T([1 4 6],:),T([2 5 4],:), T([3 6 5],:), T([4 5 6],:)]; % REPRESENTACIÓN Y ERRORES figure(1); clf; trimesh(TT’,z(:,1),z(:,2),u); % solución u figure(2); clf; trimesh(TT’,z(:,1),z(:,2),v); % solución v [ue,ve]=uex(x,y,tf); [~,~,pe]=uex(z1p(:,1),z1p(:,2),tf); ex=abs(u-ue); ey=abs(v-ve); figure(3); clf; trimesh(TT’,z(:,1),z(:,2),log10(ex)); % error u figure(4); clf; trimesh(TT’,z(:,1),z(:,2),log10(ey)); % error v nt=length(p)/3; Tp=(1:3*nt)’; Tp=reshape(Tp,3,nt); Iph=kron(ones(7,1),3*(1:nt)’)-2+[zeros(3*nt,1);ones(2*nt,1);2*ones(2*nt,1)]; Jph=kron(ones(7,1),3*(1:nt)’)-2+[zeros(nt,1);ones(nt,1);2*ones(nt,1);zeros(nt ,1);ones(nt,1);zeros(nt,1);2*ones(nt,1)]; ph=sparse(Iph,Jph, [ones(nt,1);-ones(2*nt,1);ones(4*nt,1)])*p; figure(7); clf; h=trimesh(Tp’,z1p(:,1),z1p(:,2),ph, ’FaceColor’, [0.85 0.85 0.85]); % solución p % Máximo error format short g; maxe=[max(abs(ex)),max(abs(ey)),max(abs(pe-ph))], format short e I1=reshape(T(1:6,:),6*size(T,2),1); I=kron(ones(6,1),I1); J1=kron(T(1:6,:)’,ones(6,1)); J=reshape(J1,numel(J1),1); S=sparse(I,J,Se,length(z),length(z)); ex=ue-u; ey=ve-v; ep=pe-ph; uv=[u;v]; % Error en u, v y p ERRSL2=[sqrt(ex’*(Mn*ex)), sqrt(ey’*(Mn*ey)), sqrt(ep’*(Mp*ep))] B.2 Problema de Oseen Código B.3 Cálculo de las matrices de los coeficientes, caso Scott-Vogelius. makronvq.m. function [A,M,Dx,Dy,Mp,Sp,J,Se,Me,Xi,Eta,bC,D,De]=makronvq(T,z,epsilon,bvar,t) % INPUT % T ne x nt matriz con los vértices de cada triángulo T(i,j), % i=1,2,3, j=1,...,nt es el i-ésimo vértice del % triángulo número j. nt=2*N*N % z nn x 2 matriz con las coordenadas xy de los vértices de % la triangulación. z(n,1) y z(n,2), son las 46 Capítulo B. Programas de MATLAB® -4 0 -4 0 0 8]/6; %(\partial_eta phi_i,\partial_eta phi_j). S12 =[ 30100-4 10-1-440 000000 -4 0 0 4 -4 4 0 0 4 -4 4 -4 0 0 -4 4 -4 4]/6; %(\partial_xi phi_i,\partial_eta phi_j). S4=[ ... 6 -1 -1 0 -4 0; ... -1 6 -1 0 0 -4; ... -1 -1 6 -4 0 0; ... 0 0 -4 32 16 16; ... -4 0 0 16 32 16; ... 0 -4 0 16 16 32]/360; % (phi_i,phi_j). Se=(kron(abs(J).*(xix.^2+xiy.^2),S1)+... kron(abs(J).*(etx.^2+ety.^2),S3)+... kron(abs(J).*(xix.*etx+xiy.*ety),S2)); Ae=epsilon*Se; Me=kron(abs(J),S4);Me=reshape(Me,(ne^2)*nt,1); % Índices I1=reshape(T(1:ne,:),ne*nt,1); I=I1*e’; J1=kron(T(1:ne,:)’,e); I=reshape(I,ne*ne*nt,1); J1=reshape(J1,ne*ne*nt,1); % MATRICES FEM A=sparse(I,J1,Ae,nn,nn); M=sparse(I,J1,Me,nn,nn); Me=reshape(Me,ne*nt,ne); % GRAD-DIV GD11e=(kron(abs(J).*(xix.^2),S1)+... kron(abs(J).*(etx.^2),S3)+... kron(abs(J).*(xix.*etx),S2)); GD11=sparse(I,J1,reshape(GD11e,numel(GD11e),1),nn,nn); GD22e=(kron(abs(J).*(xiy.^2),S1)+... kron(abs(J).*(ety.^2),S3)+... kron(abs(J).*(xiy.*ety),S2)); GD22=sparse(I,J1,reshape(GD22e,numel(GD11e),1),nn,nn); GD12e=(kron(abs(J).*(xix.*xiy),S1)+... kron(abs(J).*(etx.*ety),S3)+... kron(abs(J).*(xix.*ety),S12)+... kron(abs(J).*(xiy.*etx),S12’)); GD12=sparse(I,J1,reshape(GD12e,numel(GD12e),1),nn,nn); % DIVERGENCIA ne0=3; N0c=[1-cx-cy; cx; cy]; e0=ones(ne0,1); nn0=max(max(T(1:ne0,:))); Dxi=[-1 0 0 1 1 -1 0 1 0 -1 1 -1 0 0 0 0 2 -2]/6; B.3 Ecuaciones de Navier-Stokes estacionarias 47 Det=[-1 0 0 -1 1 1 0 0 0 -2 2 0 0 0 1 -1 1 -1]/6; II01=reshape(T(1:ne0,:),nt*ne0,1); I01=reshape(kron(II01,e’),nt*ne0*ne,1); J01=reshape(kron(T(1:ne,:)’,e0),ne0*ne*nt,1); Dxe=kron(abs(J).*xix,Dxi)+kron(abs(J).*etx,Det); Dye=kron(abs(J).*xiy,Dxi)+kron(abs(J).*ety,Det); Dx=sparse(I01,J01,reshape(Dxe,numel(Dxe),1),nn0,nn); Dy=sparse(I01,J01,reshape(Dye,numel(Dye),1),nn0,nn); % PRESIÓN S04=(1/24)*[2 1 1;1 2 1; 1 1 2]; M0e=kron(abs(J),S04); I10=reshape(T(1:ne0,1:nt),nt*ne0,1); J0=kron(T(1:ne0,1:nt)’,e0); IIq=reshape(I10*e0’,numel(I10)*numel(e0),1); JJq=reshape(J0,numel(J0),1); Mp=sparse(IIq,JJq,reshape(M0e,numel(M0e),1),nn0,nn0); ex0=[-1 1 0]’; ey0=[-1 0 1]’; S01=0.5*ex0*ex0’; S03=0.5*ey0*ey0’; S02=0.5*ex0*ey0’ + 0.5*ey0*ex0’; S0e=(kron(abs(J).*(xix.^2+xiy.^2),S01)+... kron(abs(J).*(etx.^2+ety.^2),S03)+... kron(abs(J).*(xix.*etx+xiy.*ety),S02)); Sp=sparse(IIq,JJq,reshape(S0e,numel(S0e),1),nn0,nn0); % TÉRMINO NO LINEAL % (phi_i,phi_j phi_kx) (phi_i, phi_j phi_ky) matxi=[-13/420 -1/140 0 4/105 1/105 -1/105; 1/280 -1/280 0 0 1/630 -1/630; 1/280 11/2520 0 -1/126 -2/315 2/315; -2/105 -4/315 0 2/63 -2/315 2/315; -1/210 -1/126 0 4/315 -2/105 2/105; -2/105 -2/315 0 8/315 -4/315 4/315; 1/280 -1/280 0 0 1/630 -1/630; 1/140 13/420 0 -4/105 1/105 -1/105; -11/2520 -1/280 0 1/126 -2/315 2/315; 4/315 2/105 0 -2/63 -2/315 2/315; 2/315 2/105 0 -8/315 -4/315 4/315; 1/126 1/210 0 -4/315 -2/105 2/105; 1/280 11/2520 0 -1/126 -2/315 2/315; -11/2520 -1/280 0 1/126 -2/315 2/315; 1/140 -1/140 0 0 1/21 -1/21; 1/126 -1/126 0 0 -2/315 2/315; 2/315 -4/315 0 2/315 2/105 -2/105; 4/315 -2/315 0 -2/315 2/105 -2/105; -2/105 -4/315 0 2/63 -2/315 2/315; 4/315 2/105 0 -2/63 -2/315 2/315; 1/126 -1/126 0 0 -2/315 2/315; -4/63 4/63 0 0 16/315 -16/315; -2/315 2/63 0 -8/315 16/315 -16/315; -2/63 2/315 0 8/315 16/315 -16/315; -1/210 -1/126 0 4/315 -2/105 2/105; 2/315 2/105 0 -8/315 -4/315 4/315; 48 Capítulo B. Programas de MATLAB® 2/315 -4/315 0 2/315 2/105 -2/105; -2/315 2/63 0 -8/315 16/315 -16/315; 4/105 4/63 0 -32/315 16/105 -16/105; -2/315 2/315 0 0 8/105 -8/105; -2/105 -2/315 0 8/315 -4/315 4/315; 1/126 1/210 0 -4/315 -2/105 2/105; 4/315 -2/315 0 -2/315 2/105 -2/105; -2/63 2/315 0 8/315 16/315 -16/315; -2/315 2/315 0 0 8/105 -8/105; -4/63 -4/105 0 32/315 16/105 -16/105]; mateta=[-13/420 0 -1/140 -1/105 1/105 4/105; 1/280 0 11/2520 2/315 -2/315 -1/126; 1/280 0 -1/280 -1/630 1/630 0; -2/105 0 -2/315 4/315 -4/315 8/315; -1/210 0 -1/126 2/105 -2/105 4/315; -2/105 0 -4/315 2/315 -2/315 2/63; 1/280 0 11/2520 2/315 -2/315 -1/126; 1/140 0 -1/140 -1/21 1/21 0; -11/2520 0 -1/280 2/315 -2/315 1/126; 4/315 0 -2/315 -2/105 2/105 -2/315; 2/315 0 -4/315 -2/105 2/105 2/315; 1/126 0 -1/126 2/315 -2/315 0; 1/280 0 -1/280 -1/630 1/630 0; -11/2520 0 -1/280 2/315 -2/315 1/126; 1/140 0 13/420 -1/105 1/105 -4/105; 1/126 0 1/210 2/105 -2/105 -4/315; 2/315 0 2/105 4/315 -4/315 -8/315; 4/315 0 2/105 2/315 -2/315 -2/63; -2/105 0 -2/315 4/315 -4/315 8/315; 4/315 0 -2/315 -2/105 2/105 -2/315; 1/126 0 1/210 2/105 -2/105 -4/315; -4/63 0 -4/105 -16/105 16/105 32/315; -2/315 0 2/315 -8/105 8/105 0; -2/63 0 2/315 -16/315 16/315 8/315; -1/210 0 -1/126 2/105 -2/105 4/315; 2/315 0 -4/315 -2/105 2/105 2/315; 2/315 0 2/105 4/315 -4/315 -8/315; -2/315 0 2/315 -8/105 8/105 0; 4/105 0 4/63 -16/105 16/105 -32/315; -2/315 0 2/63 -16/315 16/315 -8/315; -2/105 0 -4/315 2/315 -2/315 2/63; 1/126 0 -1/126 2/315 -2/315 0; 4/315 0 2/105 2/315 -2/315 -2/63; -2/63 0 2/315 -16/315 16/315 8/315; -2/315 0 2/63 -16/315 16/315 -8/315; -4/63 0 4/63 -16/315 16/315 0]; matxi=reshape(matxi’,size(matxi))’; mateta=reshape(mateta’,size(mateta))’; NLxe=kron(abs(J).*xix,matxi)+kron(abs(J).*etx,mateta); NLye=kron(abs(J).*xiy,matxi)+kron(abs(J).*ety,mateta); I0NL=reshape(T,ne*nt,1)*ones(1,ne^2); J0NL=kron(e’,kron(T’,e))+nn*kron(T’-1,ones(ne)); NLx=sparse(I0NL,J0NL,NLxe); NLy=sparse(I0NL,J0NL,NLye); B.3 Ecuaciones de Navier-Stokes estacionarias 49 Código B.6 Cálculo de la solución de las ecuaciones de Navier-Stokes estacionarias. ns2d_steady.m. function [u,v,p,rho,Mn,Mp,Se,Me,GDn]=ns2d_steady(T,z,is,ldir,epsilon,... bvar,force,uis,vis,alpha,gamma) nn=length(z); t=0; [An,Mn,Dxn,Dyn,Mp,Sp,GD11n,GD12n,GD22n,J,Se,Me,NLx,NLy]=makronvq(T,z,epsilon); % matrices [fn,fxc,fyc]=fkronq(T,z,force,J,0); ne=6; % número de nodos por elemento para velocidad cuadrática nbe=3; % número de nodos por elemento para bc Ilong=ones(length(z),1); u=zeros(length(z),1); v=u; for j=1:nbe Ilong(is(1:ldir,j))=0; u(is(1:ldir,j))=uis(:,j); v(is(1:ldir,j))=vis(:,j); end I=find(Ilong); I0=find(~Ilong); A1n=An+alpha*Mn; ff=fn; A1=A1n(I,I); Dx=Dxn(:,I); Dy=Dyn(:,I); Mpe=Mp*ones(size(Mp,2),1); % this is to obatain zero-mean presures GD=[GD11n(I,I), GD12n(I,I); GD12n(I,I)’, GD22n(I,I)]; GDn=[GD11n, GD12n; GD12n, GD22n]; A=[A1,sparse(1,1,0,size(A1,1),size(A1,2)),-Dx’,sparse(1,1,0,size(A1,1),1);... sparse(1,1,0,size(A1,1),size(A1,2)),A1,-Dy’,sparse(1,1,0,size(A1,1),1);... -Dx,-Dy,sparse(1,1,0,size(Dx,1),size(Dx,1)),Mpe; ... sparse(1,1,0,1,2*size(A1,2)),Mpe’,0]; A(1:2*length(I),1:2*length(I))=A(1:2*length(I),1:2*length(I))+gamma*GD; % Matriz NL jcut=kron(ones(length(I),1),I); mcut=kron(I,ones(length(I),1)); jcut=jcut+nn*(mcut-1); NLxc=NLx(I,jcut); NLyc=NLy(I,jcut); NL=sparse(1,1,0,size(A,1),4*length(I)^2); Zs=sparse(1,1,0,size(NLxc,1),size(NLxc,2)); NL(1:2*length(I),:)=[NLxc Zs NLyc Zs; Zs NLxc Zs NLyc]; K=I+nn; if norm(uis)~=0 || norm(vis)~=0 ff(I)=ff(I)-A1n(I,I0)*u(I0) - gamma*(GD11n(I,I0)*u(I0) + GD12n(I,I0)*v(I0)) ; ff(K)=ff(K)-A1n(I,I0)*v(I0) - gamma*(GD12n(I0,I)’*u(I0) + GD22n(I,I0)*v(I0) ); ffp=(Dxn(:,I0)*u(I0) + Dyn(:,I0)*v(I0)); else 50 Capítulo B. Programas de MATLAB® ffp=zeros(size(Dxn,1),1); end b=[ff([I;K]);ffp;0]; % Resolución método de Newton fval=@(x)value(A,NL,x,length(I),b); dfval=@(x)dvalue(A,NLxc,NLyc,x,length(I)); aux=zeros(size(A,1),1); [x,iflag,iter,normf]=newton(fval,dfval,aux); u(I)=x(1:length(I)); v(I)=x(length(I)+1:2*length(I)); p=x(2*length(I)+1:length(x)-1); rho=x(length(x)); Código B.7 Resolución ecuaciones estacionarias de Navier-Stokes, Taylor-Hood. pruste2d.m. clear, close all, clc N=24; % número de subdivisiones bvar=@bwinb; % término de convección uex=@bwinb; % solución exacta force=@forceb; % término independiente epsilon=1e-6; % parámetro de difusión alpha=10; % constante del término de reacción adicional lambda=100; % parámetro para escalar las presiones gamma=0.05; % parámetro del método de estabilización grad-div tf=0; global LANU global ELALPHA global ELLAMBDA LANU=epsilon; ELALPHA=alpha; ELLAMBDA=lambda; % Mallado [T1,z1,is1,ldir]=swne(N); % mallado regular con diagonales swne [T1b,z1b]=red_sv(T1,z1); % mallado refinado en el baricentro [T,z,is]=datosq(T1b,z1b’,is1’,is1’); z=z’; is=is’; x=z(:,1); y=z(:,2); X=x(is(1:ldir,:)); Y=y(is(1:ldir,:)); uis=zeros(size(is(1:ldir,:))); vis=uis; [u00,v00]=uex(x,y,0); u0=u00; v0=v00; tic; [u,v,p,rho,Mn,Mp,Se,Me,GDn]ns2d_steady(T,z,is,ldir,epsilon,bvar,force,uis,vis, alpha,gamma); t_gal=toc % GRÁFICAS B.3 Ecuaciones de Navier-Stokes estacionarias 51 figure(1); clf; trimesh(T1’,z1(:,1),z1(:,2),u(1:length(z1))); TT=[T([1 4 6],:),T([2 5 4],:), T([3 6 5],:), T([4 5 6],:)]; figure(2); clf; trimesh(T1’,z1(:,1),z1(:,2),v(1:length(z1))); [ue,ve]=uex(x,y,tf); [ue1,ve1,pe]=uex(z1b(:,1),z1b(:,2),tf); ex=abs(u-ue); ey=abs(v-ve); figure(3); clf; trimesh(T’,z(:,1),z(:,2),ue); figure(4); clf; trimesh(T’,z(:,1),z(:,2),ve); figure(6); clf; trimesh(T1’,z1(:,1),z1(:,2),log10(abs(p-pe))) format short g; maxe=[max(abs(ex)),max(abs(ey)),max(abs(pe-p))], format short e I1=reshape(T(1:6,:),6*size(T,2),1); I=kron(ones(6,1),I1); J1=kron(T(1:6,:)’,ones(6,1)); J=reshape(J1,numel(J1),1); S=sparse(I,J,Se,length(z),length(z)); ex=ue-u; ey=ve-v; ep=pe-p; uv=[u;v]; % Error en u, v, divergencia, presión ERRSL2=[sqrt(ex’*(Mn*ex)), sqrt(ey’*(Mn*ey)), sqrt(uv’*(GDn*uv)), sqrt(ep’*(Mp* ep))] Código B.8 Cálculo de las matrices, caso Scott-Vogelius. makronvq.m. function [A,M,Dx,Dy,Mp,Sp,GD11,GD12,GD22,J,Se,Me,NLx,NLy]=makronvq(T,z,epsilon) % INPUT y OUTPUT ídem que para el caso de Taylor-Hood %nodos de cuadratura c1=[1 1]/3; c2=(6+sqrt(15))*[1 1]/21; c3=[(9-2*sqrt(15)),(6+sqrt(15))]/21; c4=[c3(2) c3(1)]; c5=(6-sqrt(15))*[1 1]/21; c6=[(9+2*sqrt(15)),(6-sqrt(15))]/21; c7=[ c6(2) c6(1)]; c=[c1;c2;c3;c4;c5;c6;c7]’; cx=c(1,:); cy=(c(2,:)); % pesos en los nodos de cuadratura w1=0.1125; w2=(155+sqrt(15))/2400; w3=w2; w4=w2; w5=(155-sqrt(15))/2400; w6=w5; w7=w5; w=[w1 w2 w3 w4 w5 w6 w7]’; % funciones base en los nodos de cuadratura Nc=[(1-cx-cy).*(1-2*cx-2*cy); cx.*(2*cx-1); cy.*(2*cy-1); ... 4*cx.*(1-cx-cy); 4*cx.*cy; 4*cy.*(1-cx-cy)]; Ncxi=[-(1-2*cx-2*cy)-2*(1-cx-cy); 4*cx-1; zeros(size(cx));... 4*(1-cx-cy)-4*cx; 4*cy; -4*cy]; Ncet=[-(1-2*cx-2*cy)-2*(1-cx-cy); zeros(size(cx)); 4*cy-1;... -4*cx; 4*cx; 4*(1-cx-cy)-4*cy]; ne=6; ne0=3; 52 Capítulo B. Programas de MATLAB® nt=size(T,2); nn=size(z,1); x1=z(:,1);y1=z(:,2);x=x1(T’);y=y1(T’); clear x1 y1;% b=reshape(b,2,1); % Parciales de x e y con respecto a xi y eta xxi=x(:,2)-x(:,1); xet=x(:,3)-x(:,1); yxi=y(:,2)-y(:,1); yet=y(:,3)-y(:,1); % JACOBIANO J=(x(:,2)-x(:,1)).*(y(:,3)-y(:,1))-(x(:,3)-x(:,1)).*(y(:,2)-y(:,1)); % Parciales de xi y eta con respecto a x e y xix=(y(:,3)-y(:,1))./J; xiy=(x(:,1)-x(:,3))./J; etx=(y(:,1)-y(:,2))./J; ety=(x(:,2)-x(:,1))./J; % MATRICES ELEMENTALES u=[-1 1 0]’; v=[-1 0 1]’; e=ones(ne,1); S1=[ 310-400 130-400 000000 -4 -4 0 8 0 0 00008-8 0 0 0 0 -8 8]/6; %(\partial_xi phi_i,\partial_xi phi_j). S2 =[ 6 1 1 -4 0 -4 10-1-440 1 -1 0 0 4 -4 -4 -4 0 8 -8 8 0 4 4 -8 8 -8 -4 0 -4 8 -8 8]/6; % S12+S12’ S3 =[ 30100-4 000000 10300-4 0008-80 000-880 -4 0 -4 0 0 8]/6; %(\partial_eta phi_i,\partial_eta phi_j). S12 =[ 30100-4 10-1-440 000000 -4 0 0 4 -4 4 0 0 4 -4 4 -4 0 0 -4 4 -4 4]/6; %(\partial_xi phi_i,\partial_eta phi_j). S4=[ ... 6 -1 -1 0 -4 0; ... -1 6 -1 0 0 -4; ... -1 -1 6 -4 0 0; ... 0 0 -4 32 16 16; ... -4 0 0 16 32 16; ... 0 -4 0 16 16 32]/360; % (phi_i,phi_j). % DIFUSIÓN (grad(phi_i),grad(phi_j)). B.3 Ecuaciones de Navier-Stokes estacionarias 53 Se=(kron(abs(J).*(xix.^2+xiy.^2),S1)+... kron(abs(J).*(etx.^2+ety.^2),S3)+... kron(abs(J).*(xix.*etx+xiy.*ety),S2)); Ae=epsilon*Se; % MASA Me=kron(abs(J),S4);Me=reshape(Me,(ne^2)*nt,1); % Índices I1=reshape(T(1:ne,:),ne*nt,1); I=I1*e’; J1=kron(T(1:ne,:)’,e); I=reshape(I,ne*ne*nt,1); J1=reshape(J1,ne*ne*nt,1); % MATRICES FEM A=sparse(I,J1,Ae,nn,nn); M=sparse(I,J1,Me,nn,nn); Me=reshape(Me,ne*nt,ne); % GRAD-DIV GD11e=(kron(abs(J).*(xix.^2),S1)+... kron(abs(J).*(etx.^2),S3)+... kron(abs(J).*(xix.*etx),S2)); GD11=sparse(I,J1,reshape(GD11e,numel(GD11e),1),nn,nn); GD22e=(kron(abs(J).*(xiy.^2),S1)+... kron(abs(J).*(ety.^2),S3)+... kron(abs(J).*(xiy.*ety),S2)); GD22=sparse(I,J1,reshape(GD22e,numel(GD11e),1),nn,nn); GD12e=(kron(abs(J).*(xix.*xiy),S1)+... kron(abs(J).*(etx.*ety),S3)+... kron(abs(J).*(xix.*ety),S12)+... kron(abs(J).*(xiy.*etx),S12’)); GD12=sparse(I,J1,reshape(GD12e,numel(GD12e),1),nn,nn); % DIVERGENCIA ne0=3; N0c=[ones(1,size(cx,2)); 2*cx+cy-1; 2*cy+cx-1]; e0=ones(ne0,1); nn0=ne0*nt; Dxi=[-1 1 0 0 4 -4 1 1 0 -2 0 0 1 0 0 -1 1 -1]/6; Det=[-1 0 1 -4 4 0 1 0 0 -1 1 -1 1 0 1 0 0 -2]/6; II01=(1:(ne0*nt))’; I01=reshape(kron(II01,e’),nt*ne0*ne,1); J01=reshape(kron(T(1:ne,:)’,e0),ne0*ne*nt,1); Dxe=kron(abs(J).*xix,Dxi)+kron(abs(J).*etx,Det); Dye=kron(abs(J).*xiy,Dxi)+kron(abs(J).*ety,Det); Dx=sparse(I01,J01,reshape(Dxe,numel(Dxe),1),nn0,nn); Dy=sparse(I01,J01,reshape(Dye,numel(Dye),1),nn0,nn); 54 Capítulo B. Programas de MATLAB® % PRESIÓN S04=(1/24)*[12 0 0;0 2 1; 0 1 2]; M0e=kron(abs(J),S04); I10=(1:(ne0*nt))’; J0=kron(reshape(I10,ne0,nt)’,e0); IIq=reshape(I10*e0’,numel(I10)*numel(e0),1); JJq=reshape(J0,numel(J0),1); Mp=sparse(IIq,JJq,reshape(M0e,numel(M0e),1),nn0,nn0); ex0=[0 2 1]’; ey0=[0 1 2]’; S01=0.5*ex0*ex0’; S03=0.5*ey0*ey0’; S02=0.5*ex0*ey0’ + 0.5*ey0*ex0’; S0e=(kron(abs(J).*(xix.^2+xiy.^2),S01)+... kron(abs(J).*(etx.^2+ety.^2),S03)+... kron(abs(J).*(xix.*etx+xiy.*ety),S02)); Sp=sparse(IIq,JJq,reshape(S0e,numel(S0e),1),nn0,nn0); % TÉRMINO NO LINEAL % (phi_i,phi_j phi_kx) (phi_i, phi_j phi_ky) matxi=[-13/420 -1/140 0 4/105 1/105 -1/105; 1/280 -1/280 0 0 1/630 -1/630; 1/280 11/2520 0 -1/126 -2/315 2/315; -2/105 -4/315 0 2/63 -2/315 2/315; -1/210 -1/126 0 4/315 -2/105 2/105; -2/105 -2/315 0 8/315 -4/315 4/315; 1/280 -1/280 0 0 1/630 -1/630; 1/140 13/420 0 -4/105 1/105 -1/105; -11/2520 -1/280 0 1/126 -2/315 2/315; 4/315 2/105 0 -2/63 -2/315 2/315; 2/315 2/105 0 -8/315 -4/315 4/315; 1/126 1/210 0 -4/315 -2/105 2/105; 1/280 11/2520 0 -1/126 -2/315 2/315; -11/2520 -1/280 0 1/126 -2/315 2/315; 1/140 -1/140 0 0 1/21 -1/21; 1/126 -1/126 0 0 -2/315 2/315; 2/315 -4/315 0 2/315 2/105 -2/105; 4/315 -2/315 0 -2/315 2/105 -2/105; -2/105 -4/315 0 2/63 -2/315 2/315; 4/315 2/105 0 -2/63 -2/315 2/315; 1/126 -1/126 0 0 -2/315 2/315; -4/63 4/63 0 0 16/315 -16/315; -2/315 2/63 0 -8/315 16/315 -16/315; -2/63 2/315 0 8/315 16/315 -16/315; -1/210 -1/126 0 4/315 -2/105 2/105; 2/315 2/105 0 -8/315 -4/315 4/315; 2/315 -4/315 0 2/315 2/105 -2/105; -2/315 2/63 0 -8/315 16/315 -16/315; 4/105 4/63 0 -32/315 16/105 -16/105; -2/315 2/315 0 0 8/105 -8/105; -2/105 -2/315 0 8/315 -4/315 4/315; 1/126 1/210 0 -4/315 -2/105 2/105; 4/315 -2/315 0 -2/315 2/105 -2/105; -2/63 2/315 0 8/315 16/315 -16/315; -2/315 2/315 0 0 8/105 -8/105; -4/63 -4/105 0 32/315 16/105 -16/105]; mateta=[-13/420 0 -1/140 -1/105 1/105 4/105; B.3 Ecuaciones de Navier-Stokes estacionarias 55 1/280 0 11/2520 2/315 -2/315 -1/126; 1/280 0 -1/280 -1/630 1/630 0; -2/105 0 -2/315 4/315 -4/315 8/315; -1/210 0 -1/126 2/105 -2/105 4/315; -2/105 0 -4/315 2/315 -2/315 2/63; 1/280 0 11/2520 2/315 -2/315 -1/126; 1/140 0 -1/140 -1/21 1/21 0; -11/2520 0 -1/280 2/315 -2/315 1/126; 4/315 0 -2/315 -2/105 2/105 -2/315; 2/315 0 -4/315 -2/105 2/105 2/315; 1/126 0 -1/126 2/315 -2/315 0; 1/280 0 -1/280 -1/630 1/630 0; -11/2520 0 -1/280 2/315 -2/315 1/126; 1/140 0 13/420 -1/105 1/105 -4/105; 1/126 0 1/210 2/105 -2/105 -4/315; 2/315 0 2/105 4/315 -4/315 -8/315; 4/315 0 2/105 2/315 -2/315 -2/63; -2/105 0 -2/315 4/315 -4/315 8/315; 4/315 0 -2/315 -2/105 2/105 -2/315; 1/126 0 1/210 2/105 -2/105 -4/315; -4/63 0 -4/105 -16/105 16/105 32/315; -2/315 0 2/315 -8/105 8/105 0; -2/63 0 2/315 -16/315 16/315 8/315; -1/210 0 -1/126 2/105 -2/105 4/315; 2/315 0 -4/315 -2/105 2/105 2/315; 2/315 0 2/105 4/315 -4/315 -8/315; -2/315 0 2/315 -8/105 8/105 0; 4/105 0 4/63 -16/105 16/105 -32/315; -2/315 0 2/63 -16/315 16/315 -8/315; -2/105 0 -4/315 2/315 -2/315 2/63; 1/126 0 -1/126 2/315 -2/315 0; 4/315 0 2/105 2/315 -2/315 -2/63; -2/63 0 2/315 -16/315 16/315 8/315; -2/315 0 2/63 -16/315 16/315 -8/315; -4/63 0 4/63 -16/315 16/315 0]; matxi=reshape(matxi’,size(matxi))’; mateta=reshape(mateta’,size(mateta))’; NLxe=kron(abs(J).*xix,matxi)+kron(abs(J).*etx,mateta); NLye=kron(abs(J).*xiy,matxi)+kron(abs(J).*ety,mateta); I0NL=reshape(T,ne*nt,1)*ones(1,ne^2); J0NL=kron(e’,kron(T’,e))+nn*kron(T’-1,ones(ne)); NLx=sparse(I0NL,J0NL,NLxe); NLy=sparse(I0NL,J0NL,NLye); Código B.9 Resolución ecuaciones estacionarias de Navier-Stokes, caso Scott-Vogelius. pruste2d.m. clear, close all, clc N=12; % número de subdivisiones n bvar=@bwinb; % término de convección uex=@bwinb; % solución exacta force=@forceb; % término independiente epsilon=1e-6; % parámetro de difusión alpha=10; % constante del término de reacción adicional