Repositorio Institucional de Documentos

Universidad de Zaragoza Centro Polit´ ecnico Superior Proyecto Fin de Carrera Ingenier ´ ıa Industrial Estimaci´on de error local para el m´etodo de elementos finitos aplicado a las ecuaciones de Euler en un flujo alrededor de un perfil aerodin´amico Autor: Fernando Lizarraga Roncal Director: Guillermo Hauke Bernardos ´ Area de mec´anica de Fluidos Curso 2010/2011 Abril 2011 Estimaci´on de error local para el m´etodo de elementos finitos aplicado a las ecuaciones de Euler en un flujo alrededor de un perfil aerodin´amico Resumen Este proyecto investiga t´ecnicas de estimaci´on de error a posteriori para el m´etodo de elementos finitos, aplicadas a un caso de inter´es en la industria aeron´autica. El estimador de error objeto del estudio utiliza el m´etodo variacional de las multiescalas. La ventaja de este m´etodo reside en que se evita la resoluci´on de ecuaciones diferenciales para estimar el error. Este m´etodo ha resultado eficiente en el estudio de ondas de choque en flujos supers´onicos. En el presente proyecto, se pretende analizar el comportamiento de este estimador de error en r´egimen subs´onico para las ecuaciones de Euler. Para este tipo de estudios es necesario comparar el error estimado con el error real. La elecci´on de los perfiles de Joukowski para este estudio se justifica no solo en su relevancia en la historia de la aeron´autica, sino en que adem´as, haciendo uso de una transformaci´on conforme en variable compleja en combinaci´on con el m´etodo de Poggi, tienen una soluci´on anal´ıtica conocida para flujo potencial compresible. Se ha estudiado el estimador de error con respecto a gran cantidad de par´ametros como pueden ser, el n´umero de Mach del flujo, las matrices de tiempos intr´ınsecos, tando del flujo como del estimador de error y la norma elegida para el c´alculo del residuo. Tambi´en se han usado distintas geometr´ıas para asegurar la repetibilidad de los resultados. ´ Indice general 1. Introducci´on 1 1.1. Estimaci´on de error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2. Objetivos .................................... 3 1.3. Metodolog´ıa................................... 4 2. M´etodo Variacional de las Multiescalas 6 2.1. Formulaci´on matem´atica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.1.1. Formulaci´on d´ebil y formulaci´on fuerte del problema . . . . . . . . 7 2.1.2. Estimaci´on de error por el VMSEE . . . . . . . . . . . . . . . . . . 10 3. Perfil de Joukowski 12 3.1. Transformaci´on conforme . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 3.2. Transformaci´on de Joukowski . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.2.1. Obtenci´on de las coordenadas del perfil . . . . . . . . . . . . . . . . 15 3.2.2. Obtenci´on de la soluci´on anal´ıtica para flujo potencial . . . . . . . . 16 3.3. Modelo computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3.1. Creaci´on de la geometr´ıa . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.2. sec:bc .................................. 19 3.3.3. Mallado ................................. 20 4. Validaci´on y limitaciones 22 4.1. Validaci´on del c´odigo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1.1. Error en la condici´on de contorno . . . . . . . . . . . . . . . . . . . 23 4.2. Error de convergencia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 5. Resultados 26 5.1. Presentaci´on de los resultados . . . . . . . . . . . . . . . . . . . . . . . . . 26 5.2. Estudiob´asico.................................. 28 5.2.1. Repetibilidad o consistencia . . . . . . . . . . . . . . . . . . . . . . 28 5.2.2. Dependencia con el n´umero mach . . . . . . . . . . . . . . . . . . . 29 5.2.3. Dependencia con respecto a los par´ametros del m´etodo . . . . . . . 32 5.3. Estudio en profundidad del caso base . . . . . . . . . . . . . . . . . . . . . 34 6. Conclusiones y posibles trabajos futuros 37 Bibliograf´ıa 40 i Anexos 41 A. M´etodo de Poggi 42 B. Adimensionalizaci´on del problema 51 C. Figuras suplementarias 53 ii Cap´ıtulo 1 Introducci´on Los m´etodos computacionales son una de las ´ultimas herramientas a disposici´on del ingeniero especializado en mec´anica de fluidos. Se empezaron a desarrollar a mediados del siglo XX tras la irrupci´on de los computadores. En un principio se limitaron a apoyar o a extender resultados emp´ıricos. Hoy en d´ıa son capaces de ofrecer soluciones confiables a modelos f´ısicos plasmados en ecuaciones diferenciales y a pesar de su relativa novedad, se han integrado de una manera muy notable en los m´as variados procesos industriales. El r´apido desarrollo de los ordenadores1ha supuesto un aumento de su capacidad de c´alculo al mismo tiempo que un abaratamiento de su coste. La mec´anica de fluidos es un campo donde los m´etodos computacionales son de gran utilidad y al mismo tiempo suponen un gran reto. Ambas situaciones se deben a la gran complejidad de los modelos matem´aticos que describen el comportamiento de los fluidos. Hay que recordar que las soluciones anal´ıticas que existen en mec´anica de fluidos o bien sirven para geometr´ıas muy sencillas, o bien asumen hip´otesis simplificativas que desvirt´uan la soluci´on obtenida. Sin embargo, Los m´etodos de CFD2se adaptan a cualquier geometr´ıa y permiten abordar el problema con toda su complejidad a cambio de tiempo de computaci´on. Tambi´en pueden ser una poderosa herramienta de dise˜no en campos como la aerodin´amica externa de veh´ıculos, ya que permiten que s´olo en las ´ultimas iteraciones del dise˜no sea necesario el uso de maquetas o prototipos para ensayar en t´uneles aerodin´amicos, produci´endose as´ı un ahorro considerable en costes. 1La ley de Moore dice que cada dos a˜nos se duplica la cantidad de transistores que se pueden integrar por unidad de superficie, sin un aumento significativo del coste del dispositivo. Esta ley se ajusta prefectamente a los datos hist´oricos desde los a˜nos 70 hasta nuestros d´ıas. Se prev´e que siga cumpli´endose hasta el a˜no 2015. 2Fluidodin´amica Computacional 1 CAP´ ITULO 1. INTRODUCCI´ ON 1.1. ESTIMACI´ ON DE ERROR 1.1. Estimaci´on de error La estimaci´on de error es uno de los campos de investigaci´on de mayor inter´es en CFD. Hoy en d´ıa, los m´etodos computacionales son fiables, ´utiles y cada vez m´as r´apidos, de tal manera que ya est´an plenamente integrados en el proceso productivo. La estimaci´on de error surge como una mejora natural de los mismos. Como en otros campos de la ciencia o de la ingenier´ıa, es deseable conocer con qu´e margen de error se ajustan a la realidad los c´alculos o incluso los resultados de un experimento. En los m´etodos experimentales, una medida no se considera completa sin una estimaci´on del error cometido. La tendencia en los m´etodos CFD deber´ıa ser la misma, ya no s´olo por una cuesti´on de rigor cient´ıfico, sino porque, como veremos m´as adelante, una correcta evaluaci´on del error puede suponer un ahorro considerable en el coste computacional (y por tanto econ´omico) de las simulaciones. Las soluciones obtenidas por los m´etodos computacionales en general y el M´etodo de los elementos finitos en particular, dependen fuertemente de la discretizaci´on, ya sea del dominio espacial (geometr´ıa) o temporal. En este trabajo nos centraremos en el M´etodo de los elementos finitos y en la influencia de la discretizaci´on espacial. Se sabe que para una malla infinitamente densa, la soluci´on por elementos finitos converge a la soluci´on exacta del problema. En general, cuanto m´as densa es la discretizaci´on, mejor es la soluci´on obtenida, aunque su coste computacional es mayor. La elecci´on de la malla surge por tanto del compromiso entre la precisi´on de la soluci´on requerida y el coste computacional que se est´a dispuesto a asumir. La estrategia de mallado m´as interesante es aquella en la que se concentran mayor n´umero de elementos en la zona donde estos son necesarios. As´ı, se impone mayor resoluci´on en las zonas donde los gradientes son intensos. Cabe destacar que obtener el mallado adecuado ser´ıa trivial si tuviesemos ya una soluci´on donde se observen todos los fen´omenos relevantes. Es en este punto donde tradicionalmente ha entrado en juego la intuici´on y experiencia de la persona que est´a resolviendo flujo, que observando con atenci´on los resultados, tras un proceso iterativo que muchas veces es largo y tedioso, puede obtener una malla adecuada para la resoluci´on del problema. Las t´ecnicas de estimaci´on de error pueden desempe˜nar un papel crucial en esta situaci´on. Para esto es necesario que el estimador de error sea local. Si se conoce el error de la simulaci´on en cada punto con fiabilidad, y se relaciona dicho error con el tama˜no de los elementos de la malla, se podr´ıa elegir un umbral de error que se considere aceptable y elaborar una malla que obtenga ese resultado. Actuando as´ı se obtendr´ıan simulaciones con un error controlado aumentando la densidad de la malla s´olo donde sea necesario. A esta estrategia se le llama Mallado adaptativo. 2 CAP´ ITULO 1. INTRODUCCI´ ON 1.2. OBJETIVOS x y 0 0.5 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (a) Malla inicial (b) Malla refinada Figura 1.1: Mallado adaptativo Se observa en la Figura 1.1 c´omo se ha pasado de una malla inicial bastante uniforme y tosca a una bastante m´as refinada. En esta ´ultima, el tama˜no de elementos no es uniforme, intuy´endose la presencia de varias capas l´ımite muy bien resueltas por el mallado adaptativo. Es curioso comprobar el hecho de que una malla bien resuelta puede servir para “adivinar” la soluci´on del problema. Este trabajo estudia Estimadores de Error basados en el M´etodo Variacional de las Multiescalas.3A grandes rasgos, obtiene una aproximaci´on del error local cometido tomando como entradas una malla de elementos finitos y los resultados correspondientes a la simulaci´on. Ya se ha demostrado su fiabilidad para problemas de flujo de Euler bidimensional compresible y supers´onico. En este proyecto se pretende estudiar su comportamiento en flujo compresible subs´onico. Adem´as tiene un coste computacional muy bajo. Esta ´ultima cuesti´on es crucial, ya que existen m´etodos de estimaci´on de error con un coste superior al de la computaci´on del problema. A lo largo del documento se expondr´an las caracter´ısticas de este estimador. 1.2. Objetivos El objetivo de este proyecto de fin de carrera es evaluar y caracterizar un estimador de error para el m´etodo de los elementos finitos basado en el M´etodo variacional de las Multiescalas, estudiando un caso de flujo de Euler subs´onico. La evaluaci´on consiste en la comparaci´on del error estimado por el VMSEE con el error exacto. Es conveniente hacer notar en este punto que el estimador de error da como resultado una norma del error en cada elemento. En sucesivos apartados se expondr´a adecuadamente la definici´on de este residuo en el contexto de los elementos finitos. Tambi´en se busca caracterizar el VMSEE estudiando su comportamiento frente a la variaci´on de diversos par´ametros, tanto extr´ınsecos, es decir, relacionados con el flujo que se estudia; como intr´ınsecos o relacionados con las variables del m´etodo. 3A partir de ahora VMSEE 3 CAP´ ITULO 1. INTRODUCCI´ ON 1.3. METODOLOG´ IA Los par´ametros que se van a estudiar son: Dependencia con el n´umero de Mach. Dependencia con respecto al tama˜no del elemento. Dependencia de la eficiencia del error con el error exacto cometido Se estudiar´an distintas geometr´ıas, en este caso perfiles de la familia de Joukowski, para asegurar la repetibilidad de los resultados Dependencia con respecto a la matriz de tiempos intr´ınsecos del m´etodo estabilizado Dependencia con respecto a la norma elegida para el c´alculo del residuo 1.3. Metodolog´ıa Para evaluar el estimador de error es necesario establecer un patr´on de comparaci´on para saber c´omo de bueno es el m´etodo. En nuestro caso se ha elegido una simulaci´on de flujo de Euler en condiciones subs´onicas alrededor de un perfil aerodin´amico, el perfil de Joukowski. Este ´ultimo presenta la ventaja de tener una soluci´on anal´ıtica bastante f´acil de calcular y que es v´alida para n´umeros de Mach bajos en flujo subs´onico. Con objeto de poder explicar c´omo se ha desarrollado el proyecto, es necesario adelantar la definici´on de la eficiencia del estimador de error. Sea la descomposici´on de la soluci´on: u= ¯u+u′ siendo ula soluci´on exacta del problema, ¯ula soluci´on obtenida por el m´etodo de elementos finitos y u′el error exacto. El VMSEE estima el valor de la norma del error ηe. As´ı la eficiencia del estimador de error Ie eff se define como Ie eff =ηe ku−¯ukΩe donde ηees la norma del error estimado para un elemento y ku−¯ukΩees la norma del error exacto. 4 CAP´ ITULO 1. INTRODUCCI´ ON 1.3. METODOLOG´ IA La estrategia de investigaci´on en el proyecto se ha desarrollado de acuerdo con los siguientes pasos: Creaci´on de un perfil de Joukowski. Obtenci´on de la soluci´on exacta del flujo alrededor de este perfil. Creaci´on de una malla de elementos finitos y soluci´on del problema. C´alculo de las eficiencias de acuerdo con la definici´on anterior. Postproceso y an´alisis de los datos. 5 Cap´ıtulo 3 Flujo alrededor de un perfil de Joukowski Como se ha expuesto en el cap´ıtulo anterior, para el desarrollo del proyecto ser´a necesario comparar la soluci´on de un problema de flujo de Euler subs´onico alrededor de un perfil alar mediante el M´etodo de los Elementos finitos, con la soluci´on exacta del mismo, a fin de obtener el error que se est´a cometiendo. En este proyecto se utilizar´a la soluci´on para un perfil de Joukowski. Figura 3.1: Partes de un perfil sim´etrico Poggi Cilindro incompresible Cilindro compresible Joukowski Perfil Compresible Figura 3.2: Obtenci´on de la soluci´on exacta La transformaci´on de Joukowski es un mapeado conforme que permite trasladar la soluci´on de un flujo potencial alrededor de un cilindro a la de una geometr´ıa con forma de perfil alar. La soluci´on del flujo ideal (no viscoso e incompresible) alrededor del cilindro es bien conocida. Mediante el m´etodo de Poggi (Anexo A), se obtiene una soluci´on para flujo de Euler compresible. Dicha soluci´on tambi´en deriva de un potencial, por lo que aplicando sobrrta ella la transformaci´on de Joukowski, se obtiene el flujo alrededor del perfil. 12 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.1. TRANSFORMACI´ ON CONFORME 3.1. Transformaci´on conforme Desde un punto de vista matem´atico, la transformaci´on de Joukowski es una aplicaci´on o mapeado conforme. Como tal, se caracteriza porque es un cambio de coordenadas que conserva los ´angulos. Se suele estudiar en an´alisis complejo ya que, como veremos a continuaci´on, presenta propiedades muy ´utiles. Sea Uun conjunto abierto perteneciente al plano complejo C. Se dice que: w(z): U∈C−→ CEs una aplicaci´on conforme si y s´olo si ∃dw dz ydw dz 6= 0 (3.1) Tiene una propiedad muy interesante con respecto a las funciones arm´onicas. Dada la φ(x, y), se dice que es arm´onica si cumple la ecuaci´on de Laplace: d2φ dx2+d2φ dy2= 0 (3.2) O de manera m´as compacta: ∇2φ= 0 (3.3) Se cumple que la imagen a trav´es de una aplicaci´on conforme de una funci´on arm´onica, es tambi´en arm´onica. Sea z=z(w) un mapeado conforme que transforma el dominio D del plano z, en un dominio D1en el plano w. Sea φ1(g, h) una funci´on arm´onica en D1: d2φ1 dg2+d2φ1 dh2= 0 (3.4) Entonces, aplicando la transformaci´on dada por: z=z(w) = g(x, y) + ih(x, y) (3.5) Tenemos que φ(x, y) = φ1(g(x, y), h(x, y)) es tambi´en arm´onica en D, es decir: d2φ dx2+d2φ dy2= 0 Esta caracter´ıstica es particularmente ´util para resolver la ecuaci´on de Laplace,ya que permite transformar un problema de geometr´ıa muy complicada en otro de geometr´ıa m´as sencilla cuya soluci´on sea m´as asequible. Viendo la Figura 3.3, se puede observar como se transforman las condiciones de contorno. 13 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.2. TRANSFORMACI ´ ON DE JOUKOWSKI D D 1 Figura 3.3: Transformaci´on conforme. Condiciones de Dirichlet: (g0, h0)←→ (x0, y0) φ1(g0, h0) = k←→ φ(x0, y0) = k Condiciones de Neumann: (g1, h1)←→ (x1, y1) dφ dn = 0 ←→ dφ1 dn1 = 0 Donde dφ dn es la derivada en direcci´on normal a un contorno dado. 3.2. Transformaci´on de Joukowski El m´etodo para obtener la soluci´on anal´ıtica consiste en proyectar un perfil aerodin´amico del plano f´ısico z=x+ iy, sobre un cilindro en un plano ficticio w=g+ ih, para el que se conoce la soluci´on del flujo potencial. La ecuaci´on que rige el flujo potencial bidimensional es la ecuaci´on arm´onica. Por lo tanto, debido a las propiedades de la transformaci´on conforme, podremos usar la soluci´on del cilindro para obtener el campo de velocidades del perfil alar. En este apartado se explicar´a el proceso de obtenci´on, tanto de las coordenadas de un perfil de Joukowski, como del campo de velocidades que le corresponde. Cabe destacar que existen varias versiones de la transformaci´on de Joukowski. La que se va a usar a lo largo del proyecto es la que propone Katz y Plotkin [1991], porque utiliza 14 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.2. TRANSFORMACI ´ ON DE JOUKOWSKI como par´ametros dos valores que se aproximan a la espesor m´aximo y a la cuerda del perfil respectivamente. b (a) Plano del cilindro (b) Plano del perfil Figura 3.4: Transformaci´on de Joukowski Se observa en la Figura 3.4 que se est´a tratando un caso particular sim´etrico del perfil de Joukowski. La velocidad aguas arriba del perfil es paralela a la cuerda del mismo. La transformaci´on de Joukowski queda descrita por las siguientes ecuaciones z=w+a2 16w µ=−ǫa/4 b=a 4(1 + ǫ) (3.6) a: Aproximadamente la cuerda del perfil. ǫ: Aproximadamente espesor m´aximo del perfil. b: Radio de la circunferencia. µ: Coordenada del centro de la circunferencia en el eje real. 3.2.1. Obtenci´on de las coordenadas del perfil Aplicando la f´ormula 3.6 a los puntos de un cilindro de radio by centro (−µ, 0), obtendremos un perfil de cuerda aproximadamente ay espesor m´aximo aproximadamente ǫ. La Figura 3.5 representa la geometr´ıa del perfil que m´as se va a utilizar a lo largo de este estudio. Estas coordenadas se introducir´an en el preprocesador para obtener una malla que resolver por elementos finitos. Eligiendo un perfil sim´etrico sumergido en un flujo paralelo 15 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.2. TRANSFORMACI ´ ON DE JOUKOWSKI -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 (a) r= 0.275, µ =−0.025 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 (b) C≈1.0, ǫ ≈0.1 Figura 3.5: Caso base a la cuerda, se consiguen varias ventajas. En primer lugar, el coste computacional de la simulaci´on es la mitad, y en segundo lugar, se evita tener que imponer la condici´on de Kutta. 3.2.2. Obtenci´on de la soluci´on anal´ıtica para flujo potencial A continuaci´on, se obtendr´a una soluci´on anal´ıtica para el campo de velocidades del perfil de Joukowski. Sea el potencial complejo en el plano del c´ırculo φ(w) y su velocidad compleja U(w). Los resultados en el plano f´ısico son: φ(z) = φ[z(w)] U(z) = dφ dz =dφ dw dw dz Como U(w) = dφ dw (3.7) tenemos que U(z) = U(w)1 dz dw (3.8) La soluci´on del flujo potencial alrededor de un cilindro U(w) se obtiene mediante el m´etodo de Poggi, que se detalla en el Anexo A. Mediante este m´etodo se obtiene el perfil de velocidades para flujo compresible partiendo del perfil de velocidades para flujo incompresible, que es ampliamente conocido. Por tanto, para obtener la velocidad compleja U(z) de un punto del plano f´ısico de coordenadas z=x+ iy, es necesario calcular sus coordenadas w=g+ ihen el plano del cilindro. Para ello ser´a necesario calcular la transformaci´on inversa a la ecuaci´on 3.6. A trav´es de sencillas manipulaciones algebraicas obtenemos la ecuaci´on 3.9.    w(z) = −qz2−a2 4 2,si x≥0 w(z) = +qz2−a2 4 2,si x≤0 (3.9) 16 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.3. MODELO COMPUTACIONAL Por las propiedades del flujo potencial, la componente horizontal de la velocidad u(x, y) = Re(U(z)) y la vertical v(x, y) = −Im((U(z)). El campo de presiones p(x, y) se obtiene con la ecuaci´on de Bernouilli. x y -0.6 -0.4 -0.2 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 u1 1.1 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 (a) u x y -0.6 -0.4 -0.2 0 0.2 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 u2 0.6 0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 -0.05 (b) v Figura 3.6: Campo de velocidades en las proximidades del perfil El campo de velocidades se puede observar en la Figura 3.6. 3.3. Modelo computacional En este apartado se pretende ilustrar c´omo se obtuvo la soluci´on al flujo de Euler subs´onico mediante el m´etodo de los elementos finitos. En el Cap´ıtulo 2 se describe brevemente la ecuaci´on de Euler y se hace una breve introducci´on te´orica al M´etodo de los Elementos finitos. La soluci´on de un problema por elementos finitos consta de las siguientes fases: Creaci´on de la geometr´ıa. Imposici´on de las condiciones de contorno. Mallado. Elecci´on de los par´ametros de la simulaci´on. Simulaci´on. Para realizar las tres primeras fases se ha utilizado el programa comercial GiD. GiD es un preprocesador multiprop´osito para el m´etodo de elementos finitos. Proporciona una infraestructura para el dise˜no de la geometr´ıa y para el mallado. Permite programar, tan- to las condiciones de contorno necesarias para resolver el problema, como la salida del programa. As´ı, se adapta a cualquier c´odigo de elementos finitos. En nuestro caso, se han programado las condiciones de contorno de Euler el interfaz con el “solver” ENSA. Posteriormente se realiza el postproceso de los datos, consistente en transformar los datos obtenidos en informaci´on ´util. Como se ha comentado en el Cap´ıtulo 2, los resultados 17 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.3. MODELO COMPUTACIONAL de las simulaciones se utilizan como entrada para el estimador de error y para calcular la norma del error exacto, obteni´endose as´ı la eficiencia del estimador. 3.3.1. Creaci´on de la geometr´ıa Consiste en modelar el dominio fluido donde se va a aplicar el m´etodo de elementos finitos. A B C DD a Rdom Figura 3.7: Geometr´ıa y condiciones de contorno En caso de flujo de Euler alrededor de un obst´aculo, es necesario mallar una gran porci´on del espacio que lo rodea. Esto es debido a que la ausencia de fen´omenos difusivos hace que las perturbaciones aguas abajo del perfil tarden mucho en desaparecer, lo que fuerza a que las condiciones de contorno de salida (B) tengan que estar lejos. En nuestro caso se ha elegido un radio Rdom = 300 para una cuerda a= 1, ya que a esas distancias el flujo se puede considerar imperturbado. Dado que la velocidad aguas arriba del perfil es paralela a la cuerda del mismo, se tiene que tanto la geometr´ıa, como las condiciones de contorno, presentan simetr´ıa con respecto al eje x. Por tanto, s´olo es necesario simular la mitad del dominio. La geometr´ıa del perfil se ha modelado interpolando mediante NURBS1los puntos del perfil de Joukowski obtenidos por el procedimiento descrito en la secci´on 3.2.1. Este hecho introduce una fuente de error que se comentar´a en el Cap´ıtulo 4. 1NURBS (acr´onimo ingl´es de la expresi´on Non Uniform Rational B-splines) es un modelo matem´atico muy utilizado en la computaci´on gr´afica para generar y representar curvas y superficies. 18 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.3. MODELO COMPUTACIONAL 3.3.2. Imposici´on de las condiciones de contorno Como se observa en la Figura 3.7, el dominio esta delimitado por cuatro l´ıneas. Para resolver el flujo, es necesario imponer unas condiciones de contorno que determinen el problema. Entrada(A) ρ∞= 1.0 u∞= 1.0 v∞= 0.0 Salida (B) p∞ Perfil (C) −→ v·−→ n= 0 Flujo m´asico nulo. Eje de simetr´ıa (D) v= 0 Flujo m´asico nulo. Una vez tenemos la situaci´on de las condiciones de contorno que se van a imponer, necesitamos saber el valor num´erico de las mismas. Se sabe por an´alisis dimensional, que el flujo de Euler alrededor de un perfil alar depende del n´umero de Mach, My de la relaci´on de calores espec´ıficos γ(ver Anexo B). En nuestro caso, siempre vamos a simular aire, por tanto γ= 1.4. El flujo queda entonces perfectamente caracterizado mediante el n´umero de Mach. M=U∞ c ¯ R=R W c=pγ¯ RT p=ρ¯ RT Siendo: c: Velocidad del sonido en el gas. R: Constante de los gases. R= 8.314 kJ kmol·K U: Velocidad. W: Peso molecular del aire. W= 28.97 kg kmol p: Presi´on. ρ: Densidad. T: Temperatura. 19 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.3. MODELO COMPUTACIONAL Manteniendo fijas U∞= 1.0 y ρ∞= 1.0, para cada valor de M, quedan determinados los valores de pyT(Tabla 3.1). Por tanto, controlaremos el n´umero de Mach de la simulaci´on modificando la presi´on aguas abajo del perfil. La condici´on sobre la temperatura se usar´a como condici´on inicial.2 M P∞T∞ 0.1 7.14285 ×1012.4857 ×10−1 0.2 1.78571 ×1016.1943 ×10−2 0.3 7.93650 ×1002.7525 ×10−2 Tabla 3.1: Mach, temperatura y presi´on 3.3.3. Mallado El mallado es simplemente la discretizaci´on del dominio fluido donde queremos resolver el problema. En el cap´ıtulo de introducci´on, se ha hablado con cierta profundidad de la relevancia del mismo. En nuestro caso, se han seguido los siguientes criterios para realizar la malla. El mallado del perfil, particularmente en el borde de ataque, ha de ser extremadamente denso. La influencia de esta zona es muy importante en los resultados de la simulaci´on, un mal mallado no solo traer´ıa consigo los problemas inherentes a un tama˜no grande de elemento, sino que distorsiona la condici´on de contorno que se est´a imponiendo. Sobre esta cuesti´on, se abundar´a en el Cap´ıtulo 4. La densidad de mallado ir´a disminuyendo gradualmente conforme nos alejamos del perfil por dos razones. En primer lugar, para tener distintos tama˜nos de elemento en la malla, pudiendo hacer as´ı estudiar la influecia del tama˜no en la eficiencia y por supuesto, para disminuir el coste computacional de la simulaci´on. En GiD se pueden asignar tama˜nos de elemento a puntos, l´ıneas y superficies, as´ı como regular las transici´on entre distintos tama˜nos de elemento. Se ha asignado un tama˜no extremadamente peque˜no al borde de ataque(tba = 0.0005), al perfil se le ha dado un tama˜no algo m´as grande (tp= 0.0005) y para el resto de la malla se han dejado crecer libremente los elementos con transiciones no demasiado forzadas. Para las simulaciones se han elegido cuadrilateros no estructurados bilineales. El aspecto de la malla es el siguiente: 2A priori puede chocar que se hable de condici´on inicial en un flujo estacionario. La raz´on es que la implementaci´on del m´etodo de elementos finitos que se utiliza determina los estacionarios a trav´es de la simulaci´on del flujo transitorio. 20 CAP´ ITULO 3. PERFIL DE JOUKOWSKI 3.3. MODELO COMPUTACIONAL x y -300 -200 -100 0 100 200 300 0 100 200 300 (a) Malla general x y -1 -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 (b) Detalle perfil. x y -0.5 -0.49 0 0.005 0.01 (c) Detalle borde de ataque. Figura 3.8: Mallado 21 CAP´ ITULO 5. RESULTADOS 5.2. ESTUDIO B ´ ASICO x y -300 -200 -100 0 100 200 300 0 50 100 150 200 250 300 error_P 1E-05 1E-06 1E-07 (a) Dominio completo x y -0.5 0 0.5 0 0.2 0.4 0.6 0.8 error_P 0.005 0 -0.005 -0.01 -0.015 -0.02 -0.025 (b) Detalle perfil. Figura 5.1: Error exacto en P. M= 0.1,ǫ= 0.15,τf= 1 y τe= 5 A lo largo de todos los casos estudiados, la variable que peor comportamiento ha mostrado ha sido la velocidad, ya que se ha observado que en las zonas m´as cercanas al obst´aculo el VMSEE tiende a subestimar el error. En este cap´ıtulo la mayor´ıa de los ejemplos ser´an presentados sobre esta variable. 5.2. Estudio b´asico A continuaci´on se explican todos los casos que se han tratado en el proyecto. Entendemos como un caso distinto la combinaci´on de una simulaci´on de elementos finitos; variaci´on de geometr´ıa, Mach y τfy un conjunto de par´ametros del VMSEE; norma, τey modelo de estimador. Teniendo en cuenta esto, se han estudiado 486 casos, por lo que forzosamente se tendr´an que presentar los m´as representativos. VMSEE FEM τerror Modelo Norma Mach Geometr´ıa τflow τ1Standard L1 L2 0.1 Cilindro,R= 5 τ1 τ5Naive 0.2 Perfil, ǫ= 0.1τ5 τ9Upper bound 0.3 Perfil, ǫ= 0.2τ9 Tabla 5.2: Casos estudiados 5.2.1. Repetibilidad o consistencia Cuando se desarrolla un m´etodo num´erico es necesario estudiar si funciona para cualquier caso o si por contra, s´olo funciona en un n´umero limitado de casos. En este trabajo, se eval´ua el desempe˜no del VMSEE como estimador de error para flujo de Euler compresible alrededor de un obst´aculo, con objeto de determinar si ser´ıa un m´etodo ´util para estudios aerodin´amicos industriales. Para asegurar la repetibilidad, se ha simulado el flujo alrededor de tres geometr´ıas distintas que podr´ıan ser un caso industrial. Dos de ellas son perfiles de Joukowski y la tercera es un cilindro de secci´on circular. Como ya se ha se˜nalado en el Cap´ıtulo 2, poseen una soluci´on anal´ıtica para flujo de Euler. Se podr´ıan haber elegido otras geometr´ıas que no ofreciesen las dificultades en cuanto a la condici´on de contorno se˜naladas en la secci´on 28 CAP´ ITULO 5. RESULTADOS 5.2. ESTUDIO B ´ ASICO 4.2, como flujo alrededor de esquinas u otros con soluci´on anal´ıtica, pero estos no se ajustan tan bien a los casos cl´asicos en aerodin´amica de flujo alrededor de obst´aculos. Figura 5.2: Cuerda ay espesor ǫ Cuerda Espesor Radio Cilindro - - 0.5 Perfil 1.0 0.10 - 1.0 0.15 - 1.0 0.20 - Tabla 5.3: Geometr´ıas de estudio A la hora de evaluar los resultados de estas pruebas, ser´ıa necesario fijar un criterio de rechazo de los test, es decir, determinar cuando se considera que las diferencias entre el comportamiento del estimador de error para las distintas geometr´ıas son distintas. En nuestro caso el criterio ser´a cualitativo, de tal manera que se evaluar´an estos dos aspectos: Similitud de la distribuci´on frecuencial de la eficiencia. Similitud de los mapas de color. En los mapas de color de la Figura 5.3 se observa que hay problemas de baja eficiencia en los primeros elementos cercanos a la condici´on de contorno para ambas geometr´ıas. Esto es m´as o menos esperable ya que, como se ha comentado en el Cap´ıtulo 4, el VMSEE no tiene en cuenta los errores en la condici´on de contorno. En el histograma de la Figura 5.4, se observa que no hay graves disparidades en los resultados. Es importante hacer notar una vez m´as que se ha elegido un caso representativo de los muchos que se han estudiado. En todos ellos se manifiesta este tipo de comportamiento. Puede resultar llamativo que no se usen las eficiencias globales, pero se ha observado que en algunos casos, el error incluso en un elemento altera notablemente el resultado. 5.2.2. Dependencia con el n´umero mach Como ya se ha comentado en Anexo B, gracias al an´alisis dimensional sabemos que una vez fijado que el fluido es aire, el problema depende ´unicamente del n´umero de Mach del flujo imperturbado. Por consiguiente, es interesante estudiar el comportamiento del VMSEE si variamos este par´ametro. 29 CAP´ ITULO 5. RESULTADOS 5.2. ESTUDIO B ´ ASICO x y -1 -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 errL2-u 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (a) Cilindro x y -1 -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 errL2-u 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (b) Perfil ǫ= 0.15 Figura 5.3: Eficiencia de la velocidad u, norma L2,M= 0.1, modelo standard Seg´un se puede comenta en el Cap´ıtulo 3 y tal y como se desarrolla en el Anexo A, contamos con soluciones para flujo compresible bastante exactas para n´umeros de Mach bajos. Se ha simulado para M= 0.1, M= 0.2 y M= 0.3. Para poner en contexto esta limitaci´on, hay que tener en cuenta que a temperatura ambiente y presi´on atmosf´erica, a esos n´umeros de Mach corresponde una velocidad entre los 125 y los 370 Km/h. 30 CAP´ ITULO 5. RESULTADOS 5.2. ESTUDIO B ´ ASICO 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0−0.1 0.1−0.5 0.5−2.0 2.0−5 5.0 Frecuencia Buckets(Eficiencia) Círculo Perfil ε=0.10 Perfil ε=0.15 Figura 5.4: Eficiencia L2de la velocidad u,M= 0.1, τf= 1, τe= 5 para cada geometr´ıa. Radio=2.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0−0.1 0.1−0.5 0.5−2.0 2.0−5 5.0 Frecuencia Buckets(Eficiencia) M=0.1, p M=0.2, p M=0.3, p M=0.1, u M=0.2, u M=0.3, u M=0.1, v M=0.2, v M=0.3, v Figura 5.5: Eficiencias L2de la presi´on para los distintos Mach, modelo standard 31 CAP´ ITULO 5. RESULTADOS 5.2. ESTUDIO B ´ ASICO 5.2.3. Dependencia con respecto a los par´ametros del m´etodo En el Cap´ıtulo 2 se describe someramente el VMSEE y los par´ametros de los que depende, es decir, matrices de escala temporal, dise˜no del estimador y norma utilizada para el c´alculo del residuo. Se han combinado cada uno de estos factores para estudiar la eficiencia de cada uno de los estimadores resultantes. Matrices de tiempos intr´ınsecos τf τf τfyτe τe τe Como se ha comentado en el Cap´ıtulo 2, la matriz de escalas temporales es un par´ametro presente tanto en la formulaci´on del m´etodo de elementos finitos (τf) como en la del estimador de error (τe). As´ı pues, se plantea el estudio de su influencia combin´andolas por pares para el flujo y el error. Se estudiar´an tres τdistintas, la cl´asica τ1y dos dise˜nadas para flujos con Mach bajo, τ9yτ5; teniendo as´ı nueve combinaciones para cada norma, como se observa en la Figura 5.6. 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0-0.1 0.1-0.5 0.5-2.0 2.0-5 5.0 Frecuencia Buckets(Eficiencia) τf=1,τe=1 τf=9,τe=1 τf=5,τe=1 τf=1,τe=9 τf=9,τe=9 τf=5,τe=9 τf=1,τe=5 τf=9,τe=5 τf=5,τe=5 (a) p 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0-0.1 0.1-0.5 0.5-2.0 2.0-5 5.0 Frecuencia Buckets(Eficiencia) τf=1,τe=1 τf=9,τe=1 τf=5,τe=1 τf=1,τe=9 τf=9,τe=9 τf=5,τe=9 τf=1,τe=5 τf=9,τe=5 τf=5,τe=5 (b) u 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0-0.1 0.1-0.5 0.5-2.0 2.0-5 5.0 Frecuencia Buckets(Eficiencia) τf=1,τe=1 τf=9,τe=1 τf=5,τe=1 τf=1,τe=9 τf=9,τe=9 τf=5,τe=9 τf=1,τe=5 τf=9,τe=5 τf=5,τe=5 (c) v Figura 5.6: Despliegue de casos para todo τeyτf,ǫ= 0.15, modelo standard y M= 0.1,Norma L2 32 CAP´ ITULO 5. RESULTADOS 5.2. ESTUDIO B ´ ASICO Norma Se han usado dos maneras distintas de calcular la norma del residuo: L1=ZΩe |L ¯ Yj Yj Yj|dΩ L2=sZΩe (L¯ Yj Yj Yj)2dΩ Se observa en la Figura 5.7 que en general, la norma L2obtiene eficiencias m´as altas que las de la L1. Modelo de estimador En el Cap´ıtulo 2 se describen tres dise˜nos del estimador de error, naive, cl´asico y upper bound. Tal y como hemos procedido a lo largo de este trabajo, generamos otros tres casos de estudio. 0 0.1 0.2 0.3 0.4 0.5 0.6 0-0.1 0.1-0.5 0.5-2.0 2.0-5 5.0 Frecuencia Buckets (Eficiencia) Standard L1 Standard L2 Naive L1 Naive L2 Upper Tau Med L1 Upper Tau Med L2 Figura 5.7: Eficiencia de la velocidad u,M= 0.1, τf= 1, τe= 5 para los distintos modelos de estimador y normas Se observa que el modelo upper bound funciona bastante bien como cota superior mientras que el naive y el standard intentan ajustar m´as al valor exacto. 33 CAP´ ITULO 5. RESULTADOS 5.3. ESTUDIO EN PROFUNDIDAD DEL CASO BASE x y -1 -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 errL2-u 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (a) Standard x y -1 -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 errL2-u 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (b) Naive x y -1 -0.5 0 0.5 0 0.2 0.4 0.6 0.8 1 errL2-u 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 (c) Upper-bound Figura 5.8: Velocidad u, norma L2,τf= 1, τe= 5. 5.3. Estudio en profundidad del caso base Para este apartado se har´a un estudio en profundidad del caso base para los tres modelos de estimador. Se procede de esta manera porque se ha considerado que el comportamiento como cotas de los tres modelos, (normal<naive<upper bound), puede ser de mucha utilidad a la hora de desarrollar un algoritmo de remallado. Geometr´ıa Joukowski ǫ= 0.15 Mach 0.1 Norma L2 τfτ1 τeτ5 Tabla 5.4: Caso base Seg´un se coment´o en el Cap´ıtulo primero, se va a estudiar la influencia de diversos factores sobre la eficiencia de cada elemento. Se estudiar´an principalmente dos: La influencia del tama˜no del elemento. La influencia de la norma del error exacto. 34 CAP´ ITULO 5. RESULTADOS 5.3. ESTUDIO EN PROFUNDIDAD DEL CASO BASE Influencia del tama˜no del elemento 1e-05 0.0001 0.001 0.01 0.1 1 10 100 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (a) Standard 0.0001 0.001 0.01 0.1 1 10 100 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (b) Naive 0.0001 0.001 0.01 0.1 1 10 100 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (c) Upper-bound Figura 5.9: Eficiencias de p 0.0001 0.001 0.01 0.1 1 10 100 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (a) Standard 0.0001 0.001 0.01 0.1 1 10 100 1000 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (b) Naive 0.001 0.01 0.1 1 10 100 1000 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (c) Upper-bound Figura 5.10: Eficiencias de u 0.0001 0.001 0.01 0.1 1 10 100 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (a) Standard 0.001 0.01 0.1 1 10 100 1000 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (b) Naive 0.001 0.01 0.1 1 10 100 1000 1e-10 1e-08 1e-06 0.0001 0.01 1 100 10000 L2(Eficiencia) Superficie elemento (c) Upper-bound Figura 5.11: Eficiencias de v 35 CAP´ ITULO 5. RESULTADOS 5.3. ESTUDIO EN PROFUNDIDAD DEL CASO BASE Influencia del tama˜no del norma del error de la variable 1e-05 0.0001 0.001 0.01 0.1 1 10 100 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (a) Standard 0.0001 0.001 0.01 0.1 1 10 100 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (b) Naive 0.0001 0.001 0.01 0.1 1 10 100 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (c) Upper-bound Figura 5.12: Eficiencias de p 0.0001 0.001 0.01 0.1 1 10 100 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (a) Standard 0.0001 0.001 0.01 0.1 1 10 100 1000 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (b) Naive 0.001 0.01 0.1 1 10 100 1000 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (c) Upper-bound Figura 5.13: Eficiencias de u 0.0001 0.001 0.01 0.1 1 10 100 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (a) Standard 0.001 0.01 0.1 1 10 100 1000 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (b) Naive 0.001 0.01 0.1 1 10 100 1000 1e-08 1e-07 1e-06 1e-05 0.0001 0.001 L2(Eficiencia) L2(Error exacto) (c) Upper-bound Figura 5.14: Eficiencias de v Podemos destacar varias conclusiones: Cada uno de los modelos de estimador funciona pr´acticamente como un offset. La forma de las gr´aficas se mantienen, solo cambian los niveles. La forma de las graficas de pyvse parecen bastante entre s´ı, la de upresenta variaciones bastante significativas. En el caso de las Figuras 5.9, 5.10 y 5.11, las zonas con m´as variabilibad de la eficiencia son para los extremos inferior y superior del tama˜no de elemento. Es dif´ıcil obtener una relaci´on clara entre las variables. 36 Cap´ıtulo 6 Conclusiones y posibles trabajos futuros A tenor de los resultados obtenidos a lo largo del estudio y expuestos en el apartado anterior, podemos hacer las siguientes consideraciones con respecto a las cuestiones expuestas en los objetivos del proyecto. Seg´un se ha visto en el apartado 5.2.1, el VMSEE funciona de manera muy parecida para el flujo aerodin´amico alrededor de todas las geometr´ıas analizadas. Se ha observado un mejor rendimiento conforme va a aumentando el espesor del perfil, sobre todo en la velocidad horizontal u. No todas las variables se comportan de la misma manera frente al VMSEE: •El error de la variable presi´on suele estar sobrestimado. Este comportamiento es bueno ya que nos deja del lado de la seguridad. •En las dos variables de la velocidad, hay zonas donde el m´etodo subestima el error. Esto ocurre principalmente en la zona cercana al obst´aculo. Este comportamiento seguramente es debido a que el VMSEE no tiene en cuenta el error debido a las condiciones de contorno. Esta podr´ıa ser una futura l´ınea de trabajo. •A consecuencia de las dos cuestiones inmediatamente anteriores, el comportamiento de vy sobre todo use va a considerar cr´ıtico, por lo que se le prestar´a especial atenci´on. El m´etodo obtiene bajos rendimientos en las zonas muy alejadas del perfil, donde el error es extremadamente bajo. Esto podr´ıa ser debido a cuestiones de precisi´on num´erica. En cualquier caso, debido a que este es un m´etodo orientado al remallado, esta cuesti´on no tiene demasiada importancia. El error estimado seg´un la norma L2suele ser algo mayor que el de la norma L1. Por tanto se puede considerar la L2como la norma m´as recomendable. Segun se puede ver en el apartado 5.2.3 Las eficiencias del error son pr´acticamente independientes de la matriz de escalas temporales (τf) con las que se ha calculado el 37 ANEXO A. M´ ETODO DE POGGI De esta manera, la ecuaci´on A.5 queda: ∂2φ ∂2x+∂2φ ∂2y= ∂φ ∂x ∂q2 ∂x +∂φ ∂y ∂q2 ∂y 2c2 ∞h1−γ−1 2q2 q2 ∞−1q2 ∞ c2 ∞i(A.8) Sabemos adem´as que la velocidad del sonido del flujo imperturbado es M=q∞ c∞.∇2es el operador de Laplace. 1−γ−1 2M2q2 q2 ∞ −1∇2φ=1 2M2∂φ ∂x ∂ ∂x q2 q2 ∞+∂φ ∂y ∂ ∂y q2 q2 ∞ (A.9) Se desarrolla el campo φen una serie ascendente de potencias pares de M. φ=φ0+M2φ1+M4φ2+··· (A.10) Si definimos una velocidad compleja como: wi=∂φi ∂x −i∂φi ∂y (A.11) Se sigue que: w=wo+M2w1+M4w2+··· (A.12) Si adem´as: q2=ww= (u+iv)(u−iv) (A.13) Desarrollando la f´ormula anterior e igualando los t´erminos con la misma potencia de M, nos quedan las siguientes expresiones: ∇2φ0= 0 (A.14) ∇2φ1=1 2∂φ0 ∂x ∂ ∂x w0w0 q2 ∞+∂φ0 ∂y ∂ ∂y w0w0 q2 ∞ (A.15) . . . En este caso, φ0es el flujo potencial incompresible. Si nos fijamos en la ecuaci´on A.15, tiene un t´ermino fuente dependiente de φ0. El m´etodo de Poggi consiste en considerar el flujo compresible como una distribuci´on continua de fuentes en la regi´on externa a la frontera del s´olido. La intensidad de esta fuente vendr´a dada por: −1 4πq2 ∞∂φ0 ∂x ∂ ∂x w0w0 q2 ∞+∂φ0 ∂y ∂ ∂y w0w0 q2 ∞dxdy (A.16) Se introducen unas nuevas variables independientes, z=x+iy yz=x−iy. Por tanto x=1 2(z+z) e y=i 2(z−z). Utilizando la regla de la cadena queda:    ∂ ∂x =∂ ∂z +∂ ∂z ∂ ∂y =i∂ ∂z −∂ ∂z  (A.17) 44 ANEXO A. M´ ETODO DE POGGI As´ı tenemos: −1 2πq2 ∞∂φ0 ∂z ∂ ∂z w0w0 q2 ∞+∂φ0 ∂z ∂ ∂z w0w0 q2 ∞dxdy (A.18) Y la intensidad de la fuente quedar´a: −1 4πq2 ∞w2 0 dw0 dz +w2 0 dw0 dz dxdy (A.19) Conocida la intensidad de la fuente en el dominio fluido en Q(z), debemos calcular el efecto en Zpde su interacci´on con un cilindro de secci´on recta circular de radio R1con centro en el origen. Para ello, enunciamos el siguiente teorema, que podemos encontrar en Milne-Thomson [1962]. Teorema 1 (Teorema del c´ırculo) Sea un flujo irrotacional bidimensional no viscoso e incompresible en el plano complejo sin fronteras r´ıgidas. Sea el potencial complejo del flujo f(z), cuyas singularidades est´an a una distancia del origen mayor que a. Si se introduce en el flujo un cilindro circular de secci´on recta C,|z|=a, el potencial complejo se transforma en: w=f(z) + fa2 z(A.20) Una fuente de intensidad unidad situada en btiene como funci´on potencial: φs(z) = −ln(z−b) (A.21) Aplicando el Teorema 1 y reglas del algebra en variable compleja tenemos φt(z) = −ln(z−b) + ln a2 z−b= =−ln(z−b)−ln a2 z−b!= =−ln(z−b)−ln a2 z−b(A.22) A continuaci´on sumamos la constante ln −1/ba˜nadiendo las modificaciones pertinentes a la expresi´on anterior. φt=−ln(z−b)−ln a2 z−b=−ln(z−b)−ln 1−a2 bz+ ln −1 b(A.23) Procedemos de igual manera con ln(z) φt(z) = −ln(z−b)−ln z−a2 b+ ln(z) + ln −1 b(A.24) 45 ANEXO A. M´ ETODO DE POGGI Para calcular la velocidad compleja, procedemos a derivar la expresi´on anterior, teniendo en cuenta que el ´ultimo sumando es una constante. dφt dt =−1 z−b−1 z−a2 b +1 z=1 b−z+1 a2 b−z+1 z(A.25) Por comodidad, tendremos en cuenta los siguientes cambios de variable: z=zp b=z a=R1(A.26) Quedando por fin: dφt dt =1 z−zp +1 R2 1 z +1 zp (A.27) Figura A.1: Imagen de una fuente unidad con respecto a un c´ırculo 46 ANEXO A. M´ ETODO DE POGGI Sabemos que R(R2 1 z) es el punto inverso de zcon respecto a una circunferencia de centro en el origen y radio R1, por lo que podemos enunciar: El potencial creado por un t´ermino fuente unidad Q(z) en presencia de una circunferencia centrada en el origen, es equivalente al creado por dos fuentes, una en zy otra en el punto inverso de zcon respecto a dicha circunferencia y una tercera fuente unidad en el origen de coordenadas. Agrupando t´erminos por comodidad, tenemos: dφt dt =1 z−zp +R2 1/zp R2 1−zpz(A.28) En el caso que nos ocupa, en cada punto externo al cilindro hay un t´ermino fuente cuaya intensidad viene dada por la f´ormula A.19 y cuyo efecto sobre un punto Pes como se ha expresado en A.28. Para calcular el efecto de todas estas fuentes ser´a necesario integrar en el espacio cada una de sus contribuciones. Figura A.2: Integral de los t´erminos fuente: teorema de Stokes w1(P) = −1 4πq2 ∞RR d dz 1 z−zpw2 0w0dxdy −1 4πq2 ∞RR d dz R2 1/zp R2 1−zpzw2 0w0dxdy −1 4πq2 ∞RR 1 z−zp d dz (w2 0w0)dxdy −1 4πq2 ∞RR R2 1/zp R2 1−zpz d dz (w2 0w0)dxdy (A.29) Esta ecuaci´on resulta bastante complicada, pero mediante teorema de Stokes se puede transformar en una ecuaci´on m´as manejable. Hay que destacar que hasta el momento, la ´unica simplificaci´on que se ha propuesto es despreciar los terminos de M4y superiores. Teorema 2 (Teorema de Stokes en el plano complejo) Si f(z, z)es una funci´on de z=x+iy yz=x−iy, que es continua y diferenciable en un ´area Sencerrada por un 47 ANEXO A. M´ ETODO DE POGGI contorno C, entonces: Ic f(z, z)dz = 2iZZS ∂f ∂z dS (A.30) Ic f(z, z)dz =−2iZZS ∂f ∂z dS (A.31) Consideremos un contorno en el plano Z consistente en una curva cerrada C1que contiene el obst´aculo,una peque˜na curva C2que contiene el punto zp, una gran curva C3 que contiene a las dos anteriores y sendas rectas AB yCD que conectan c2yc1con c3 (Figura A). El contorno as´ı descrito se recorre de manera que la superficie squeda a la izquierda. Ahora supongamos f(z, z) una funci´on compleja definida sobre este dominio. Entonces, por el Teorema 2, para transformar las integrales de superficie en integrales de l´ınea, es necesario expresar los integrandos como derivadas con respecto a zoz. Los dos primeros integrandos de la ecuaci´on A.29 ya est´an en esta forma. Ser´a necesario manipular los dos siguientes. w2 0=d dz Zw2 0dz (A.32) w2 0=d dz Zw2 0dz (A.33) (A.34) Por medio del Teorema 2 y de las dos ecuaciones anteriores, tenemos que la Ecuaci´on A.29 queda w1(p) = −1 8πiq2 ∞Z1 z−zp F(z, z)dz −R2 1/Z2 p 8πiq2 ∞Z1 Z−R2 1/zp F(z, z)dz (A.35) donde F(z, z) = w2 0w0+d dz (w0)Zw02dz (A.36) F(z, z) = w02w0+d dz (w0)Zw2 0dz (A.37) Se toman las integrales de la Ecuaci´on A.35 sucesivamente alrededor de cada uno de los c´ırculos de la Figura . Cada una de ellas se toma en sentido antihorario, por lo que se tendr´an en cuenta los signos del Teorema de Stokes, negativo para C1yC2y positivo para C3. Para hayar el resultado de las integrales con comodidad se usar´a el Teorema del Residuo de Cauchy. Teorema 3 (Teorema del residuo) Si f(z)es anal´ıtica en un contorno Cy en su interior excepto en un n´umero de polos dentro del contorno, entonces: Ic f(z)dz = 2πiM (A.38) 48 ANEXO A. M´ ETODO DE POGGI donde Mdenota la suma de los residuos de f(z)interiores a cde igual forma: Ig(z)dz =−2πiN (A.39) donde Ndenota la suma de los residuos de g(z)interiores a c. Se sabe que el valor del residuo de un polo simple de multiplicidad nsituado en ade la funci´on f(z) es: Res(f, a) = 1 (n−1)! dn−1 dzn−1l´ım z→a[(z−a))nf(z))] (A.40) Consideremos la integral de la ecuaci´on A.35 . En ´el z=R2/z y por tanto F(z, z) = F(z, R2/z). De esta manera ya tenemos una funci´on anal´ıtica depediente de zsobre la que aplicar la Ecuaci´on (A.38). Con la segunda integral se proceder´ıa de igual forma, teniendo en cuenta que F(z, z) = F(z, R2/z), lo que nos habilita para aplicar la Ecuaci´on (A.38). Adem´as, dado que F(z, z) es el conjugado de F(z, R2/z), la segunda integral se puede escribir en t´erminos de la primera. Consideremos la primera integral de l´ınea de (A.35) en torno al contorno C1: Los polos asociados con terminos que solo dependen de zest´an dentro del contorno, mientras que los asociados exclusivamente a zest´an fuera, igual que el polo en zp. De esta manera, solo contribuyen a la integral los polos de z. Por tanto, seg´un la Ecuaci´on (A.38): ZC1 1 z−zp F(z, z)dz = 2πiS(zp) (A.41) donde S(zpdenota la suma de los residuos de 1 z−zpF(z, z) dentro del contorno C1. Integrando en el c´ırculo de radio infinitesimal centrado en p,C2; se observa que s´olo hay un polo simple. Como en el l´ımite R2→0, z→zpyz→zp, se sigue que: ZC2 1 z−zp F(z, z)dz = 2πiF(zp, zp) (A.42) Finalmente, tenemos la integral alrededor del c´ırculo centrado en el origen con R3→ ∞. Para conseguir que el integrando sea anal´ıtico se reemplaza zpor su conjugado R2 3/z. Para realizar la integral ser´a necesario desarrollar la integral cuando R3→ ∞. En la lejan´ıa de un obst´aculo circular, el campo de velocidades es el correspondiente al del flujo imperturbado: w∞ 0=q∞(eiα) (A.43) teniendo en cuenta que α= 0 w∞ 0=q∞(A.44) obtenemos: ZC3 1 z−zp F(z, z)dz = 2πiq3 ∞(A.45) 49 ANEXO A. M´ ETODO DE POGGI A continuaci´on se procede de igual manera con la segunda integral de A.35. Integrando alrededor de C1tenemos Z1 Z−R2 1/zp F(z, z)dz =−2πiS(R2 1/Zp)−2πiF (R2 1/zp, zp) (A.46) La integral C2no contribuye, dado que no hay polos de Zen el interior de su contorno. El segundo termino de la ecuaci´on anterior es an´alogo al calculado en (A.42). Finalmente, sabemos que el resultado de la integral alrededor de C3es el conjugado de lo obtenido en la Ecuaci´on (A.45), por lo que Z1 Z−R2 1/zp F(z, z)dz =−2πiq3 ∞(A.47) En resumen, tenemos que: w1=1 4q2 ∞S(z)−R2 1 z2S(R2 1 zi)−R2 1 z2F(R2 1 z, z)+1 4q2 ∞ F(z, z)−q∞ 41−R2 1 z2(A.48) As´ı, recordando la Ecuaci´on (A.12),la velocidad alrededor del c´ırculo ser´ıa: w=w0+M2w1 Ahora, sustituyendo la expresi´on de la velocidad w0=q∞(1−R2 1 Z2) en la Ecuaci´on (A.48) y aplicando (A.40) y abandonando el sub´ındice p: w1 q∞ =−1 21−R2 1 z2 1+R2 1 2 (1−R2 1)(z2+1) z2−1+R2 1 4 (z2−R2 1)3 z2(z2−1)(z2−R4 1) −(R2 1−1)(1−R2 1) 8hz2+1 (z2−1)2(R2 1−1) + R2 1 z2+R4 1 (z2−R4 1)2(1 −R2 1)ilog R2 1+1 R2 1−1 +R4 1 4 z(1−R2 1)2(R2 1−1) (z2−R4 1)2log z+1 z−1 −1 4n1−R2 1 Z2R2 1 z2−1z2z2 (z2−1)(z2−1) + 2 (1 −R2 1)z (z2−1)2hz+R4 1 z+R2 1−1)2 2log z−1 z+1io (A.49) 50 Anexo B Adimensionalizaci´on del problema En este anexo se pretende establecer con cuantos par´ametros adimensionales se caracteriza el flujo sobre el que se est´a realizando el estudio. En este caso se trata de flujo de Euler compresible alrededor de un perfil aerodin´amico. Para realizar el an´alisis dimensional es necesario identificar las variables relevantes del problema. Esto no ser´a problem´atico, ya que conocemos las ecuaciones que lo rigen. Podemos escribir las ecuaciones de conservaci´on de masa, momento y energ´ıa para un fluido compresible no viscoso como: ∂ρ ∂t +∇ · (ρu u u) = 0 (B.1) ∂ρu u u ∂t +∇ · (ρu u uu u u) + ∇p= 0 (B.2) ∂ρetot ∂t +∇ · [(ρetot +p)u u u] = 0 (B.3) donde ρes la densidad del fluido, u u ues la velocidad, pes la presi´on y etot la energ´ıa total por unidad de masa. Esta ´ultima se puede descomponer en dos contribuciones: energ´ıa interna y cin´etica. etot =e+u u u·u u u/2. La ecuaci´on (B.3) se puede simplificar a´un m´as teniendo en cuenta que: ∇ · [(ρetot +p)u u u] = ∇ · ρu u uetot +u u u·u u u/2 + p ρ=∇ · ρu u uh +ρu u uu u u·u u u 2(B.4) Si asumimos las ecuaciones para gas ideal del fluido tenemos que p=ρ¯ RT (B.5) y adem´as h= ¯cpT(B.6) con ¯ Rconstante del gas y ¯cpcalor espec´ıco del gas, ambos en base m´asica. Por tanto nos quedan cinco ecuaciones y cinco inc´ognitas. Sabiendo esto, tenemos que por ejemplo que el campo de velocidades u u uviene determinado por U∞,T∞,p∞, y dos par´ametros 51 ANEXO B. ADIMENSIONALIZACI´ ON DEL PROBLEMA M L T Θ u u u0 1 −1 0 a0 1 0 0 T∞0 0 0 −1 U∞0 1 −1 0 p∞−1 1 −2 0 ¯ R0 2 −2−1 ¯cp0 2 −2−1 dependientes del gas, ¯ Ry ¯cp; as´ı como a, la cuerda del perfil. Por el teorema Pi, sabemos que habr´a 7 −4 = 3 par´ametros adimensionales Por tanto sabemos que: Π0=u u u U∞ (B.7) Π1=¯ R ¯cp =γ−1 γ(B.8) Π2=u u u2 γ¯ RT∞ =M2(B.9) Por lo tanto tenemos que: u u u U∞ =f(γ, M) (B.10) siendo Mel n´umero de Mach del flujo imperturbado y γla relaci´on de calores espec´ıficos. Procediendo an´alogamente para el campo de presiones p p p∞ =f(γ, M) (B.11) y para el de temperaturas TT T∞ =f(γ, M) (B.12) Como en este proyecto, el fluido es en todo momento aire, γqueda fijo, por lo que el flujo solo depende de M. 52 Anexo C Figuras suplementarias En este anexo se presentan los mapas de color de la zona del perfil relacionados con lo que en el Cap´ıtulo 5 se ha denominado caso base. Se van a examinar las soluciones de elementos finitos, la soluci´on exacta y las eficiencias seg´un los tres modelos de estimador para la norma L2. Destacan los siguientes resultados: Es muy dif´ıcil distinguir a simple vista la soluci´on exacta de la soluci´on calculada por elementos finitos. Mientras que en la presi´on y la velocidad vertical se ve que la zona con mayor error es el borde de ataque del perfil, en el caso de la velocidad horizontal no est´a tan claro. A lo largo de todo el perfil se observa que la velocidad obtenida por elementos finitos es alrededor de un 2 % m´as baja que la soluci´on exacta. Como se ha comentado anteriormente, todos los modelos de estimador y todas las variables presentan en mayor o menor medida una zona de eficiencia baja en las proximidades del perfil. 53