Desarrollo e implementación en Matlab de un código numérico para la resolución de problemas de Aerodinámica potencial linealizada no estacionaria

i Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 Desarrollo e implementación en Matlab de un código numérico para la resolución de problemas de Aerodinámica potencial linealizada no estacionaria Autor: Fernando Moreno Pino Tutor: José Manuel Gordillo Arias de Saavedra Trabajo Fin de Master Master Universitario en Ingeniería Aeronáutica ii iii Trabajo Fin de Master Master Universitario en Ingeniería Aeronáutica Desarrollo e implementación en Matlab de un código numérico para la resolución de problemas de Aerodinámica potencial linealizada no estacionaria Autor: Fernando Moreno Pino Tutor: José Manuel Gordillo Arias de Saavedra Dep. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2019 iv v Trabajo Fin de Master: Desarrollo e implementación en Matlab de un código numérico para la resolución de problemas de Aerodinámica potencial linealizada no estacionaria Autor: Fernando Moreno Pino Tutor: José Manuel Gordillo Arias de Saavedra El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2019 El Secretario del Tribunal vi vii A mi familia A mis maestros A mis amigos viii ix Agradecimientos En primer lugar quiero dar las gracias a mi padre, a mi madre y a mi hermano, pues han sido la base en la que apoyarme desde mucho antes que estudiará el Grado en Ingeniería Aeroespacial y por supuesto durante el Master en Ingeniería Aeronáutica. Dar gracias también a mi compañera de viaje, Bea, que desde que la conozco ha estado ahí para apoyarme, calmarme, y sacarme una sonrisa cada vez que lo necesitaba y sin esfuerzo aparente. Muchas gracias por todo Quiero agradecer también a todas las personas que me han demostrado que me quieren y sé que puedo contar con ellas para cualquier cosa: Guille, Violeta, Dani, Elena, Álvaro, Fran, Alex, Pablo… la lista podría contar con una gran cantidad de nombres de personas cuya forma de ser consigue que siempre piense que todo va a salir bien. Gracias por hacerme reír y por haber compartido conmigo tantos momentos. Es un honor formar parte del grupo “Aeromafia” con todos ustedes. Dar las gracias a las personas que me esperaban en el pueblo cada vez que he vuelto: mis abuelas Gonzala y Antonia, mi abuelo “Juani”, mi abuelo Paco, que en paz descanse, y al resto de mi familia. A mis amigos, los “empresarios” porque aunque estén locos y sean unos payasos, siempre serán mis payasos. Finalmente, agradecer a mi tutor y acompañante en este proyecto, José Manuel Gordillo, pues con él empezó mi andadura en el mundo de la aerodinámica y con él cierro mi etapa universitaria, con un proyecto que trata dicho campo. Gracias por haber estado ahí todos y cada uno de los días que he tenido dudas y me he acercado a consultarle y por el trato que ha tenido conmigo. Ha sido un placer. xvi 10. Montaje............................................................................................................................... 110 xvii ÍNDICE DE FIGURAS Figura 1: Geometría Típica de un Ala en Flecha ..................................................................................... 23 Figura 2: Superficies de Integración....................................................................................................... 37 Figura 3: Normales en el entorno de la Singularidad ............................................................................. 38 Figura 4: Discretización de la geometría de Ala y Estela ......................................................................... 45 Figura 5: Mallado triangular simétrico en Ala y Estela ........................................................................... 45 Figura 6: Cambio de variable para la integración. Gráfico explicativo I ................................................... 47 Figura 7: Ángulos empleados en la Integración...................................................................................... 47 Figura 8: Ángulos empleados en la discretización. Gráfico explicativo II................................................. 51 Figura 9: Definición de 𝑟(𝜃) .................................................................................................................. 51 Figura 10: Mallado en el entorno de la singularidad .............................................................................. 52 Figura 11: Potencial en el borde de salida y en la estela en régimen estacionario para un valor de y dado ............................................................................................................................................................. 56 Figura 12: Mallado triangular en régimen estacionario .......................................................................... 56 Figura 13: Condiciones de contorno para el Potencial en régimen estacionario ..................................... 57 Figura 14: Notación Matricial Usada...................................................................................................... 58 Figura 15: Condición de Kutta para el caso no estacionario. .................................................................. 61 Figura 16: Propagación de los potenciales en la Estela no estacionaria .................................................. 63 Figura 17: Discretización en alas rectangulares ..................................................................................... 65 Figura 18: Potencial sobre el extradós de alas rectangulares ................................................................. 66 Figura 19: Velocidad perturbada en el eje x sobre el extradós de alas rectangulares. ............................ 66 Figura 20: Velocidad perturbada en el eje x sobre el extradós de alas rectangulares. Refinamiento ...... 67 Figura 21: Discretización en alas con flecha .......................................................................................... 68 Figura 22: Potencial sobre el extradós de alas con flecha ...................................................................... 69 Figura 23: Velocidad perturbada en el eje x sobre el extradós de alas con Flecha. ................................. 69 Figura 24: Evolución de la Pendiente de la Curva de Sustentación con el tiempo. Ala Rectangular ......... 71 figura 25: Evolución de la pendiente de la curva de sustentación con el tiempo. ala con flecha ............. 72 Figura 26: Estela no Estacionaria en caso de Ala en Flecha .................................................................... 73 Figura 27; Evolución de la Pendiente de la Curva de Sustentación con el Alargamiento. Ala Rectangular 77 Figura 28: Evolución de la Pendiente de la Curva de Sustentación con el Alargamiento. Ψ=30º ......... 77 Figura 29: Evolución de la Pendiente de la Curva de Sustentación con el Alargamiento.Ψ=45º .......... 78 Figura 30: Movimiento del Fluido sobre el extradós de aeronaves con flecha ........................................ 79 Figura 31: Velocidad perturbada en el eje x sobre el extradós de alas con Flecha Negativa ................... 79 Figura 32: Respuesta ante un escalón ................................................................................................... 81 Figura 33: Evolución del Coeficiente de Sustentación en un movimiento oscilatorio .............................. 82 Figura 34: Evolución del coeficiente de sustentación en un movimiento oscilatorio. (KyP) [2] ............... 82 xviii xix ÍNDICE DE TABLAS Tabla 1; Resultados para Ala Rectangular en Régimen Estacionario ......................................... 67 Tabla 2: Resultados para Ala Rectangular en Régimen Estacionario. Mallado Refinado ........... 67 Tabla 3: Resultados para ala con Flecha 30º en régimen estacionario. .................................... 69 Tabla 4: Resultados para Ala Rectangular en Régimen No Estacionario ................................... 71 Tabla 5: Resultados para Ala con Flecha 45º en Régimen Estacionario .................................... 72 Tabla 6: Resultados obtenidos por distintos métodos para alas rectangulares con distinto alargamiento ............................................................................................................................................... 75 Tabla 7: Resultados obtenidos por distintos métodos para alas con Ψ= 30º con distinto alargamiento ............................................................................................................................................... 75 Tabla 8: Resultados obtenidos por distintos métodos para alas con Ψ= 45º con distinto alargamiento ............................................................................................................................................... 75 Tabla 9: Error Cometido en Alas Rectangulares en Régimen Estacionario ................................ 80 Tabla 10: Error Cometido en Alas con Ψ=30º en Régimen Estacionario ................................ 80 Tabla 10: Error Cometido en Alas con Ψ=45º en Régimen Estacionario ................................ 80 xx 21 1. INTRODUCCIÓN Desde la invención de la aviación y hasta nuestros días, la tecnología empleada en este campo ha evolucionado considerablemente. Si se comparan los aviones de hoy en día con sus ancestros de la segunda guerra mundial o incluso con los anteriores biplanos, se observa que el avance es abrumador. Los campos en los que más se nota esta diferencia a simple vista son la aviónica que llevan a bordo las aeronaves actuales, la capacidad que son capaces de transportar, la propulsión, etc… Sin embargo y aunque a simple vista no se observe, a todas las aeronaves se les pide 2 cosas fundamentalmente: que vuelen y que aguanten las cargas a las que se ven sometidas. La disciplina que se encarga de la segunda petición que se les hace a las aeronaves es la elasticidad y la resistencia de materiales, la cual estudia el comportamiento de la estructura de la aeronave frente a la acción de las distintas fuerzas a las que se ve sometida durante el vuelo y establece las cargas límite que esta es capaz de aguantar sin poner en riesgo a la seguridad de los pasajeros, tripulación y estructura. Por otro lado, la disciplina que se encarga de que una aeronave vuele es la mecánica de fluidos y más específicamente la aerodinámica. Esta disciplina se encarga de estudiar cómo es el movimiento del fluido alrededor de la aeronave. Esto permite modelar cómo es el campo de presiones sobre la aeronave y por tanto conocer las fuerzas que tiene que soportar la estructura. Estas dos disciplinas son la base para lo que fundamentalmente se les pide a las aeronaves. Sin embargo actualmente se ha conseguido ver de forma tan natural el hecho de que una aeronave despega y aterriza, que la preocupación de las compañías aéreas recae en incorporar y desarrollar sistemas que consigan una mayor eficiencia del vuelo. En consecuencia, la aeroelasticidad, que es la disciplina que estudia la interacción tanto de la parte elástica como de la aerodinámica, es una disciplina fundamental en el diseño de cualquier tipo de aeronaves. Aquí entran en juego dos conceptos: eficacia y eficiencia. La eficacia se define cómo la capacidad de lograr el efecto que se desea, mientras que la eficiencia implica hacerlo con el mínimo desaprovechamiento de los recursos disponibles. Una vez que la aviación ya se consiguió que fuese eficaz, ahora se desarrollan conceptos y se investiga para mejorar su eficiencia. Si se analiza la eficiencia de un vuelo de un punto a otro, se puede observar que el vuelo más eficiente es el que se desarrolla en régimen de crucero y a ser posible a un valor de coeficiente de sustentación constante (cruise climb). 22 Este trabajo tiene como objetivo desarrollar un método de cálculo que permita obtener valores reales del coeficiente de sustentación de alas para distinta geometría y compararlos con los resultados experimentales de la literatura de forma que se puedan extraer conclusiones acerca de las ventajas e inconvenientes de las diferentes geometrías. Además, en este trabajo no sólo se va a estudiar el movimiento en régimen estacionario, sino también el movimiento de la aeronave en régimen no estacionario. Para ello en los distintos puntos del trabajo se van a tratar los siguientes temas: - En primer lugar se va a llevar a cabo la deducción de las ecuaciones del problema aerodinámico y las simplificaciones e hipótesis necesarias para plantear la ecuación a resolver con el método que se va a plantear. - A continuación se presentará la discretización del problema y resolución numérica. Se explicarán los aspectos más importantes a tener en cuenta para resolver los problemas numéricamente. - Tras esto, se exponen los resultados obtenidos para alas de distinta geometría tanto para el caso estacionario como para el caso no estacionario y se compararán con algunos obtenidos de la literatura así como otros obtenidos haciendo uso de otro software comercial. - Finalmente se expondrán las conclusiones obtenidas del trabajo y se propondrán proyectos de investigación futuros que permitan mejorarlo y complementarlo. 23 2. GEOMETRÍA En este punto se van a analizar las características de la geometría sobre la que se va a realizar el estudio. Dicha geometría consiste principalmente en un ala con forma en planta rectangular, pero posteriormente se realiza la extensión a alas con geometría trapezoidal genérica. Se considerará que el ala tiene su forma en planta contenida en el plano Z = 0. Las características principales que definen a este tipo de alas son: - La envergadura (b): Definida como la distancia entre los bordes marginales del ala. - La cuerda en la raíz (cr): Definida como la distancia entre el borde de ataque y el borde de salida del ala medida en la raíz de la misma. En nuestro caso es la cuerda del ala en el plano Y=0. - Superficie alar (S): Superficie de la forma en planta del ala, definida como la proyección del ala sobre el plano Z=0. - ct: es la cuerda en las puntas. - 𝛹: es el ángulo de flecha del ala medido desde el borde de ataque. - α: es el ángulo de ataque. Se recuerda que para alas con flecha, la geometría viene dada por: 𝑏=√AR·𝑆 ; 𝑐𝑟=2 √𝑆 /AR 1+𝐸 ; 𝑐𝑡=𝐸𝑐𝑟 FIGURA 1: GEOMETRÍA TÍPICA DE UN ALA EN FLECHA 24 𝑥𝑎(𝑦)=|𝑦|∙tan(𝛹) 𝑥𝑠(𝑦)=𝑐𝑟+𝑏2tan(Ψ)+𝑐𝑡−𝑐𝑟 𝑏2|𝑦| Donde se definen los parámetros alargamiento y estrechamiento como sigue: AR=𝑏2 𝑆 𝐸=𝑐𝑟 𝑐𝑡 Y se ha utilizado la siguiente nomenclatura: - xa: define los puntos x del borde de ataque del ala. - xs: define los puntos x del borde de salida del ala. La Forma en planta del ala se puede describir de forma general con la siguiente función: 𝐹(𝑥,𝑦,𝑧,𝑡)≡𝑍w(𝑥,𝑦,𝑡)−𝑧=0 El ala estará compuesta por perfiles aerodinámicos esbeltos, lo cual será necesario para la posterior simplificación y linealización de las ecuaciones generales de la mecánica de fluidos, con el fin de obtener unas ecuaciones más sencillas y resolubles. La condición de esbeltez de los perfiles aerodinámicos implica que las variaciones de la geometría a lo largo de las direcciones longitudinal y transversal del ala son pequeñas. Matemáticamente esto es: 𝜕𝑍𝑤 𝜕𝑥~∆𝑍𝑐 ∆𝑥𝑐~𝑧0 𝑐𝑟≪1 𝜕𝑍𝑤 𝜕𝑦~∆𝑍𝑐 ∆𝑦𝑐~𝑧0 𝑏≪1 Siendo z0 un valor de espesor característico del ala. De forma general, la ecuación que define la geometría del ala se puede expresar como la suma de 2 contribuciones. 𝑍𝑤(𝑥,𝑦,𝑡)=𝑍𝑐(𝑥,𝑦,𝑡)±𝑍𝑒(𝑥,𝑦) Donde Zc(x,y,t) representa los puntos de las líneas de curvatura de cada perfil en cada instante de tiempo y Ze(x,y) es la función que representa el espesor de cada perfil. Que el valor de la función de espesor de cada perfil no dependa del tiempo, es indicativo de que el perfil no sufrirá ninguna deformación a lo largo del problema. 25 3. ECUACIONES Y CONDICIONES DE CONTORNO 3.1. Ecuaciones Generales En este punto se plantearán las ecuaciones que gobiernan el problema a partir de las ecuaciones generales de la Mecánica de Fluidos y se impondrán las condiciones de contorno del problema. El movimiento del aire alrededor del ala viene regido por las ecuaciones de Navier-Stokes, las cuales se pueden escribir de forma general como se muestra a continuación.  Ecuación de Continuidad: 𝜕𝜌 𝜕𝑡+∇·(𝜌𝑣)=0 (3.1.1)  Ecuación de Cantidad de Movimiento: 𝜌𝜕𝑣 𝜕𝑡+𝜌𝑣·∇𝑣=−∇𝑝+∇·𝜏′+𝜌𝑓𝑚 󰇍 󰇍 󰇍 󰇍  (3.1.2)  Ecuación de la Energía: 𝜌𝑐𝑣𝜕𝑇 𝜕𝑡+𝜌𝑐𝑣𝑣·∇𝑇=−𝑝∇·𝑣+𝜏′:∇𝑣+𝑄𝑟+𝑄𝑞+∇·(𝑘∇𝑇) (3.1.3) Este sistema conforma un conjunto de 5 ecuaciones con 6 incógnitas: - 𝜌: Densidad del fluido - 𝑣: Velocidad del fluido - P: Presión del fluido - T: Temperatura del fluido Para cerrar el problema se debe completar el sistema con la ecuación de estado del fluido bajo estudio, el cual se va a considerar que es un gas perfecto.  Ecuación de Estado: 𝑝 𝜌=𝑅𝑔𝑇 (3.1.4) Donde 𝑅𝑔=𝑅 𝑀𝑚=𝑐𝑝−𝑐𝑣=𝑐𝑣(𝛾−1). Siendo 𝑅=8,314 𝐽 𝑚𝑜𝑙·𝐾 la constante universal de los gases y 𝑀𝑚 la masa molar del gas. Con esto, se tiene un sistema no lineal de ecuaciones en derivadas parciales con 6 ecuaciones y 6 incógnitas al que habrá que proporcionarle las condiciones de contorno e iniciales adecuadas para cada problema. 32 ℎ=𝐶𝑝𝑇=𝐶𝑝𝑝 𝑅𝑔𝜌=(𝐶𝑝 𝐶𝑣)𝑝 (𝑅𝑔 𝐶𝑣)𝜌=𝛾 𝛾−1𝑝 𝜌=𝑎2 𝛾−1 Se tiene que: 𝜕𝜙 𝜕𝑡+|∇𝜙|2 2+𝑎2 𝛾−1=𝐶(𝑡) En este proyecto se considerará únicamente el límite 𝑀∞≪1 lo que permite simplificar las ecuaciones aún más. En este límite 1 𝜌Δ𝜌<< 1. Lo que implica que el fluido puede ser considerado como incompresible. En este caso la ecuación de la continuidad se simplifica a: ∇·𝑣′ 󰇍 󰇍 󰇍  =0 Los términos siguiente pueden ser despreciados. Véase Gordillo y Riboux, 2012 [1] ∂ρ ∂t ; 𝑣′ 󰇍 󰇍 󰇍  ·∇𝜌≪𝜌∇·𝑣′ 󰇍 󰇍 󰇍  Para el caso 𝜌≈𝑐𝑡𝑒, la ecuación escalar para el potencial queda: 𝜕𝜙 𝜕𝑡+|∇𝜙|2 2+𝑝 𝜌=𝐶(𝑡) (3.4.2) Se observa que la ecuación de cantidad de movimiento ha pasado a ser una ecuación escalar para el potencial de velocidades. 𝐶(𝑡) es una constante que es únicamente una función del tiempo. Conocido el potencial de velocidades, se puede utilizar la ecuación anterior para determinar la distribución de presiones. En el caso de 𝜌≈𝑐𝑡𝑒, la sustitución de la velocidad en función del potencial en la ecuación de continuidad proporciona que el potencial de velocidades satisface la ecuación de Laplace. ∇2𝜙=0 Con la hipótesis de flujo incompresible se estarían cometiendo errores de 𝑂(𝑀2), donde 𝑀=𝑣/𝑎 denota el número de Mach. Dicha demostración se puede encontrar en [1]. Esta hipótesis es válida para alas volando a un Mach<0.3. Debido a que la ecuación de Laplace es lineal, se puede aplicar el método de superposición de soluciones elementales para hallar la solución del problema general. El problema general queda definido por: ∇2𝜙=0 |𝑥|→∞, 𝜙→𝑈∞𝑥 (3.4.3) 33 |𝑥|∈Σ𝑠, 𝑛 󰇍  ·(∇𝜙−𝑣𝑝 󰇍 󰇍 󰇍 󰇍  )=0 𝐶𝑜𝑛𝑑𝑖𝑐𝑖𝑜𝑛 𝑑𝑒 𝐾𝑢𝑡𝑡𝑎−𝐽𝑜𝑢𝑘𝑜𝑤𝑠𝑘𝑖 Donde 𝑣𝑝 󰇍 󰇍 󰇍 󰇍  es la velocidad vertical del perfil en caso de encontarnos en régimen no estacionario y la condición de Kutta-Joukowski es la que permite determinar la solución real a nuestro problema de entre las infinitas soluciones existentes que tiene el problema del flujo potencial alrededor de un objeto. 3.5. Ecuaciones Linealizadas A continuación se va a proceder a linealizar las ecuaciones a resolver. Se va a considerar que el ala está formada por perfiles muy esbeltos y se mueve a ángulos de ataque pequeños. Esto permite suponer que la capa límite no se va a desprender en ningún punto del ala y en consecuencia, dicho ala perturba poco la corriente incidente. Esto quiere decir, que las perturbaciones en las variables del fluido, originadas por la presencia del ala, son de un orden de magnitud inferior a su valor aguas arriba. Las hipótesis a considerar son: 𝛼≪1 ℎ0≪𝑐 𝜕𝑧𝑝 𝜕𝑥≪1, 𝜕𝑧𝑝 𝜕𝑦≪1 (3.5.1) Donde ℎ0 es el espesor característico del perfil. Estas hipótesis junto con la de número de Reynolds muy elevado, nos permiten asegurar que la capa límite no se desprende y que la perturbación en las variables fluidas es pequeña. Con estas hipótesis los campos de velocidades, presión, densidad y temperatura de las variables fluidas se puede reescribir como la suma de su valor aguas arriba más una perturbación originada por la presencia del ala (3.5.2): 𝑣(𝑥)=𝑈∞ 󰇍 󰇍 󰇍 󰇍 󰇍  +𝑣′(𝑥), |𝑣′(𝑥)|≪𝑈∞ 𝑝(𝑥)=𝑝∞+𝑝′(𝑥), 𝑝′(𝑥)≪𝑝∞ (3.5.2) Donde 𝑣′(𝑥), 𝑝′(𝑥), 𝜌′(𝑥), 𝑇′(𝑥) son los campos perturbados de velocidad, presión, densidad y temperatura. El potencial de velocidades se puede descomponer igualmente como la suma del potencial de velocidades en el infinito más el potencial de velocidades perturbadas (3.5.3), que será la incógnita de nuestro problema. 𝜙(𝑥)=𝜙∞+𝜙′(𝑥) (3.5.3) 34 Debido a la linealidad de la ecuación de Laplace, el potencial de velocidades perturbadas debe de cumplir (3.5.4): ∇2𝜙′=0 (3.5.4) La ecuación de cantidad de movimiento en términos del potencial de velocidades perturbadas debe de cumplir (3.5.5) 𝜌𝜕𝜙′ 𝜕𝑡+𝜌𝑈∞𝜕𝜙′ 𝜕𝑥+𝑝′=0 (3.5.5) La condición de impenetrabilidad puede simplificarse gracias a la linealización de las ecuaciones: 𝑛 󰇍  ·(∇𝜙−𝑣𝑝 󰇍 󰇍 󰇍 󰇍  )=∇𝐹𝑒,𝑖 |∇𝐹𝑒,𝑖|·(∇𝜙−𝑣𝑝 󰇍 󰇍 󰇍 󰇍  )=0→∇𝐹𝑒,𝑖·(∇𝜙−𝑣𝑝 󰇍 󰇍 󰇍 󰇍  ) Donde 𝐹𝑒,𝑖 es la función de los puntos del extradós y del intradós del ala. Se puede reescribir pues la ecuación de la siguiente forma: (𝑒3 󰇍 󰇍 󰇍  −𝜕𝑧𝑝 𝜕𝑥𝑒1 󰇍 󰇍 󰇍  −𝜕𝑧𝑝 𝜕𝑦𝑒2 󰇍 󰇍 󰇍  )·((𝜕𝜙′ 𝜕𝑧−𝜕𝑧𝑝 𝜕𝑡)𝑒3 󰇍 󰇍 󰇍  +(𝑈∞+𝜕𝜙′ 𝜕𝑥)𝑒1 󰇍 󰇍 󰇍  +𝜕𝜙′ 𝜕𝑦𝑒2 󰇍 󰇍 󰇍  )=0 Despreciando los términos de segundo orden, se llega a (3.5.6): 𝜕𝜙′ 𝜕𝑧(𝑥,𝑦,𝑧=0)=𝑤′(𝑥,𝑦,𝑧=0)=𝑈∞𝜕𝑧𝑝 𝜕𝑥(𝑥,𝑦,𝑡)+𝜕𝑧𝑝 𝜕𝑡(𝑥,𝑦,𝑡) (3.5.6) Esta condición se debería de imponer estrictamente sobre el ala, pero debido a que los perfiles son esbeltos y los ángulos que vamos a manejar son pequeños, dicha condición se impone sobre la forma en planta del ala (FP). Gracias a estas hipótesis, esta aproximación no supone un error importante. Véase Gordillo y Riboux, 2012 [1]. El problema linealizado viene pues descrito por las siguientes ecuaciones: ∇2𝜙′=0 |𝑥|→∞, 𝜙′→0 |𝑥|∈𝐹𝑃,𝜕𝜙′ 𝜕𝑧(𝑥,𝑦,𝑡)=𝑈∞𝜕𝑧𝑝 𝜕𝑥(𝑥,𝑦,𝑡)+𝜕𝑧𝑝 𝜕𝑡(𝑥,𝑦,𝑡) 𝐶𝑜𝑛𝑑𝑖𝑐𝑖𝑜𝑛 𝑑𝑒 𝐾𝑢𝑡𝑡𝑎−𝐽𝑜𝑢𝑘𝑜𝑤𝑠𝑘𝑖 (3.5.7) 3.6. Condición de Kutta-Joukowky En esta apartado se explica la condición de Kutta-Joukowsky y sus implicaciones en el problema aerodinámico. 35 Dicha condición establece que de las infinitas soluciones que tiene el movimiento de un fluido alrededor del ala, es aquella para la cual el flujo no rebordea el borde de salida. Esta condición fija el valor de la circulación (Γ), permitiendo elegir de entre todas las soluciones aquella que se da en la realidad. La condición de contorno que se debe de imponer en el borde de salida se obtiene de esta condición, junto con la condición de que la diferencia de presiones entre estrados e intradós debe de ser nula a lo largo del borde de salida porque no existe sólido que soportase dicha diferencia en el caso de que la hubiese. En el extradós y en el intradós del borde de salida se debe de cumplir que: 𝜌𝜕𝜙+ 𝜕𝑡 +𝜌𝑈∞𝜕𝜙+ 𝜕𝑥 +𝑝+=0 𝜌𝜕𝜙− 𝜕𝑡 +𝜌𝑈∞𝜕𝜙− 𝜕𝑥 +𝑝−=0 Si se restan ambas ecuaciones y se tienen cuenta que la presión en el extradós y en el intradós son iguales: 𝜕(𝜙+−𝜙−) 𝜕𝑡 +𝑈∞𝜕(𝜙+−𝜙−) 𝜕𝑥 =0 Teniendo en cuenta la antisimetría del campo de velocidades en el problema sustentador, Gordilllo y Riboux, 2012, [1], 𝜙+=−𝜙−, se puede reescribir la condición anterior para cualquiera de las superficies extradós o intradós. Suprimiendo el superíndice + en el potencial, se tiene finalmente. 𝜕𝜙 𝜕𝑡+𝑈∞𝜕𝜙 𝜕𝑥=𝐷𝜙 𝐷𝑡=0 (3.6.1) Esta ecuación es la que se va a imponer en el borde de salida. Esto quiere decir que la variación por unidad de tiempo del potencial de velocidades siguiendo a la partícula fluida es nula. En el caso estacionario, esta condición se traduce en que en el borde de salida hay un punto de remanso (en el caso de borde de salida anguloso), o bien que la velocidad en el borde de salida tiene la misma dirección y sentido tanto en el extradós como en el intradós (caso de borde de salida en retroceso). En el caso estacionario la condición anterior queda: 𝜕𝜙 𝜕𝑥=0 (3.6.2) 36 37 4. SOLUCIÓN GENERAL DEL PROBLEMA Se define el potencial 𝜓0 de una solución elemental de la ecuación de Laplace para el caso tridimensional correspondiente a una fuente localizada en el punto (𝑥0,𝑦0,𝑧0) como: 𝜓0=1 √(𝑥−𝑥0)2+(𝑦−𝑦0)2+(𝑧−𝑧0)2 (4.1) Se puede comprobar que este tipo de solución cumple la ecuación de Laplace ∇2𝜓0=0. El potencial de velocidades perturbadas, debe de cumplir también la ecuación de Laplace ∇2𝜙′=0. Recordemos que se tiene que cumplir que: ∇2𝜙′=0, 𝑣′=∇𝜙′, → ∇2𝑣′=0 Donde se ha hecho uso de la intercambiabilidad de las derivadas. Debido a esto se debe de cumplir la siguiente ecuación: 𝜓0∇2𝑣′−𝑣′ ∇2𝜓0=0 (4.2) A continuación se procede a integral esta ecuación en un volumen de control Ω𝑐 limitado por la superficie Σ𝑐=Σ∞∪Σ𝜖∪Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− . Dichas superficies se observan en 2D en la siguiente figura: ∫(𝜓0∇2𝑣−𝑣 ∇2𝜓0)𝑑𝜔= Ω𝑐∫∇·(𝜓0∇2𝑣−𝑣 ∇2𝜓0)𝑑𝜔= Ω𝑐 FIGURA 2: SUPERFICIES DE INTEGRACIÓN 38 =∫𝜓0𝑛 󰇍  ·∇𝑣′−𝑣′𝑛 󰇍  ·∇𝜓0 𝑑𝜎=0 Σ𝑐 La integral a lo largo de Σ∞ es nula, y las integrales a lo largo de Σ𝜖 y Σ𝑎𝑙𝑎 pueden escribirse de la siguiente forma: ∫𝜓0𝑛 󰇍  ·∇𝑣′−𝑣′𝑛 󰇍  ·∇𝜓0𝑑𝜎= Σ𝜖∪Σ𝑎𝑙𝑎 =∫𝜓0𝑛 󰇍  ·∇𝑣′−𝑣′𝑛 󰇍  ·∇𝜓0 𝑑𝜎+ Σ𝑎𝑙𝑎 ∫𝜓0𝑛 󰇍  𝑒−·∇𝑣′−𝑣′𝑛 󰇍  𝑒−·∇𝜓0 𝑑𝜎+ Σ𝜖− +∫𝜓0𝑛 󰇍  ·∇𝑣′−𝑣′𝑛 󰇍  ·∇𝜓0 𝑑𝜎+∫𝜓0(−𝑛 󰇍  𝑒−)·∇𝑣′−𝑣′(−𝑛 󰇍  𝑒−)·∇𝜓0𝑑𝜎= Σ𝜖− Σ𝜖 =𝐼𝑎𝑙𝑎+𝐼𝜖 Donde Σ𝜖 es la semiesfera superior que rodea la singularidad y cuya normal positiva se define hacia el interior de la esfera, y Σ𝜖− es la semiesfera inferior que rodea la singularidad y cuya normal positiva se define hacia el exterior de dicha esfera. La integral 𝐼𝜖 se puede calcular de la siguiente manera. 𝐼𝜖=lim ϵ→0∫(1ϵ(−𝑒𝑟)·∇𝑣′(𝑥)−𝑣′(𝑥)(𝑒𝑟)·(−𝑒𝑟) 𝜖2)𝜖2𝑑Ω=−4𝜋𝑣′(𝑥0 󰇍 󰇍 󰇍 󰇍  ) 4𝜋 0 (4.3) Donde en este caso 𝑑Ω=𝑑S/r2 es el diferencial de ángulo sólido. Tras este cálculo la ecuación de Green puede escribirse de la siguiente forma: 4𝜋𝑣′(𝑥0)=∫𝜓0𝑛 󰇍  ·∇𝑣′−𝑣′𝑛 󰇍  ·∇𝜓0𝑑𝜎 Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− (4.4) Proyectando sobre el vector unitario 𝑘 󰇍  . 4𝜋𝑤 󰇍 󰇍  ′(𝑥0 󰇍 󰇍 󰇍 󰇍  )=∫𝜓0𝑛 󰇍  ·∇𝑤 󰇍 󰇍  ′−𝑤 󰇍 󰇍  ′𝑛 󰇍  ·∇𝜓0𝑑𝜎 Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− (4.5) Debido a que la antisimetría 𝑤′(𝑥,𝑦,𝑧=0+)=𝑤′(𝑥,𝑦,𝑧=0−) y a que las normales al ala y a la estela son aproximadamente 𝑘 󰇍  y -𝑘 󰇍  para el extradós y para el intradós, se tiene pues que: ∫𝑤′(𝑥,𝑦,𝑧=0 )𝑛 󰇍  ·∇𝜓0𝑑𝜎=0 Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− FIGURA 3: NORMALES EN EL ENTORNO DE LA SINGULARIDAD 39 ∫𝑤′(𝑥,𝑦,𝑧=0 )𝑛 󰇍  ·∇𝜓0𝑑𝜎=0 Σ𝑎𝑙𝑎+∪Σ𝑎𝑙𝑎− Luego se llega a la siguiente ecuación 4𝜋𝑤 󰇍 󰇍  ′(𝑥0)=∫𝜓0𝑛 󰇍  ·∇𝑤 󰇍 󰇍  ′𝑑𝜎= Σ𝑎𝑙𝑎∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− =−∫𝜓0𝜕𝑤′+ 𝜕𝑧 𝑑𝜎++∫𝜓0𝜕𝑤′− 𝜕𝑧 𝑑𝜎− Σ𝑎𝑙𝑎−∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+ Aplicando la ecuación de continuidad tanto en el extradós como en el intradós se tiene: ∇·𝑣′+ 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  =𝜕𝑢′+ 𝜕𝑥 +𝜕𝑣′+ 𝜕𝑦 +𝜕𝑤′+ 𝜕𝑧 =0 → 𝜕𝑤′+ 𝜕𝑧 =−𝜕𝑢′+ 𝜕𝑥 −𝜕𝑣′+ 𝜕𝑦 ∇·𝑣′− 󰇍 󰇍 󰇍 󰇍 󰇍 󰇍  =𝜕𝑢′− 𝜕𝑥 +𝜕𝑣′− 𝜕𝑦 +𝜕𝑤′− 𝜕𝑧 =0 → 𝜕𝑤′− 𝜕𝑧 =−𝜕𝑢′− 𝜕𝑥 −𝜕𝑣′− 𝜕𝑦 Dado que la placa parte del reposo, la circulación Γ es idénticamente nula en cualquier superficie cerrada que englobe al ala como a la estela según el Teorema de Bjerknes-Kelvin, el cual dice: 𝐷Γ 𝐷𝑡=0 (4.6) Tomando como superficie, una compuesta por el ala y la estela, se tiene: Γ=∫ ∫ 𝛾(𝑥,𝑦,𝑡) 𝑑𝑥 𝑑𝑦=∫ ∫ (𝑢′+−𝑢′−) 𝑑𝑥 𝑑𝑦=∫ ∫ 𝜕(𝜙+−𝜙−) 𝜕𝑥 𝑑𝑥 𝑑𝑦 𝑥𝑒𝑠𝑡 0 𝑏 2 −𝑏 2 𝑥𝑒𝑠𝑡 0 𝑏 2 −𝑏 2 𝑥𝑒𝑠𝑡 0 𝑏 2 −𝑏 2 Γ=∫(𝜙+(𝑥𝑒𝑠𝑡,𝑦,𝑡)−𝜙− 𝑏 2 −𝑏 2(𝑥𝑒𝑠𝑡,𝑦,𝑡))𝑑𝑦=0 Donde el subíndice “est” hace referencia a la estela y 𝑥𝑒𝑠𝑡 es el final de la estela y la variable 𝛾(𝑥,𝑦,𝑡), es la densidad de circulación. Se ha utilizado que el potencial es 0 en el borde de ataque, condición que se explicará posteriormente, cuando se analicen las condiciones de contorno del problema. Como esta integral es nula para todo instante de tiempo, de ello se deduce que el integrando no puede depender del tiempo. Esto implica que la variación con el tiempo del potencial en la coordenada donde acaba la estela es nula, y por tanto las velocidades horizontales son iguales. Utilizando las siguientes igualdades obtenidas de la derivación de un producto: 𝜓0𝜕𝑢′ 𝜕𝑥=𝜕(𝜓0𝑢′) 𝜕𝑥 −𝑢′𝜕𝜓0 𝜕𝑥 40 𝜓0𝜕𝑣′ 𝜕𝑦=𝜕(𝜓0𝑣′) 𝜕𝑦 −𝑣′𝜕𝜓0 𝜕𝑦 Y sustituyendo en la ecuación de Green, se tiene que: 4𝜋𝑤 󰇍 󰇍  ′(𝑥0)=∫(𝜕(𝜓0𝑢′+) 𝜕𝑥 −𝑢′+𝜕𝜓0 𝜕𝑥 +𝜕(𝜓0𝑣′+) 𝜕𝑦 −𝑣′+𝜕𝜓0 𝜕𝑦)𝑑𝜎+− Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+ −∫(𝜕(𝜓0𝑢′−) 𝜕𝑥 −𝑢′−𝜕𝜓0 𝜕𝑥+𝜕(𝜓0𝑣′−) 𝜕𝑦 −𝑣′−𝜕𝜓0 𝜕𝑦)𝑑𝜎− Σ𝑎𝑙𝑎−∪Σ𝑒𝑠𝑡𝑒𝑙𝑎− Para este razonamiento se va a tener en cuenta que las integrales de superficie se pueden hacer indistintamente primero en la dirección x y luego en la dirección y o viceversa. Tal como se ha escrito la expresión anterior se observa que los términos son iguales en ambas integrales, simplemente cambiando la notación de extradós a intradós. Entonces: - La suma de las integrales del primer término es 0, debido a que al hacer la integral en dirección x se tiene que en el borde de ataque dichas velocidades son 0, y a que al final de la estela se cumple 𝑢′+=𝑢′− . - La suma de las integrales del tercer término es 0 debido a que al hacer la integral en dirección y se debe de tener en cuenta que las velocidades 𝑣′ son antisimétricas respecto al plano z=0, siendo entonces dichas velocidades idénticas pero de signo contrario en los bordes marginales. Eliminando estas integrales de la expresión anterior y haciendo uso de nuevo de la antisimetría del campo de velocidades se llega a: 2𝜋𝑤 󰇍 󰇍  ′(𝑥0)=−∫(𝑢′+𝜕𝜓0 𝜕𝑥+𝑣′+𝜕𝜓0 𝜕𝑦)𝑑𝜎+ Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+ Ecuación que finalmente queda: 𝑤 󰇍 󰇍  ′(𝑥0)=𝑈∞𝜕𝑧𝑝 𝜕𝑥(𝑥0,𝑦0,𝑡)+𝜕𝑧𝑝 𝜕𝑡(𝑥0,𝑦0,𝑡)=−1 2𝜋∫𝑣′(𝑥)·∇𝜓0 𝑑𝜎+ Σ𝑎𝑙𝑎+∪Σ𝑒𝑠𝑡𝑒𝑙𝑎+ (4.7) Esta es la solución general del problema aerodinámico tridimensional. Se observa que dicha expresión se corresponde con la ecuación de Biot y Savart. La variable 𝑣′(𝑥) indica las componentes de la velocidad en el plano Z=0. El lector interesado puede encontrar más información en Gordillo y Riboux, 2012, [1]. De esta forma se observa que la componente vertical del campo de velocidades perturbado en cualquier punto del dominio fluido puede ser expresada como la superposición del campo de velocidades creado por una distribución continua de torbellinos situados sobre la superficie del ala y de la estela, y orientados según la envergadura, con intensidad por unidad de longitud 41 2𝑢′+(𝑥,𝑦), y una distribución continua de torbellinos situados sobre ala y estela, de intensidad por unidad de longitud 2𝑣′+(𝑥,𝑦) orientados en la dirección del eje x. Trabajos anteriores tratan sobre la resolución de esta ecuación integral. Sin embargo, dichos estudios demuestran que el sistema que se obtiene es inestable desde el punto de vista numérico. Véase Flores Caballero, 2016 [4]. El método que se va a implementar corrige esta inestabilidad mediante la transformación de dicha ecuación integral, ya que en el integrando aparecerá como incógnita el valor del potencial perturbado en lugar de su gradiente. 48 De esta forma el denominador de la ecuación se puede desarrollar como sigue: |𝑅 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  | =|(𝑅 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )+(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  )| =|𝑟𝑒𝑟 󰇍 󰇍 󰇍  +(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  )| = =|𝑟𝑒𝑟 󰇍 󰇍 󰇍  +(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  )| =(𝑟2+|𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  |2+2𝑟𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  ))1/2= =((𝑟+𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  ))2+|𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  |2−(𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  ))2)1/2= =(𝑈2+𝑎2)1/2 Donde x y a vienen dadas por las siguientes expresiones (5.8) y (5.9): 𝑈=𝑟+𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  ) (5.8) 𝑎2=|𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  |2−(𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅1 󰇍 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  ))2 (5.9) Esta nomenclatura será utilizada más adelante. Utilizando las nuevas variables 𝜃 y r, se llega a la siguiente conclusión. 𝑥−𝑥1=𝑟𝑐𝑜𝑠(𝜃) 𝑦−𝑦1=𝑟𝑠𝑖𝑛(𝜃) Sustituyendo en la expresión anterior, desarrollando el denominador de la integral y utilizando la definición de las variables U y 𝑎2 se llega a lo siguiente. 𝑤𝑖=∫𝜙(𝑥𝑖) |𝑅𝑖 󰇍 󰇍 󰇍  −𝑅0 󰇍 󰇍 󰇍 󰇍  |3𝑑𝜎𝑖 Σi= 𝜙1∫[ 1 (𝑈2+𝑎2)3/2+𝐶1𝑟𝑐𝑜𝑠(𝜃) (𝑈2+𝑎2)3/2+𝐷1𝑟𝑠𝑖𝑛(𝜃) (𝑈2+𝑎2)3/2]𝑑𝜎𝑖 Σi+ + 𝜙2∫[ 𝐶2𝑟𝑐𝑜𝑠(𝜃) (𝑈2+𝑎2)3/2+𝐷2𝑟𝑠𝑖𝑛(𝜃) (𝑈2+𝑎2)3/2]𝑑𝜎𝑖 Σi+ + 𝜙3∫[ 𝐶3𝑟𝑐𝑜𝑠(𝜃) (𝑈2+𝑎2)3/2+𝐷3𝑟𝑠𝑖𝑛(𝜃) (𝑈2+𝑎2)3/2]𝑑𝜎𝑖 Σi Donde U=U(r), y los límites de la integral van desde r=0, hasta r=r(𝜃). Desarrollando la integral doble y sacando los términos que no dependen de r fuera de la integral se llega a lo siguiente. 𝜙1∫ [∫ 𝑟𝑑𝑟 (𝑈2+𝑎2)3/2 𝑟(𝜃) 0+(𝐶1𝑐𝑜𝑠(𝜃)+𝐷1𝑠𝑖𝑛(𝜃))∫𝑟2𝑑𝑟 (𝑈2+𝑎2)3/2 𝑟(𝜃) 0]𝑑𝜃+ 𝜃3 θ2 49 + 𝜙2∫ [(𝐶2𝑐𝑜𝑠(𝜃)+𝐷2sin (𝜃))∫𝑟2𝑑𝑟 (𝑈2+𝑎2)3 2 𝑟(𝜃) 0]𝑑𝜃+ 𝜃3 θ2 + 𝜙3∫ [(𝐶3𝑐𝑜𝑠(𝜃)+𝐷3𝑠𝑖𝑛(𝜃))∫𝑟2𝑑𝑟 (𝑈2+𝑎2)3/2 𝑟(𝜃) 0]𝑑𝜃 𝜃3 θ2 Como se puede observar, las integrales a realizar son similares entre ellas. Se va a pasar a analizarlas una a una individualmente. La primera se muestra a continuación: 𝑆1=∫𝑟𝑑𝑟 (𝑈2+𝑎2)3/2 𝑟(𝜃) 0=∫𝑈−𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ) (𝑈2+𝑎2)3/2 𝑑𝑈 𝑈(𝜃) 𝑈0 Donde se ha hecho el cambio de variable mostrado anteriormente. Operando con esta expresión se obtiene lo siguiente: ∫𝑈−𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ) (𝑈2+𝑎2)3/2 𝑑𝑈 𝑈(𝜃) 𝑈0=∫𝑈 (𝑈2+𝑎2)3/2𝑑𝑈 𝑈(𝜃) 𝑈0−𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )∫ 𝑑𝑈 (𝑈2+𝑎2)3/2 𝑈(𝜃) 𝑈0 𝑆1=−𝑆0−𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )𝑆3 Donde: 𝑆0=−∫𝑈 (𝑈2+𝑎2)3 2𝑑𝑈 𝑈(𝜃) 𝑈0=1 ((𝑈(𝜃))2+𝑎2)1 2−1 (𝑈02+𝑎2)1 2 𝑆3=∫𝑑𝑈 (𝑈2+𝑎2)3/2 𝑈(𝜃) 𝑈0=1 𝑎2[𝑈(𝜃) (𝑈(𝜃)2+𝑎2)1 2−𝑈0 (𝑈02+𝑎2)1 2] A continuación se va a analizar el otro tipo de expresión que aparece en la integral. 𝑆2=∫𝑟2𝑑𝑟 (𝑈2+𝑎2)3/2 𝑟(𝜃) 0=∫(𝑈−𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ))2 (𝑈2+𝑎2)3/2 𝑑𝑈= 𝑈(𝜃) 𝑈0 =∫𝑈2−2𝑈𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )+(𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ))2 (𝑈2+𝑎2)3/2 𝑑𝑈= 𝑈(𝜃) 𝑈0 =∫𝑈2−2𝑈𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )+(𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ))2±𝑎2 (𝑈2+𝑎2)3/2 𝑑𝑈= 𝑈(𝜃) 𝑈0 Donde en el último paso se ha sumado y restado el termino de 𝑎2. A continuación, desarrollando la integral y utilizando las integrales que se han calculado anteriormente, se obtiene. 50 𝑆2=∫𝑑𝑈 (𝑈2+𝑎2)1/2+2𝑈𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )𝑆0+ 𝑈(𝜃) 𝑈0[2(𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ))2−|𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  |2]𝑆3 Por comodidad y al igual que se ha hecho anteriormente, se va a utilizar la siguiente nomenclatura. 𝑆4=∫𝑑𝑈 (𝑈2+𝑎2)1/2=𝐿𝑛(√𝑎2+𝑈(𝜃)2+𝑈(𝜃) √𝑎2+𝑈02+𝑈0) 𝑈(𝜃) 𝑈0 𝑆2=𝑆4+2𝑈𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  )𝑆0+[2(𝑒𝑟 󰇍 󰇍 󰇍  ·(𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  ))2−|𝑅0 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  |2]𝑆3 Sustituyendo las definiciones de 𝑆1 y 𝑆2 en la ecuación integral para el triángulo bajo estudio se tendría que: 𝑤𝑖= 𝜙1∫[𝑆1+(𝐶1𝑐𝑜𝑠(𝜃)+𝐷1sin (𝜃))𝑆2]𝑑𝜃+ 𝜃3 θ2 + 𝜙2∫[(𝐶2𝑐𝑜𝑠(𝜃)+𝐷2sin (𝜃))𝑆2]𝑑𝜃+ 𝜃3 θ2 + 𝜙3∫[(𝐶3𝑐𝑜𝑠(𝜃)+𝐷3sin (𝜃))𝑆2]𝑑𝜃 𝜃3 θ2 A continuación es importante recordar que 𝑈=𝑈(𝑟(𝜃)), donde 𝑈(𝑟) es conocida, pero aún hay que determinar 𝑟(𝜃). Para ello hay que utilizar el siguiente triángulo. 51 Donde se puede ver que 𝑟(𝜃) es la distancia desde un punto del lado 3-2 al vértice 1 del triángulo. Aplicando relaciones entre los distintos ángulos, se pueden obtener los ángulos interiores de la región indicada. Aplicando el Teorema del seno y despejando se obtiene la expresión 𝑟(𝜃) buscada. 𝑟(𝜃) sin (𝜋+𝜃2−𝜃4)=|𝑅2 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  | sin(𝜃4−𝜃) 𝑟(𝜃)=sin (𝜃4−𝜃2)|𝑅2 󰇍 󰇍 󰇍 󰇍  −𝑅1 󰇍 󰇍 󰇍 󰇍  | sin(𝜃4−𝜃) FIGURA 8: ÁNGULOS EMPLEADOS EN LA DISCRETIZACIÓN. GRÁFICO EXPLICATIVO II FIGURA 9: DEFINICIÓN DE 𝒓(𝜽) 52 Una vez se haya hecho la integral en la variable 𝑟, se procede a la integración en la variable 𝜃, la cual se hace de forma numérica mediante un método de colocación. Una vez hecha la integración para la primera ecuación, ya se tienen los valores de los coeficientes de influencia que acompañan a los potenciales de velocidades. Estos potenciales van colocados en los vértices de los triángulos en los que se ha dividido el ala y la estela. En el caso en que se tenga que integrar en el entorno de la singularidad, se realizará el siguiente mallado: Y se supondrá que el potencial es de la siguiente forma: 𝜙′(𝑥)=𝜙1′=𝜙0′ Es decir, en el entorno de la singularidad (8 triángulos) el potencial es constante e igual al potencial de la singularidad. Esto evita la enorme inestabilidad numérica que tendrían las ecuaciones si se dejase evolucionar el potencial como en el resto de los casos. La integral no es divergente en el entorno de la singularidad sólo si el potencial es constante. La condición de contorno de impenetrabilidad discretizada queda de la siguiente forma: 𝑤(𝑥𝑖 󰇍 󰇍 󰇍  )=[𝐴1,1,𝐴1,2,…𝐴1,𝑖,…𝐴1,(𝑁+2)(𝑀+2)]· [ 𝜙1,1 𝜙1,2 … 𝜙1,𝑖 … 𝜙1,(𝑁+2)(𝑀+2) ] + +[𝐸1,1,𝐸1,2,…𝐸1,𝑖,…𝐸1,(𝑁+2)(𝑀 𝑒𝑠𝑡𝑒𝑙𝑎)]· [ 𝜙𝐸1,1 𝜙𝐸1,2 … 𝜙𝐸1,𝑖 … 𝜙𝐸1,(𝑁+2)(𝑀 𝑒𝑠𝑡𝑒𝑙𝑎) ] FIGURA 10: MALLADO EN EL ENTORNO DE LA SINGULARIDAD 53 Donde 𝜙𝑖,𝑗 indica el potencial en el nodo i,j y 𝐴𝑖,𝑗 y 𝐸𝑖,𝑗 son los correspondientes coeficientes de inflencia, calculados empleando las expresiones analíticas de las integrales halladas anteriormente. 54 55 6. PROBLEMA ESTACIONARIO. RESOLUCIÓN Y CONDICIONES DE CONTORNO Una vez realizada la formulación de la condición de impenetrabilidad en su forma general, en este apartado se va a pasar a exponer las condiciones de contorno para el problema estacionario y su resolución. Para cerrar el problema, se deben de exponer un total de ecuaciones igual al número de incógnitas, con el objetivo de que la matriz de coeficientes sea cuadrada e invertible. Cómo se ha visto anteriormente que en cualquier punto del dominio fluido donde no haya sólido, se debe cumplir que el potencial de velocidades perturbadas debe de ser 0 al serlo en el infinito aguas arriba, en particular, sobre el borde de ataque y sobre los bordes marginales se debe de cumplir que, en virtud de la ecuación de Euler-Bernouilli linealizada, que expresa que la derivada sustancial del potencial perturbado debe de ser 0, se concluye que dicho potencial es 0 en los bordes de ataque y marginales.  Borde de ataque y bordes marginales: 𝜙′(𝑥𝑏𝑎,𝑦)=𝜙′(𝑥,𝑏2)=𝜙′(𝑥,−𝑏2)=0 Por otra parte, en el borde de salida, la condición de Kutta-Joukowski expresa en el caso estacionario que:  Borde de salida: 𝜕𝜙′ 𝜕𝑥=0→𝜙(𝑥𝑏𝑠−1,𝑦)=𝜙(𝑥𝑏𝑠,𝑦)=𝜙(𝑥𝑏𝑠+1=𝑥𝐸1,𝑦) Esta condición implica que para una línea de coordenada y, el valor del potencial de velocidades en el borde de salida en esa línea y es el mismo en el borde de salida, en la discretización en x anterior al borde de salida y a lo largo de toda la estela. Gráficamente se puede ver en la figura siguiente: 56 Se observa que en como el potencial es el mismo tanto en el borde de salida, como en el punto anterior y como en toda la estela a lo largo de cada línea y, no tiene sentido realiza un mallado a lo largo de la estela para el caso estacionario. Luego en este caso, el mallado es mucho más simple: FIGURA 11: POTENCIAL EN EL BORDE DE SALIDA Y EN LA ESTELA EN RÉGIMEN ESTACIONARIO PARA UN VALOR DE Y DADO FIGURA 12: MALLADO TRIANGULAR EN RÉGIMEN ESTACIONARIO 57 A la hora de montar las ecuaciones este recinto se puede subdividir de la siguiente forma Únicamente se debe de imponer la condición de contorno de impenetrabilidad en la zona interior al recuadro amarillo, es decir, un total de N*M ecuaciones. El potencial en el borde de salida y en la estela (recuadro azul) debe de ser el mismo que en los puntos anteriores al borde de salida, y el potencial en los bordes de ataque y marginales (recuadros verdes) es 0. Así pues se tiene un total de N*(M+1) incógnitas sobre el ala, y N*(M+1) ecuaciones divididas en N*M condiciones de impenetrabilidad, más N condiciones de contorno de Kutta sobre el borde de salida. Esto permite formar una matriz cuadrada e invertible con la que resolver el problema. La ecuación quedaría de la siguiente forma general para el problema estacionario (se omiten los términos que acompañan a potenciales nulos por la condición de contorno). 𝑤′(𝑥𝑖 󰇍 󰇍 󰇍  )=𝑤𝑖′=[𝐴𝑖,2,𝐴𝑖,3,…𝐴𝑖,𝑀+2,…𝐴𝑖,(𝑁+1)(𝑀+2)]· [ 𝜙2,2 𝜙2,3 … 𝜙2,𝑀+2 … 𝜙2,(𝑁+1)(𝑀+2) ] + +[𝐸𝑖,2,𝐸𝑖,3,…𝐸𝑖,𝑖,…𝐸𝑖,(𝑁+1)·]· [ 𝜙𝐸2 𝜙𝐸2 … 𝜙𝐸𝑖 … 𝜙𝐸(𝑁+1) ] FIGURA 13: CONDICIONES DE CONTORNO PARA EL POTENCIAL EN RÉGIMEN ESTACIONARIO 64 instante de tiempo simplemente mediante la “convección” aguas debajo de los potenciales obtenidos en el borde de salida. En el anexo se proporciona el programa de Matlab que resuelve el problema estacionario y no estacionario para geometrías genéricas. 65 8. SOLUCIÓN EN RÉGIMEN ESTACIONARIO. 8.1. Alas Rectangulares A continuación se va a exponer la solución para alas rectangulares en régimen estacionario, y para distintos alargamientos. Para obtener una solución asequible es necesario alcanzar un equilibrio entre el número de divisiones del mallado y el tiempo computacional, de manera que se intente alcanzar la mejor solución (o una solución con precisión suficiente) en un tiempo de cálculo razonable. Para alcanzar este objetivo es preciso recordar que la distribución de presiones en el extradós del ala sufre unos gradientes muy acentuados en las proximidades del borde de ataque y de los bordes marginales. Este fenómeno lleva a pensar que es necesario refinar el mallado en estas zonas, con objetivo de obtener unos valores razonables del coeficiente de sustentación y de la pendiente de la curva de sustentación (Se supondrá en todo momento variación lineal del coeficiente de sustentación con el ángulo de ataque). El mallado en el eje x e y se obtiene con las siguientes expresiones (teniendo en cuenta el número de divisiones): 𝑚𝑦=0.5𝑏·cos((𝑗−1)𝜋 𝑁+1) 𝑐𝑜𝑛 𝑗∈[1,𝑁+2] 𝑚𝑥=(1−cos((𝑖−1)𝜋 2(𝑀−1))) 𝑐𝑜𝑛 𝑖∈[1,𝑀] El mallado para el caso de alas en flecha es similar. Se muestra un ejemplo de cómo serían las divisiones del mallado. Se ha presentado el ejemplo para que se observe el alargamiento del ala y se pueda intuir la longitud de la estela (la estela FIGURA 17: DISCRETIZACIÓN EN ALAS RECTANGULARES 66 llegaría más lejos, pero la idea con esta imagen es mostrar al lector las dimensiones del problema). En la figura anterior se observa cómo el mallado es más fino en los puntos del borde de ataque y marginales para captar bien los gradientes. Resolviendo el problema para un ángulo de ataque de 5º, un alargamiento de 5 (b=5, c=1) y estrechamiento nulo, se obtienen las siguientes distribuciones del potencial sobre la superficie del ala. Si se realiza la derivación del potencial respecto a la coordenada x, se puede obtener el campo de velocidades sobre el extradós. Estas dos figuras muestran el enorme gradiente de velocidades que existe tanto en los bordes marginales, como en el borde de ataque y por tanto justifican el mallado realizado. La pendiente de la curva de sustentación obtenida para esta configuración es de: FIGURA 18: POTENCIAL SOBRE EL EXTRADÓS DE ALAS RECTANGULARES FIGURA 19: VELOCIDAD PERTURBADA EN EL EJE X SOBRE EL EXTRADÓS DE ALAS RECTANGULARES. 67 𝐶𝐿𝛼 Nx Ny AR 3,8741 15 21 5 TABLA 1: RESULTADOS PARA ALA RECTANGULAR EN RÉGIMEN ESTACIONARIO A continuación se ha realizado el mismo cálculo pero aumentando el mallado hasta un valor de Nx=25 y Ny=41 para observar las diferencias. Se muestra a continuación una figura del campo de velocidades para este caso, de forma que se aprecie mejor el gradiente que se produce en el borde de ataque y en los bordes marginales. Los resultados para este caso son: 𝐶𝐿𝛼 Nx Ny AR 3,9096 25 41 5 TABLA 2: RESULTADOS PARA ALA RECTANGULAR EN RÉGIMEN ESTACIONARIO. MALLADO REFINADO El tiempo computacional ha sido mayor para este caso, y la mejora se nota en la tercera cifra significativa. Se puede observar que la pendiente de la curva de sustentación tiende a aumentar conforme el mallado se hace más fino. Esto se debe a que se capta mejor las evoluciones de las variables fluidas en el entorno de los bordes marginales y del borde de ataque, permitiendo que la integración de lugar a un resultado más preciso. De aquí en adelante se realizará el mallado con los siguientes valores de Nx=15 y Ny=8*b+1. La razón de estos valores es porque cuando se resuelva el problema no estacionario, el tiempo FIGURA 20: VELOCIDAD PERTURBADA EN EL EJE X SOBRE EL EXTRADÓS DE ALAS RECTANGULARES. REFINAMIENTO 68 computacional aumenta enormemente para mallados muy grandes si se quieren obtener valores que mejoran poco en precisión con respecto a estos. De forma que con este método acotamos el error a la tercera cifra significativa mientras mantenemos un tiempo razonable. La razón de la ecuación Ny=8*b+1, estriba en que en el análisis de alas con gran alargamiento, se deberán hacer más divisiones en la coordenada y debido a que al aumentar la envergadura, también lo hacen los paneles en los que se divide el ala en esa dirección y por tanto se perdería precisión si no se aumenta el mallado. 8.2. Alas con Flecha En este apartado se trata la solución del problema estacionario para el caso de alas con flecha. De esta forma, dando valores de S, AR, 𝛹 y E, queda definida el ala. Para el caso de alas en flecha, se probó a utilizar un mallado similar al caso de alas rectangulares, sin embargo la precisión adoptada en comparación con los valores típicos de la literatura no era suficiente. La razón estriba en que las alas en flecha poseen un gradiente de velocidades en las proximidades de la cuerda raíz del ala, que resultaba ser la zona donde el mallado anterior es menos preciso (elementos más grandes, captaban mal la evolución). Para la comprobación de que este código era efectivo, puede comprobar el lector interesado que si se utiliza el código para el ala en flecha con las particularidades de E=1 y ángulo de flecha nulo, se vuelve a obtener los valores que se obtuvieron para el ala con forma en planta rectangular. Para solventar este problema se recurrió a un refinado también por dicha zona similar a lo que se había hecho para los bordes de ataque y marginales. Se muestra el mallado particular realizado para las alas en flecha: FIGURA 21: DISCRETIZACIÓN EN ALAS CON FLECHA 69 El potencial de velocidades sobre este tipo de alas tiene la siguiente forma: Derivando el potencial, se puede obtener el campo de velocidades perturbadas sobre la superficie de extradós, el cual se muestra en la siguiente figura: En esta figura se puede observar el gradiente de velocidades que se produce en las proximidades de la cuerda raíz del ala. 𝐶𝐿𝛼 Nx Ny AR 3,5368 15 41 5 TABLA 3: RESULTADOS PARA ALA CON FLECHA 30º EN RÉGIMEN ESTACIONARIO. Se observa que el valor de la pendiente de la curva de sustentación para un ala en flecha es menor que para el ala rectangular, como era de esperar. FIGURA 22: POTENCIAL SOBRE EL EXTRADÓS DE ALAS CON FLECHA FIGURA 23: VELOCIDAD PERTURBADA EN EL EJE X SOBRE EL EXTRADÓS DE ALAS CON FLECHA. 70 71 9. SOLUCIÓN EN RÉGIMEN NO ESTACIONARIO. 9.1. Alas Rectangulares A continuación se va a mostrar los resultados obtenidos para alas rectangulares en régimen no estacionario. Para ello, se va a realizar el mismo análisis, con el mismo mallado en el ala, pero para régimen no estacionario. El incremento de tiempo Δ𝑡 debe de ser inferior a 𝑐/𝑈∞, para que el efecto no estacionario se observe correctamente. El tiempo de cálculo se ha restringido a 200 𝑐/𝑈∞, con un Δ𝑡 de 0.25𝑐/𝑈∞, entre cada iteración. Las ecuaciones han sido integradas en el tiempo utilizando un método de Euler de primer orden. El valor obtenido se muestra a continuación junto con la curva que indica la evolución que ha seguido el valor de la pendiente de la curva de sustentación desde el instante inicial hasta el final para el caso de un ala rectangular con relación de aspecto AR=5 y suponiendo que la velocidad adimensional en el infinito es una función escalón en el tiempo. 𝐶𝐿𝛼 Nx Ny AR 3,8788 15 41 5 TABLA 4: RESULTADOS PARA ALA RECTANGULAR EN RÉGIMEN NO ESTACIONARIO Cabe resaltar el hecho de que la pendiente de la curva de sustentación alcanza el valor obtenido en régimen estacionario anteriormente. FIGURA 24: EVOLUCIÓN DE LA PENDIENTE DE LA CURVA DE SUSTENTACIÓN CON EL TIEMPO. ALA RECTANGULAR 72 9.2. Alas con Flecha Para las alas en flecha se ha seguido un procedimiento similar al de las alas rectangulares, con la particularidad de hacer un mallado más fino en la zona cercana a la cuerda raíz del ala. Una vez se tiene implementado el código para la resolución no estacionaria de alas rectangulares, simplemente se modifica la geometría para realizar el estudio del comportamiento no estacionario para alas en flecha. Se ha realizado de nuevo una restricción en el tiempo de cálculo de 200 𝑐/𝑈∞, con un Δ𝑡 de 0.25𝑐/𝑈∞. Esto da paneles de longitud 0.25 en la estela y da valores lo suficientemente precisos aunque con un tiempo computacional elevado. 𝐶𝐿𝛼 Nx Ny AR 3,5303 15 41 5 TABLA 5: RESULTADOS PARA ALA CON FLECHA 45º EN RÉGIMEN ESTACIONARIO Para que el lector tenga una idea de las dimensiones que tiene la estela no estacionaria analizada para el case de un ala en flecha se adjunta la figura siguiente: FIGURA 25: EVOLUCIÓN DE LA PENDIENTE DE LA CURVA DE SUSTENTACIÓN CON EL TIEMPO. ALA CON FLECHA 73 Cada división en el eje x corresponde a un instante de tiempo. En total se tienen 200 instantes de tiempo en divisiones de 0.25 s. FIGURA 26: ESTELA NO ESTACIONARIA EN CASO DE ALA EN FLECHA 80 El error obtenido entre el cálculo estacionario y el obtenido del gráfico se muestran a continuación: ALA RECTANGULAR AR 𝐶𝐿𝛼𝐸𝑠𝑡𝑎𝑐𝑖𝑜𝑛𝑎𝑟𝑖𝑜 KyP Error (%) 4 3,5656 3,6 0,95 5 3,8741 3,9 0,66 6 4,1321 4,2 1,61 7 4,3236 4,3 0,54 TABLA 9: ERROR COMETIDO EN ALAS RECTANGULARES EN RÉGIMEN ESTACIONARIO ALA FLECHA 30º AR 𝐶𝐿𝛼𝐸𝑠𝑡𝑎𝑐𝑖𝑜𝑛𝑎𝑟𝑖𝑜 KyP Error (%) 4 3,2911 3,35 1,75 5 3,5368 3,7 4,41 6 3,7222 3,85 3,31 7 3,8712 4 3,22 TABLA 10: ERROR COMETIDO EN ALAS CON 𝚿=𝟑𝟎º EN RÉGIMEN ESTACIONARIO ALA FLECHA 45º AR 𝐶𝐿𝛼𝐸𝑠𝑡𝑎𝑐𝑖𝑜𝑛𝑎𝑟𝑖𝑜 KyP Error (%) 4 2,922 3 2,6 5 3,0942 3,2 3,30625 6 3,2197 3,35 3,8895 7 3,3231 3,5 5,054 TABLA 11: ERROR COMETIDO EN ALAS CON 𝚿=𝟒𝟓º EN RÉGIMEN ESTACIONARIO Por lo que se puede observar que el método desarrollado tiene un error máximo del 5% para el caso de alas con un alto valor del ángulo de flecha. Esto es debido a que el mallado para este tipo de alas debería de estar mejor adaptado y refinado. Pero hacer estos refinamientos aumenta el tiempo computacional. 81 11. CONCLUSIONES A la vista de los resultados se puede concluir que el método presentado en este trabajo supera ampliamente las expectativas iniciales. Se ha desarrollado un método que consigue una precisión de cálculo alta, que arroja valores similares a los obtenidos experimentalmente y a los que se pueden obtener con otros programas. El método desarrollado presenta una gran robustez, permitiendo el análisis de alas de cualquier geometría y en ambos regímenes: estacionario y no estacionario. También se han realizado pruebas cuando el ala se somete a distintos movimientos en pleno vuelo. Por ejemplo las siguientes situaciones físicas. - Ala arranca y pasa de estar en reposo a un ángulo de ataque 𝛼0, a moverse con una velocidad 𝑈∞. En un instante t0, el ángulo de ataque pasa a ser instantáneamente el doble. (Comportamiento ante un escalón) - Ala vuela a una velocidad 𝑈∞y tiene además un movimiento periódico en la dirección vertical Z(t). Angulo de ataque permanece nulo. Para el caso de respuesta ante un escalón en un instante determinado, en el cual el ángulo de ataque pasa de valer α=5º a un valor de α=10º, se ha obtenido la siguiente evolución de la pendiente de la curva de sustentación. FIGURA 32: RESPUESTA ANTE UN ESCALÓN Se observa que tras la perturbación se produce un descenso acusado del valor de la pendiente de la curva de sustentación. Esto es lógico debido a que el ángulo de ataque pasa a ser el doble de forma instantánea. La razón de que la pendiente de la curva de sustentación no sea exactamente la mitad del valor para el ala en flecha, es porque el coeficiente de sustentación aumenta, pero lo hace evidentemente de forma menos abrupta que el ángulo de ataque. AR=4 E=1 Ψ=0 82 Tras la perturbación se puede ver como el valor de la pendiente de la curva de sustentación vuelve a tender al valor que tendría en régimen estacionario como era de esperar (se puede comprobar que el valor límite coincide con el de la tabla, ya que se recuerda que este valor es constante en régimen estacionario según la teoría potencial linealizada. En cuanto a la evolución para el caso de movimiento oscilatorio a ángulo de ataque dado (5º) se ha obtenido la siguiente gráfica FIGURA 33: EVOLUCIÓN DEL COEFICIENTE DE SUSTENTACIÓN EN UN MOVIMIENTO OSCILATORIO Donde: Ω= 𝑤𝑐 2𝑈∞ Es la frecuencia adimensional y 𝑤 la frecuencia dimensional. Estas pruebas ya no se han podido estudiar a fondo por falta de tiempo y porque quedaban fuera del objetivo principal del trabajo. Se deja para estudios posteriores la validación de estos resultados. Las condiciones mostradas anteriormente son las que se muestran en la figura siguiente, la cual se puede encontrar en Katz y Plotkin, 1991 [2] FIGURA 34: EVOLUCIÓN DEL COEFICIENTE DE SUSTENTACIÓN EN UN MOVIMIENTO OSCILATORIO. (KYP) [2] AR=4 E=1 Ψ=0 H0=0.1 Ω=0.5 α =-5º 83 12. DESARROLLOS FUTUROS Tras este proyecto se proponen los siguientes desarrollos: - Validación del método con resultados obtenidos de la literatura para el caso de movimientos verticales Z(t) y de cabeceo 𝛼0(t) y estudio del caso de deflexión de superficies hipersustentadoras simples como flaps o slats - Implementación de las ecuaciones de la elasticidad para resolver el problema aeroelástico para alas con distinta geometría y calcular velocidades de Flutter y de divergencia. Estudio del efecto de la flecha en dichas velocidades. - Optimización de la posición de los alerones para obtener un momento de balance deseado dentro de los límites de la normativa para alas de distinta geometría. - Implementación de las ecuaciones de la elasticidad para resolver el problema aeroelástico para alas con distinta geometría y calcular velocidades de inversión de mando. 84 85 REFERENCIAS [1] José Manuel Gordillo Arias de Saavedra y Guillaume Riboux Acher, "Introducción a la Aerodinámica Potencial" Paraninfo,2012. [2] Joseph Katz y Allen Plotkin, “Low-Speed Aerodynamics: From Wing Theory to Panel Methods”, McGraw-Hill, 1991. [3] Dowell, E. H., E. F. Crawley, H. C. Curtiss, Jr., D. A. Peters, R. H. Scanlan, and F. Sisto, “A Modern Course in Aeroelasticity”, 3rd ed., Kluwer Academic Publishers, 1995. [4] Flores Caballero, Míguel Ángel, “Aplicación de la Teoría Potencial Linealizada para el estudio de las fuerzas ejercidas por una corriente incidente sobre placas rectangulares oscilantes”, Sevilla 2016 86 87 ANEXOS 1. CÓDIGO PRINCIPAL DE CÁLCULO PARA EL CASO NO ESTACIONARIO DE ALA EN FLECHA. %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% CÁLCULO DE COEFICIENTE DE SUSTENTACIÓN PARA ALAS EN FLECHA %%%%%%%%%%% %% %%%%%%%%% RÉGIMEN NO ESTACIONARIO %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% AUTOR: Fernando Moreno Pino %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% TUTOR : José Manuel Gordillo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; close all; clc %% %%% INTRODUCCIÓN DE LOS PARÁMETROS alpha0 = 5*pi/180; vecx=[-0.238619186083197; -0.661209386466265; -0.932469514203152;0.238619186083197; 0.661209386466265; 0.932469514203152]; vecw=[0.467913934572691; 0.360761573048139; 0.171324492379170; 0.467913934572691; 0.360761573048139; 0.171324492379170]; dospi=1/2/pi; for i=1:6 vecthetak(i)=0.5*pi*(vecx(i)+1); end %% GEOMETRÍA DEFINICION DE NUMERO DE PANELES Y DE INCOGNITAS AR=5; % Alargamiento S=AR; % Superficie E=1; % Estrechamiento b=sqrt(S*AR); % envergadura psi = 0*pi/180; % Angulo de flecha c = 2*sqrt(S/AR)/(1+E); % Cuerda iala = 15; % Numero de líneas a lo largo del eje x en el ala Nyincog = 4*b+1; % Numero de incógnitas a lo largo del eje y Nyincog2 = (Nyincog+1)/2; ct=c*E; % Cuerda en los extremos ifin = iala+1; % Primera division de la estela dthetax = pi/(2*(iala-1)); % División de angulo que barre las coordenadas X dthetay = pi/(Nyincog+1); % División de ángulo que barre las coordenadas Y dthetay2= pi/(Nyincog2); deltay = b/(Nyincog+1); % Separación en el mallado deltax = c/(iala-1); Nxec = iala-2; Nyec = Nyincog; 88 NTec = Nxec*Nyec; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% GEOMETRÍA. DEFINICION DE LA FORMA EN PLANTA Y MALLADO %%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% MALLADO DOBLE EN Y y(1:Nyincog2+1)=(0.25*b+0.25*b*cos(((1:Nyincog2+1)-1)*dthetay2)); y(Nyincog2+1:Nyincog+2)=(-0.25*b+0.25*b*cos(((1:Nyincog2+1)-1)*dthetay2)); xa = abs(y)'*tan(psi); % puntos del borde de ataque xs = c+(b/2*tan(psi)+ct-c)/(b/2)*abs(y)'; % borde de salida ifin = iala+1; % Final de la estela Nxec = iala-2; Nyec = Nyincog; NTec = Nxec*Nyec; % El numero total de ecuaciones donde imponer la velocidad vertical mx = zeros (Nyincog+2, iala+1); % Matriz que almacenará todas las componentes x de los puntos del mallado my = zeros (Nyincog+2, iala+1); % Matriz que almacenará todas las componentes y de los puntos del mallado UNSTEADY=1; if UNSTEADY ==0 iestela=1; finT=0; Dt=20*b/iestela; else Dt=0.25; iestela=200; % longitud de la estela finT=iestela; end conT=0; for j=1:iala my(:,j)=y'; end for j= 1:Nyincog+2 for i=1:iala mx(j,i) = xa(j)+(xs(j)-xa(j))*(1-cos((i-1)*dthetax)); end end for j=1:Nyincog+2 for i=1:iestela-1 if UNSTEADY ==0; vecy(j) = my(Nyincog+3-j, iala); 89 mx (j,iala+1)=mx(j,iala)+100*b; my(j,iala+1) = (0.5*b*cos((j-1)*dthetay)); else mx (j,iala+i)=mx(j,iala)+i*Dt; my(j,iala+i) = my(j,iala); end end end mxS=mx(:,iala:iala+1); % Primera franja de estela en x myS=my(:,iala:iala+1); % Primera franja de estela en y mxE=mx(:,iala+1:end); %matrices de la estela en x myE=my(:,iala+1:end); % matrices de la estela en y Dxbs=mx(1,iala)-mx(1,iala-1); %% PINTAR ALA %% figure (10) for i = 1:iala % Este bucle se encarga de pintar las lineas en vertical de los paneles plot( mx(:,i), my(:, i),'r') hold on end for j = 1:Nyincog+2 % Este bucle pinta las lineas horizontales de los paneles plot(mx(j,1:iala), my(j,1:iala)) hold on end %% PINTAR ESTELA %% figure(11) for i = iala: iestela+iala-1 % Este bucle se encarga de pintar las lineas en vertical de los paneles plot( mx(:,i), my(:, i),'r') hold on end for j = 1:Nyincog+2 % Este bucle pinta las lineas horizontales de los paneles plot(mx(j,iala:iestela+iala-1 ), my(j,iala:iestela+iala-1)) hold on end %%%%%%%%%%%%%%%%%%%%%%%%%%% %% CREACION DE MATRICES %%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%% sid = ((iala-2)*(Nyincog)+Nyincog); %Dimensiones de la matriz de coeficientes en ala %(en "sid" puntos se impone la condicion de contorno de %impenetrabilidad) sidE=((iala-2)*(Nyincog)); % Dimensión de la matriz de coeficientes en estela wEIm1=zeros(sid,1); % Velocidad vertical creada por estela % debido a torbellinos inyectados en la estela hasta el instante t-1 for j=1:Nyincog 96 vfila(1) = j; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i; vfila(3)=j-1; vcolumna(3)=i-1; vfila(4)=j; vcolumna(4)=i-1; end end for triangulo=1:2 x1=mx(vfila(1),vcolumna(1)); y1=my(vfila(1),vcolumna(1)); x2=mx(vfila(triangulo+1),vcolumna(triangulo+1)); y2=my(vfila(triangulo+1),vcolumna(triangulo+1)); x3=mx(vfila(triangulo+2),vcolumna(triangulo+2)); y3=my(vfila(triangulo+2),vcolumna(triangulo+2)); if i<=iala %No estoy en la estela columnaphi1=vcolumna(1)-1+(vfila(1)-2)*(iala-1); columnaphi2=vcolumna(triangulo+1)-1+(vfila(triangulo+1)-2)*(iala-1); columnaphi3=vcolumna(triangulo+2)-1+(vfila(triangulo+2)-2)*(iala-1); end D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1); c1=(y2-y3)/D; d1=(x3-x2)/D; c2=(y3-y1)/D; d2=(x1-x3)/D; c3 =(y1-y2)/D; d3 =(x2-x1)/D; fs facphi1=dospi*(vI(1)+c1*vI(3)+d1*vI(2)); facphi2=dospi*(c2*vI(3)+d2*vI(2)); facphi3=dospi*(c3*vI(3)+d3*vI(2)); if vfila(1)>1 && vfila(1)<Nyincog+2 && vcolumna(1)>1 Mat(f0,columnaphi1)=Mat(f0,columnaphi1)+facphi1; End if vfila(triangulo+1)>1 && vfila(triangulo+1)<Nyincog+2 && vcolumna(triangulo+1)>1 Mat(f0,columnaphi2)=Mat(f0,columnaphi2)+facphi2; End if vfila(triangulo+2)>1 && vfila(triangulo+2)<Nyincog+2 && vcolumna(triangulo+2)>1 Mat(f0,columnaphi3)=Mat(f0,columnaphi3)+facphi3; end end 97 3. Montaje_SOLAPE %% FUNCION el montaje en el borde de salida y comienzo de la estela x0=mx(l, k); y0=my(l,k); % Coordenadas de la singularidad f0 =k-1+(l-2)*(iala-2); Nyincog2=(Nyincog-1)/2; vfila=zeros(1,4); vcolumna=zeros(1,4); if j<Nyincog2+3 if mod(i,2)==0 tiporectangulo=1; vfila(1) = j; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i; vfila(3)=j-1; vcolumna(3)=i-1; vfila(4)=j; vcolumna(4)=i-1; else tiporectangulo=2; vfila(4) = j; vcolumna(4)=i; vfila(1)=j-1; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i-1; vfila(3)=j; vcolumna(3)=i-1; end else if mod(i,2)==0 tiporectangulo=2; vfila(4) = j; vcolumna(4)=i; vfila(1)=j-1; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i-1; vfila(3)=j; vcolumna(3)=i-1; else tiporectangulo=1; vfila(1) = j; vcolumna(1)=i; 98 vfila(2)=j-1; vcolumna(2)=i; vfila(3)=j-1; vcolumna(3)=i-1; vfila(4)=j; vcolumna(4)=i-1; end end for triangulo=1:2 x1=mxS(vfila(1),vcolumna(1)); y1=myS(vfila(1),vcolumna(1)); x2=mxS(vfila(triangulo+1),vcolumna(triangulo+1)); y2=myS(vfila(triangulo+1),vcolumna(triangulo+1)); x3=mxS(vfila(triangulo+2),vcolumna(triangulo+2)); y3=myS(vfila(triangulo+2),vcolumna(triangulo+2)); if i<=2 %No estoy en la estela columnaphi1=vcolumna(1)+(vfila(1)-2)*(2); columnaphi2=vcolumna(triangulo+1)+(vfila(triangulo+1)-2)*(2); columnaphi3=vcolumna(triangulo+2)+(vfila(triangulo+2)-2)*(2); end D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1); c1=(y2-y3)/D; d1=(x3-x2)/D; c2=(y3-y1)/D; d2=(x1-x3)/D; c3 =(y1-y2)/D; d3 =(x2-x1)/D; fs facphi1=dospi*(vI(1)+c1*vI(3)+d1*vI(2)); facphi2=dospi*(c2*vI(3)+d2*vI(2)); facphi3=dospi*(c3*vI(3)+d3*vI(2)); if vfila(1)>1 && vfila(1)<Nyincog+2 && vcolumna(1)>1 MatSOLAPE(f0,columnaphi1)=MatSOLAPE(f0,columnaphi1)+facphi1; End if vfila(triangulo+1)>1 && vfila(triangulo+1)<Nyincog+2 && vcolumna(triangulo+1)>0 MatSOLAPE(f0,columnaphi2)=MatSOLAPE(f0,columnaphi2)+facphi2; End if vfila(triangulo+2)>1 && vfila(triangulo+2)<Nyincog+2 && vcolumna(triangulo+2)>0 MatSOLAPE(f0,columnaphi3)=MatSOLAPE(f0,columnaphi3)+facphi3; end end 99 4. Montaje_Estela_2 %% FUNCION PARA COLOCAR EN LA MATRIZ DE LA ESTELA x0=mx(l, k); y0=my(l,k); % Coordenadas de la singularidad f0 =k-1+(l-2)*(iala-2); Nyincog2=(Nyincog-1)/2; vfila=zeros(1,4); vcolumna=zeros(1,4); if j<Nyincog2+3 if mod(i,2)==0 tiporectangulo=1; vfila(1) = j; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i; vfila(3)=j-1; vcolumna(3)=i-1; vfila(4)=j; vcolumna(4)=i-1; else tiporectangulo=2; vfila(4) = j; vcolumna(4)=i; vfila(1)=j-1; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i-1; vfila(3)=j; vcolumna(3)=i-1; end else if mod(i,2)==0 tiporectangulo=2; vfila(4) = j; vcolumna(4)=i; vfila(1)=j-1; vcolumna(1)=i; vfila(2)=j-1; vcolumna(2)=i-1; vfila(3)=j; vcolumna(3)=i-1; else tiporectangulo=1; vfila(1) = j; vcolumna(1)=i; 100 vfila(2)=j-1; vcolumna(2)=i; vfila(3)=j-1; vcolumna(3)=i-1; vfila(4)=j; vcolumna(4)=i-1; end end for triangulo=1:2 x1=mxE(vfila(1),vcolumna(1)); y1=myE(vfila(1),vcolumna(1)); x2=mxE(vfila(triangulo+1),vcolumna(triangulo+1)); y2=myE(vfila(triangulo+1),vcolumna(triangulo+1)); x3=mxE(vfila(triangulo+2),vcolumna(triangulo+2)); y3=myE(vfila(triangulo+2),vcolumna(triangulo+2)); if i<=iestela-1 %No estoy en la estela columnaphi1=vcolumna(1)+(vfila(1)-2)*(iestela-1); columnaphi2=vcolumna(triangulo+1)+(vfila(triangulo+1)-2)*(iestela-1); columnaphi3=vcolumna(triangulo+2)+(vfila(triangulo+2)-2)*(iestela-1); end D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1); c1=(y2-y3)/D; d1=(x3-x2)/D; c2=(y3-y1)/D; d2=(x1-x3)/D; c3 =(y1-y2)/D; d3 =(x2-x1)/D; fs facphi1=dospi*(vI(1)+c1*vI(3)+d1*vI(2)); facphi2=dospi*(c2*vI(3)+d2*vI(2)); facphi3=dospi*(c3*vI(3)+d3*vI(2)); if vfila(1)>1 && vfila(1)<Nyincog+2 && vcolumna(1)>1 MatE(f0,columnaphi1)=MatE(f0,columnaphi1)+facphi1; End if vfila(triangulo+1)>1 && vfila(triangulo+1)<Nyincog+2 && vcolumna(triangulo+1)>0 MatE(f0,columnaphi2)=MatE(f0,columnaphi2)+facphi2; end if vfila(triangulo+2)>1 && vfila(triangulo+2)<Nyincog+2 && vcolumna(triangulo+2)>0 MatE(f0,columnaphi3)=MatE(f0,columnaphi3)+facphi3; end end 101 5. Montaje_0 %% FUNCION QUE COLOCA LOS VALORES DE LA MATRIZ CUANDO ANALIZAMOS LA SINGULARIDAD x0=mx(l,k); y0=my(l,k); x1=x0; y1=y0; vfila(1)=l; vcolumna(1)=k; vfila(2)=l; vcolumna(2)=k+1; vfila(3)=l-1; vcolumna(3)=k+1; vfila(4)=l-1; vcolumna(4)=k; vfila(5)=l-1; vcolumna(5)=k-1; vfila(6)=l; vcolumna(6)=k-1; vfila(7)=l+1; vcolumna(7)=k-1; vfila(8)=l+1; vcolumna(8)=k; vfila(9)=l+1; vcolumna(9)=k+1; vfila(10)=l; vcolumna(10)=k+1; fila0=k-1+(l-2)*(iala-2); columna0=k-1+(l-2)*(iala-1); for triangulo=1:8 x2=mx(vfila(triangulo+1),vcolumna(triangulo+1)); y2=my(vfila(triangulo+1),vcolumna(triangulo+1)); x3=mx(vfila(triangulo+2),vcolumna(triangulo+2)); y3=my(vfila(triangulo+2),vcolumna(triangulo+2)); columnaphi2=vcolumna(triangulo+1)-1+(vfila(triangulo+1)-2)*(iala-1); columnaphi3=vcolumna(triangulo+2)-1+(vfila(triangulo+2)-2)*(iala-1); D=(x2-x1)*(y3-y1)-(y2-y1)*(x3-x1); c1=(y2-y3)/D; d1=(x3-x2)/D; c2=(y3-y1)/D; d2=(x1-x3)/D; c3=(y1-y2)/D; d3=(x2-x1)/D; fs_0 facphi1=dospi*(vI(1)); Mat(fila0,columna0)=Mat(fila0,columna0)+facphi1; end 102 6. FS %% PROGRAMA QUE REALIZA LA INTEGRACIÓN EN CASO DE ESTAR EN UN PUNTO QUE NO ES LA %%SINGULARIDAD %%% r1r0x=x1-x0; r1r0y=y1-y0; r2r1x=x2-x1; r2r1y=y2-y1; r3r1x=x3-x1; r3r1y=y3-y1; r3r2x=x3-x2; r3r2y=y3-y2; L12=sqrt(r2r1x*r2r1x+r2r1y*r2r1y); L10=sqrt(r1r0x*r1r0x+r1r0y*r1r0y); theta2=angulo(r2r1x,r2r1y); theta3=angulo(r3r1x,r3r1y); if(abs(theta3)<1e-12 && theta3<theta2) theta3=2*pi; end theta4=angulo(r3r2x,r3r2y); if theta2<0||theta3<0||theta4<0 disp('para') pause end vI1k=0; vI2k=0; vI3k=0; for p=1:6 vthetak=0.5*(theta3-theta2)*(vecx(p)+1)+theta2; vrk=sin(theta4-theta2)*L12/(sin(theta4)*cos(vthetak)-cos(theta4)*sin(vthetak)); vx0k=cos(vthetak)*r1r0x+sin(vthetak)*r1r0y; vxk=vrk+vx0k; va2k=L10*L10-vx0k*vx0k; S0k=1/sqrt(vxk*vxk+va2k)-1/L10; if(abs(va2k)>1e-8) S3k=(1/va2k)*(vxk/sqrt(va2k+vxk*vxk)-vx0k/sqrt(va2k+vx0k*vx0k)); S4k=log((sqrt(va2k+vxk*vxk)+vxk)/(sqrt(va2k+vx0k*vx0k)+vx0k)); S1k=-S0k-vx0k*S3k; S2k=S4k-2*vx0k*(-S0k-vx0k*S3k)-L10*L10*S3k; else S3k=-0.5*(1/((vrk+L10)*(vrk+L10))-1/(L10*L10)); S4k=log((L10+vrk)/L10); S1k=-S0k-vx0k*S3k; S2k=S4k-2*vx0k*(-S0k-vx0k*S3k)-L10*L10*S3k; 103 end vI1k=vI1k+S1k*vecw(p); vI2k=vI2k+S2k*sin(vthetak)*vecw(p); vI3k=vI3k+S2k*cos(vthetak)*vecw(p); end vI(1)=0.5*(theta3-theta2)*vI1k; vI(2)=0.5*(theta3-theta2)*vI2k; vI(3)=0.5*(theta3-theta2)*vI3k; 104 7. FS_0 %% INTEGRACION EN CASO DE SINGULARIDAD r2r1x= x2-x1; r2r1y= y2-y1; r3r1x= x3-x1; r3r1y=y3-y1; r3r2x=x3-x2; r3r2y=y3-y2; L12 = sqrt(r2r1x*r2r1x+r2r1y*r2r1y); theta2=angulo(r2r1x,r2r1y); theta3=angulo(r3r1x,r3r1y); if(abs(theta3)<1e-12 && theta3<theta2) theta3=2*pi; end theta4=angulo(r3r2x,r3r2y); vI1k=0; for p=1:6 vthetak=0.5*(theta3-theta2)*(vecx(p)+1)+theta2; vrk=sin(theta4-theta2)*L12/(sin(theta4)*cos(vthetak)-cos(theta4)*sin(vthetak)); vI1k=vI1k+(-1/vrk)*vecw(p); end vI(1)=0.5*(theta3-theta2)*vI1k; 105 8. Ángulo %% FUNCION QUE COGE EL ÁNGULO CORRECTO PARA LA INTEGRACION %% function [theta]=angulo(vx,vy) theta0=atan(abs(vy)/abs(vx)); if (abs(vx)>1e-15) theta0=atan(abs(vy)/abs(vx)); else if(vy>0) theta=pi/2.; else theta=3*pi/2.; end end if (abs(vy)<1e-15) if(vx>0) theta=0.; else theta=pi; end else if(vy>1e-15 && vx>1e-15) theta=theta0; end if(vy>1e-15 && vx<-1e-15) theta=pi-theta0; end if(vy<1e-15 && vx<-1e-15) theta=pi+theta0; end if(vy<1e-15 && vx>1e-15) theta=2*pi-theta0; end end end 112 facphi2=dospi*(c2*vI(3)+d2*vI(2)); facphi3=dospi*(c3*vI(3)+d3*vI(2)); if vfila(1)>1 && vfila(1)<Nyincog+2 && vcolumna(1)>1 Mat(f0,columnaphi1)=Mat(f0,columnaphi1)+facphi1; End if vfila(triangulo+1)>1 && vfila(triangulo+1)<Nyincog+2 && vcolumna(triangulo+1)>1 Mat(f0,columnaphi2)=Mat(f0,columnaphi2)+facphi2; End if vfila(triangulo+2)>1 && vfila(triangulo+2)<Nyincog+2 && vcolumna(triangulo+2)>1 Mat(f0,columnaphi3)=Mat(f0,columnaphi3)+facphi3; end end