Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible

Equation Chapter 1 Section 1 Trabajo Fin de Grado Grado de Ingeniería Aeroespacial Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Autor: Francisco Muñoz Soler Tutor: Miguel Pérez-Saborid Sánchez Pastor Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 iii Trabajo Fin de Grado Grado en Ingeniería Aeroespacial Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Autor: Francisco Muñoz Soler Tutor: Miguel Pérez-Saborid Sánchez Pastor Profesor titular Dep. de Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016 v Trabajo Fin de Grado: Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Autor: Francisco Muñoz Soler Tutor: Miguel Pérez-Saborid Sánchez Pastor El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2016 El Secretario del Tribunal vii A mi familia, por ese apoyo incondicional que representan. A mis amigos, por hacer el largo camino mucho más ameno. A mis maestros, en especial a mi tutor Miguel Pérez-Saborid por su implicación y apoyo en este proyecto. A todos los que hicieron esto posible, muchas gracias. ix Resumen El objetivo de este proyecto es la implementación del método Vortex-Lattice para el estudio de la Aerodinámica no estacionaria e incompresible de alas. Las ecuaciones que rigen el comportamiento de los fluidos no permiten la obtención de una solución analítica, por ello, deben ser resueltas numéricamente a partir de métodos como el expuesto en la presente memoria. La principal ventaja que presenta este método es que ofrece la posibilidad de entender la física del problema al mismo tiempo que se obtienen resultados para problemas reales, los cuales no serían abordables desde un punto de vista analítico; por esta razón, resulta idóneo para su uso didáctico en las escuelas de ingeniería. De hecho, el método Vortex-Lattice para el caso estacionario es ampliamente tratado en las clases de Aerodinámica II de la Escuela Superior de Ingeniería de Sevilla. En este proyecto, se propone extender el método desde el caso estacionario al no estacionario. Figura 4-8. Comparación geométrica ala elíptica y hexagonal 49 Figura 4-9. Formas en planta para distintos estrechamientos incluyendo ángulo de flecha 50 Figura 4-10. Coeficiente de sustentación incluyendo ángulo de flecha 50 Figura 4-11. Coeficiente de sustentación para diferentes ángulos de flecha 51 Figura 4-12. Plano en el que está contenida la componente normal de la velocidad que ve cada perfil 51 Figura 4-13. Problema de Theodorsen 53 Figura 4-14. Coeficiente de sustentación para problema de Theodorsen con diferentes valores de Λ 54 Figura 5-1. Placa plana con dos grados de libertad ℎ,𝛼 56 Figura 5-2. Placa rectangular flexible semiempotrada. Imagen adaptada de [5] 60 Figura 5-3. Respuesta para la velocidad de flameo 𝑈𝐹 del ala de alargamiento Λ=100 67 Figura 5-4. Respuesta para una velocidad superior a la de flameo 𝑈𝐹 del ala de alargamiento Λ=100 67 Figura 5-5. Respuesta para una velocidad inferior a la de flameo 𝑈𝐹 del ala de alargamiento Λ=100 68 Figura 5-6. Respuesta para 𝑈∞=0.55 para un ala de alargamiento Λ=10 69 Figura 5-7. Evolución de 𝑈𝐹 frente al alargamiento Λ 70 Figura 5-8. Evolución de 𝑈𝐹 frente al alargamiento Λ según [5] 70 Figura 5-9. Vista de perfil de la vibración de una placa flexible 71 Figura 5-10. Vista experimental de perfil de la vibración de una placa flexible 72 Figura 5-11. Resultados del método de integración para frecuencias bajas 73 Figura 5-12. Resultados del método de integración para frecuencias altas 73 Figura A-1. Segmento 𝑃𝑄 78 1 INTRODUCCIÓN a compleja maquinaria del mundo moderno no sería posible sin el conveniente aprovechamiento que se ha hecho históricamente del fluir del agua. No es casualidad que las primitivas civilizaciones florecieran a orillas de los grandes ríos: el Nilo, el Tigris, el Éufrates… y es que aprender a controlar y dirigir el curso del agua ha sido un ingrediente critico en el desarrollo de las grandes civilizaciones. Para poder prosperar, cada sociedad tuvo que desarrollar medios para manipular, controlar y distribuir las corrientes de agua. Por tanto, mucho antes de que Newton estableciera sus principios para la Mecánica Clásica, la humanidad tuvo que adquirir un conocimiento de las características fundamentales de algunos fluidos y aprender a manipularlos, es decir, tuvo la necesidad de comenzar a estudiar la ‘Mecánica de Fluidos’. Esto significa que al principio, esta disciplina fue algo plenamente experimental y fue conocida como ‘Hidráulica’ debido a su preocupación central: el agua. El avance en la comprensión del comportamiento de los fluidos fue muy lento y la ausencia de una verdadera teoría sobre el comportamiento de los fluidos y, sobre todo, de las matemáticas y ecuaciones que describieran ese comportamiento, hizo que nuestro conocimiento fuera eminentemente cualitativo. La dificultad en avanzar hacia conocimientos más precisos utilizando estos procedimientos y sin un aparato teórico más avanzado provocó finalmente un estancamiento en casi todo lo relacionado con el conocimiento de los fluidos. No fue hasta Leonardo da Vinci cuando se comenzó a atacar el problema desde un punto de vista más científico. El italiano llevó a cabo numerosos experimentos sobre el flujo de agua y aire alrededor de objetos, y documentó sus descubrimientos en detallados diagramas. L “Las matemáticas son el lenguaje en el que Dios ha escrito el Universo” - Galileo Galilei - Introducción 2 Figura 1-1. Diagrama de Turbulencia realizado por Leonardo da Vinci. En la época de Leonardo da Vinci, la Física aún no requería de un instrumento matemático para su estudio, fue con Galileo Galilei cuando dicha relación comenzó a llevarse a cabo. De este modo, una nueva teoría de fluidos surgió: ‘La Hidrodinámica’. Fueron dos de los discípulos de Galileo, Benedetto Castelli y Evangelista Torricelli, los primeros en establecer las bases de esta nueva disciplina que resultaba la contrapartida teórica de la hidráulica. El problema principal de esta nueva teoría radicaba en la complejidad del comportamiento de los fluidos, lo cual supuso que la hidrodinámica solo fuese útil en casos muy concretos, fuera de ellos las predicciones no guardaban relación alguna con el comportamiento real. La limitación se encontraba en el instrumento matemático del que se disponía, hizo falta el desarrollo del cálculo infinitesimal para poder llevar a cabo una adecuada descripción del comportamiento de los fluidos que sería ideada paralelamente por Isaac Newton y Gottfried Leibniz. Gracias a la poderosa herramienta del cálculo infinitesimal pudo producirse el mayor salto en nuestra comprensión de la Mecánica de Fluidos, llevado a cabo por Leonhard Euler: el desarrollo de unas ecuaciones en derivadas parciales que se creía que permitirían describir y predecir de forma teórica el comportamiento general de cualquier fluido. Sin embargo, estas ecuaciones sólo funcionaron en algunos casos llevando en ocasiones a contradicciones tan conocidas como la Paradoja de D’Alembert. Por esta razón, la hidrodinámica quedó relegada al puesto de mera curiosidad teórica ya que los ingenieros seguían obteniendo resultados mucho mejores acudiendo a métodos empíricos antes que a las ecuaciones de Euler. Se había perdido una batalla, pero la hidrodinámica tenía aún mucha guerra por dar y es que en el siglo XIX, acudieron al rescate el inglés Sir George Stokes y el francés Claude-Louis Navier, quienes establecieron en 1822 unas nuevas ecuaciones que describían adecuadamente el comportamiento de los fluidos. Desde entonces, el estudio de la Mecánica de Fluidos no volvió a ser igual, el estudio teórico predecía correctamente las observaciones experimentales, pasando de ser una mera curiosidad a alcanzar el puesto que ocupa actualmente: una poderosa herramienta sin la cual no seríamos capaces de manejar la compleja maquinaria de la que se habló al principio. La diferencia entre ‘hidráulica’ e ‘hidrodinámica’ comenzó a desaparecer y se formó lo que se conoce como ‘Mecánica de Fluidos’. Es decir, desde mediados del siglo XIX los ingenieros comienzan a utilizar más las ecuaciones diferenciales y es en el siglo XX donde se encuentra un nuevo obstáculo, las matemáticas nos daban una descripción muy buena de la realidad, pero en muchos de los casos esta descripción presentaba un comportamiento caótico, es decir, prácticamente imposible de calcular con exactitud para tiempos relativamente alejados del actual. Las matemáticas nos ofrecían el modo en que funcionaban los fluidos, pero nuestra capacidad de cálculo no alcanzaba a hacer frente a la complejidad de estas ecuaciones diferenciales. 3 3 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Si se ha comentado que en el siglo XIX el rescate fue llevado a cabo por Navier y Stokes, en nuestro siglo el rescate lo ha realizado la informática, y es que sin la potencia de cálculo que nos ofrecen los programas informáticos, tendríamos las ecuaciones que gobiernan a los fluidos, pero seríamos incapaces de resolverlas para dar lugar a predicciones precisas del comportamiento real. Dentro de la Mecánica de Fluidos se encuentra la rama de la ‘Aerodinámica’, que resulta de especial interés para los ingenieros aeroespaciales, ya que nos permite comprender las interacciones que el flujo del aire realiza sobre la geometría de una aeronave, permitiendo que el vuelo controlado sea posible. Sin embargo, las ecuaciones que describen correctamente el flujo del aire son como ya se ha dicho anteriormente especialmente difíciles de resolver, es por ello, que es necesaria la implementación de un método numérico para la obtención de predicciones precisas de la realidad. La dificultad en la resolución de dichas ecuaciones representa un problema no sólo en el ámbito ingenieril, sino también en el didáctico, puesto que a la hora de enseñar esta disciplina los problemas abordables en clase suelen ser demasiado artificiales buscando el poder obtener una solución analítica. Esta artificialidad da como resultado pocas aplicaciones prácticas de las soluciones obtenidas, sin embargo, con el desarrollo informático experimentado en los últimos años, es posible abordar situaciones más reales mediante el uso de métodos numéricos implementados en clase, lo cual ofrece al alumno la posibilidad de entender la física del problema y al mismo tiempo obtener soluciones que concuerdan con bastante aproximación con los hechos experimentales, y todo esto sin necesidad de una excesiva manipulación analítica de las expresiones algebraicas que rigen el comportamiento de los fluidos. Por tanto, el objetivo de este Trabajo Fin de Grado es doble: por un lado demostrar que la potencia de los métodos numéricos tienen hoy en día una importancia vital en el ámbito de la ingeniería y por otro lado, que dicha potencia puede ser ampliamente aprovechada para enseñar a los futuros ingenieros cual es la realidad del comportamiento de los fluidos. Para ello, se va a extrapolar el método Vortex-Lattice desarrollado en las clases de Aerodinámica II de la Escuela Superior de Ingeniería de Sevilla desde el caso estacionario al no estacionario e incompresible. 1.1 Aerodinámica no estacionaria Dentro de la disciplina de la Aerodinámica, hay que distinguir una rama de gran interés práctico: la Aerodinámica Potencial. La Aerodinámica potencial es la parte de la Aerodinámica que estudia el flujo de gases alrededor de objetos fuselados despreciando los esfuerzos de viscosidad, ya que en los casos para los que se aplica, éstos se encuentran confinados en regiones muy estrechas del dominio fluido. En las escuelas de Ingeniería se suele comenzar estudiando la Aerodinámica potencial para el caso estacionario, esto es así porque se pueden eliminar términos de variación con el tiempo en las ecuaciones de Navier-Stokes lo cual simplifica enormemente el estudio y compresión de los fenómenos aerodinámicos. Sin embargo, el tener en cuenta los términos de variación con el tiempo nos ofrece la posibilidad de obtener resultados de gran interés para la ingeniería. Gracias a un estudio de la Aerodinámica no estacionaria, se puede dar explicación a fenómenos como el vuelo de insectos, el funcionamiento de las cuerdas vocales e incluso prevenir efectos perjudiciales para las estructuras que se ven sometidas a cargas ocasionadas por ráfagas de viento, como puede ocurrirle a puentes, edificios y por supuesto, a aeronaves. Poniendo especial interés en la aplicación al campo de la Aeronáutica, citar que el estudio de la Aerodinámica no estacionaria nos proporciona la posibilidad de comprender y advertir los fenómenos del flameo y la divergencia a los que pueden verse sometidos los elementos estructurales de un aeronave. El fenómeno del flameo aparece cuando un sistema mecánico comienza a oscilar sin presentar amortiguamiento en respuesta a una corriente fluida. En general, a mayor velocidad de la corriente incidente, más tarda en amortiguarse la respuesta, hasta que se alcanza un punto tal que ésta no llega a amortiguarse nunca y se mantiene en una oscilación armónica. En esta situación, el sistema mecánico está extrayendo energía del fluido, la cantidad de energía extraída es igual a la que se disipa por el amortiguamiento del sistema. Este fenómeno se puede observar en una bandera ondeando (por ello, un sinónimo de ondear es flamear), además, también puede aparecer en determinados momentos para el ala de un aeronave o incluso en puentes en los que no se haya tenido especial cuidado en evitar que puedan alcanzarse las frecuencias de Introducción 4 resonancia, como ocurrió en el puente de Tacoma-Narrows en 1940. Figura 1-2. Puente de Tacoma-Narrows flameando. Una de las aplicaciones más importantes que tendrá el método Vortex-Lattice que se busca implementar, será la predicción de las velocidades de flameo con el objetivo de prevenir la aparición de dichos fenómenos tan desfavorables para la estructura. 1.2 Objetivos y estructura del trabajo Por todo lo explicado anteriormente, los objetivos propuestos para el presente TFG son los siguientes:  Exponer la base teórica sobre la que se fundamenta el método Vortex-Lattice que se busca implementar para el caso particular de Aerodinámica no estacionaria e incompresible de un ala.  Extender con la ayuda del método Katz-Plotkin explicado en la referencia [9] y [3] el programa implementado en las clases de Aerodinámica II desde el caso estacionario al no estacionario.  Obtener con la ayuda del programa resultados para alas con diferentes geometrías realizando comparaciones entre ellas y el caso 2D ya estudiado en [3] incluyendo una interpretación física de los resultados.  Exponer una aplicación del método Vortex-Lattice sobre alas en régimen no estacionario, en este caso, será el flameo de una placa ante una corriente incidente. Con estos objetivos por cumplir, el trabajo se estructurará en cinco capítulos más además del presente y un anexo. En el siguiente capítulo se van a exponer los fundamentos teóricos sobre los que se basa el método VortexLattice a implementar. Seguidamente, se explicará en qué consta dicho método en el capítulo 3, además, se incluirá el código MatLab realizado para la implementación. En el capítulo 4, se analizarán los resultados a los que puede llegarse mediante la ejecución del programa MatLab desarrollado en el capítulo 3 y tras esto, se estudiará el flameo de un ala finita como aplicación directa del método en el capítulo 5. Finalmente, se cerrará la memoria con un último capítulo que incluirá una serie de conclusiones finales. 2 FUNDAMENTOS TEÓRICOS l objetivo del presente capítulo es el de exponer los fundamentos teóricos necesarios para el desarrollo del método Vortex-Lattice. Con el método se busca obtener las interacciones que se producen entre un ala y el fluido en el que se halla inmersa cuando ambos se mueven a velocidades diferentes. Esto es, la fuerza 𝑭𝑓𝑠y el momento 𝑴𝑓𝑠 que el fluido ejerce sobre el ala. La integración sobre la superficie del ala de los esfuerzos generados por el fluido nos darán las expresiones para dichas interacciones, dichos esfuerzos se pueden obtener con gran precisión gracias a las ecuaciones de Navier-Stokes, lo que hace posible expresar en forma de ecuaciones uno de los principales objetivos de nuestro método: el cálculo aproximado de las siguientes integrales de superficie 𝑭𝑓𝑠=∫∫ (𝑝−𝑝∞)(−𝒏𝑠)𝑑𝜎 Σ𝑎𝑙𝑎 +∫∫ 𝒏𝑠⋅𝝉′𝑑𝜎 Σ𝑎𝑙𝑎 𝑴𝑓𝑠=∫∫ (𝒙−𝒙0)×(𝑝−𝑝∞)(−𝒏𝑠)𝑑𝜎 Σ𝑎𝑙𝑎 +∫∫ (𝒙−𝒙0)×(𝒏𝑠⋅𝝉′)𝑑𝜎 Σ𝑎𝑙𝑎 (2–1) donde 𝑝 es la presión, 𝑝∞ es la presión de referencia, 𝒏𝑠 es el vector normal exterior al sólido, 𝝉′=2𝜇𝜸+𝑰(𝜇𝑣−2 3𝜇)∇⋅𝐯 (2–2) es el tensor de esfuerzos que, para un fluido Newtoniano está dado por la ley de Navier-Stokes y 𝒙0 es el vector posición del punto respecto al que tomamos momentos. No entraremos en más detalle sobre el tensor de esfuerzos puesto que será despreciado para nuestras aproximaciones y no será objeto de estudio, para una explicación detallada de éste, se puede acudir a la referencia [2]. Una vez hayamos calculado la fuerza que el fluido ejerce sobre el sólido, podemos descomponer dicha fuerza en dos componentes: sustentación y resistencia aerodinámica. Para ello, vamos a hacer previamente la elección del sistema de coordenadas. El sistema de coordenadas que tomemos será de coordenadas cartesianas y lo consideraremos inercial, dado por las coordenadas (x, y, z). Su origen se encontrará en el borde de ataque del perfil central, la dirección x estará dada según la dirección de la corriente incidente, la dirección y según la dirección de la envergadura y la dirección z completará el triedro a derechas. Este sistema de referencia viaja con el ala a velocidad constante 𝑈∞,pero el ala no se encuentra E “La verdad es demasiado complicada como para permitir nada más allá de meras aproximaciones” - John Von Neumann - Fundamentos teóricos 6 fijada a él, puesto que el estudio de los fenómenos no estacionarios se realizará mediante cambios en las posiciones que ocupa la superficie alar con respecto a este sistema de referencia. Figura 2-1. Ejes cartesianos elegidos sobre la superficie del ala. Estos ejes cartesianos tienen de base ortonormal los vectores (𝒙 ,𝒚 ,𝒛) en las direcciones x, y, z respectivamente. De este modo, puede definirse la sustentación y la resistencia aerodinámica como: 𝐷=𝑭𝑓𝑠⋅𝒙  ; 𝐿=𝑭𝑓𝑠⋅𝒛 (2–3) Los valores L y D suelen calcularse a través de los coeficientes adimensionales de la resistencia y la sustentación aerodinámica que se definen, respectivamente, como 𝐶𝐷=𝐷 1 2𝜌∞𝑈∞ 2𝐴 ; 𝐶𝐿=𝐿 1 2𝜌∞𝑈∞ 2𝐴 (2–4) donde 𝜌∞ es la densidad del aire lejos del ala, 𝑈∞ es la velocidad relativa de la aeronave con respecto al medio en el que se desplaza y 𝐴 representa una de las áreas características, para nuestro caso particular de un ala, suele ser el área total de dicha ala. Dicho esto, estamos en condiciones de abordar nuestro problema:  En primer lugar, presentaremos las ecuaciones generales de Navier-Stokes en el apartado 2.1.  En el apartado 2.2, simplificaremos estas ecuaciones para hacer posible la resolución del problema.  Posteriormente, introduciremos algunas definiciones que nos llevarán al Teorema de BjerknessKelvin, lo cual implica necesariamente tres conclusiones de vital importancia que serán analizadas detalladamente.  En el apartado 2.4, se introducirá la hipótesis de pequeñas perturbaciones sobre el flujo de la corriente incidente por parte del ala y se sacarán las conclusiones de tomar esta suposición.  En el apartado 2.5, se simplificarán las condiciones de contorno introducidas en el apartado 2.1 y se impondrán nuevas que nos ayuden a cerrar el sistema de ecuaciones del problema simplificado.  Finalmente, se hará una descomposición del problema para destacar que el objetivo de esta memoria es centrarse principalmente en los fenómenos no estacionarios. x y z 𝑼∞ 7 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 2.1 Ecuaciones de Navier-Stokes Para poder realizar las integrales 2-1 y calcular así los coeficientes 𝐶𝐿 y 𝐶𝐷 hay que conocer previamente los valores de la presión y del tensor de esfuerzos sobre la superficie del ala. Para ello, tenemos que resolver las ecuaciones de Navier-Stokes demostradas en la referencia [2] sujetas a unas determinadas condiciones de contorno. El problema completo a resolver podemos encontrarlo enunciado en 2-5, 2-6 y 2-7 a partir de [8]: Las ecuaciones de continuidad, cantidad de movimiento y de la energía que forman el sistema no lineal de ecuaciones en derivadas parciales de Navier-Stokes: 𝜕𝜌 𝜕𝑡+∇⋅(𝜌𝒗)=0 𝜌𝜕𝒗 𝜕𝑡+𝜌𝒗⋅∇𝒗=−∇𝑝+∇⋅𝝉′+𝜌𝒇𝑚 𝜌𝑐𝑣𝜕𝑇 𝜕𝑡+𝜌𝑐𝑣𝒗⋅∇𝑇=−𝑝∇⋅𝒗+𝝉′:∇𝒗+𝑄𝑟+𝑄𝑞+∇⋅(𝑘∇𝑇) (2–5) donde 𝜌,𝑝,𝑇,𝒗 son los campos de densidad, presión, temperatura y velocidad respectivamente y 𝑡 denota el tiempo, 𝒇𝑚 es el vector de fuerzas másicas, 𝑐𝑣 es la capacidad calorífica del gas a volumen constante, 𝑘 es su conductividad térmica, 𝝉′:∇𝒗>0 es el término de disipación de energía cinética en energía interna y 𝑄𝑟 y 𝑄𝑞 son las potencias caloríficas que, por unidad de volumen, recibe el fluido por radiación y por reacción química respectivamente. S 𝒙∈Σ𝑎𝑙𝑎: 𝒗=0 𝑇=𝑇𝑎𝑙𝑎 𝑘𝜕𝑇 𝜕𝑛𝑎𝑙𝑎1=𝑘𝑎𝑙𝑎𝜕𝑇𝑎𝑙𝑎 𝜕𝑛𝑎𝑙𝑎 (2–6) Si 𝒙→∞: 𝒗=𝑈∞𝒙  𝑝−𝑝∞=0 𝑇−𝑇∞=0 (2–7) Este sistema de cinco ecuaciones y seis incógnitas (𝒗,𝜌,𝑝,𝑇) ha de ser complementado con la ecuación de estado de los gases perfectos: 𝑝 𝜌=𝑅𝑔𝑇 (2–8) Puede parecer que necesitemos una ecuación adicional para calcular 𝑇𝑎𝑙𝑎, sin embargo, el cálculo aproximado de las integrales en 2-1 es independiente del flujo de calor a través de las paredes del sólido, tal y como se demuestra en [8]. Nos queda aún por imponer condiciones iniciales a nuestro problema, éstas las impondremos para cada problema particular que resolveremos una vez hayamos implementado el método numérico para la resolución del problema. 1 𝑛𝑎𝑙𝑎 es el vector normal exterior a la superficie del ala. Fundamentos teóricos 8 2.2 Ecuaciones simplificadas Para comenzar, conviene recordar que nuestro problema objeto de estudio se iba a restringir al caso incompresible, el cual impone que Δ𝜌/𝜌 ≪1. Si partimos de la ecuación de estado 2-8, podemos llevar a cabo el siguiente desarrollo: ln(𝑝 𝜌)=ln(𝑅𝑔𝑇) ln𝑝−ln𝜌=ln𝑅𝑔+ln𝑇 𝑑(ln𝑝−ln𝜌)=𝑑(ln𝑅𝑔+ln𝑇) 𝑑𝑝 𝑝−𝑑𝜌 𝜌=𝑑𝑇 𝑇 Δ𝑝 𝑝−Δ𝜌 𝜌~Δ𝑇 𝑇 Si suponemos ahora flujo incompresible Δ𝜌/𝜌 ≪1, para Δ𝑇/𝑇 ≪1: Δ𝜌 𝜌~Δ𝑝 𝑝~𝜌𝑈∞ 2 𝑝=𝛾𝜌𝑈∞ 2 𝛾𝑝 =𝛾𝑈∞ 2 𝑎∞ 2=𝛾𝑀∞ 2~𝑀∞ 2 Por tanto, la condición para que el flujo sea incompresible puede imponerse como 𝑀∞ 2≪1. Una vez hemos supuesto que el flujo es incompresible, la ecuación de continuidad puede simplificarse: 𝜕𝜌 𝜕𝑡+∇⋅(𝜌𝒗)=0 → ∇⋅𝒗=0 Con esta nueva ecuación, el término de la viscosidad que aparecía en la ecuación de cantidad de movimiento del sistema de Navier-Stokes también puede simplificarse: ∇⋅𝝉′=∇⋅(2𝜇𝜸+𝑰(𝜇𝑣−2 3𝜇)∇⋅𝐯)=𝜇∇2𝒗 Dejando a un lado por ahora la ecuación de la energía, y suponiendo que tanto 𝜇 como 𝜌 son constantes, el sistema de ecuaciones que nos queda es el siguiente: ∇⋅𝒗=0 𝜌(𝜕𝒗 𝜕𝑡+𝒗⋅∇𝒗)=−∇𝑝+𝜇∇2𝒗+𝜌𝒇𝑚 (2–9) El primer paso para simplificar estas ecuaciones es adimensionarlas definiendo previamente algunas magnitudes características:  𝑐𝑟: Longitud característica, cuerda en el perfil central del ala.  𝑈∞: Velocidad característica.  𝑇𝑐: Tiempo característico.  𝑝0: Presión característica.  𝑓0: Fuerza másica característica. De este modo, pueden definirse las siguientes variables adimensionales: 9 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 𝑥∗=𝑥 𝑐𝑟, 𝑦∗=𝑦 𝑐𝑟, 𝑧∗=𝑧 𝑐𝑟 𝒗∗=𝒗 𝑈∞ 𝑡∗=𝑡 𝑇𝑐 𝑝∗=𝑝 𝑝0 𝒇∗=𝒇𝒎 𝑓0 (2–10) Introduciendo (2-10) en (2-9), nos queda: ∇⋅𝒗∗=0 𝑐𝑟 𝑈∞𝑇𝑐𝜕𝒗∗ 𝜕𝑡 +𝒗∗⋅∇𝒗∗=−(𝑝0 𝜌𝑈∞ 2)∇𝑝∗+( 𝜇 𝜌𝑈∞𝑐𝑟)∇2𝒗∗+(𝑐𝑟𝑓0 𝑈∞ 2)𝒇∗ (2–11) Definimos los siguientes números adimensionales: 𝑆𝑡=𝑐𝑟 𝑈∞𝑇𝑐 𝐹𝑟=𝑈∞ √𝑐𝑟𝑓0 𝐸𝑢=𝑝0 𝜌𝑈∞ 2 𝑅𝑒=𝜌𝑈∞ 2𝑐𝑟 𝜇 (2–12) donde:  𝑆𝑡 es el número de Strouhal que marca cuando un proceso es estacionario o cuasi-estacionario, para ello se compara la variación temporal de la velocidad con el término convectivo 2 , 𝑂(𝜌𝜕𝒗 𝜕𝑡) 𝑂(𝜌𝒗⋅∇𝒗)~𝜌𝑈∞ 𝑇𝑐 𝜌𝑈∞𝑈∞ 𝑐𝑟=𝑐𝑟 𝑈∞𝑇𝑐=𝑆𝑡  𝐹𝑟 es el número de Froude que aparece cuando se plantea la importancia de las fuerzas másicas frente a las convectivas, 𝑂(𝜌𝒗⋅∇𝒗) 𝑂(𝜌𝒇𝑚)~𝜌𝑈∞𝑈∞ 𝑐𝑟 𝜌𝑓0=𝑈∞ 2 𝑐𝑟𝑓0=𝐹𝑟2  𝐸𝑢 es el número de Euler y representa la relación entre las fuerzas de presión y el término convectivo: 𝑂(∇𝑝) 𝑂(𝜌𝒗⋅∇𝒗)~𝑝0 𝑐𝑟 𝜌𝑈∞𝑈∞ 𝑐𝑟=𝑝0 𝜌𝑈∞ 2=𝐸𝑢 2 El término convectivo de la ecuación de cantidad de movimiento es 𝜌𝒗⋅∇𝒗 y representa la variación de la velocidad debido al cambio en la posición que ocupa una partícula fluida. Fundamentos teóricos 16 potencial de perturbación: 𝜕𝜙 𝜕𝑡+((𝑈∞+𝑣𝑥′)2+𝑣𝑦′2+𝑣𝑧′2) 2+𝑝 𝜌∞=𝜕𝜙∞ 𝜕𝑡 +𝑈∞2 2+𝑝∞ 𝜌∞ teniendo en cuenta que (𝑈∞+𝑣𝑥′)2+𝑣𝑦′2+𝑣𝑧′2≅𝑈∞ 2+2𝑈∞𝑣𝑥: 𝜕𝜙 𝜕𝑡+𝑈∞ 2+2𝑈∞𝑣𝑥 2+𝑝 𝜌∞=𝜕𝜙∞ 𝜕𝑡 +𝑈∞2 2+𝑝∞ 𝜌∞ 𝑝′=−𝜌∞(𝜕 𝜕𝑡+𝑈∞𝜕 𝜕𝑥)𝜙′ (2–27) 2.5 Condiciones de Contorno En el apartado 2.1 se adelantaron las condiciones de contorno del problema general para resolver las ecuaciones de Navier-Stokes, sin embargo, hemos realizado diferentes simplificaciones en estas ecuaciones que requieren una reformulación de las condiciones de contorno 2-6 y 2-7 ya formuladas. Para nuestro problema, supondremos un ala dada por las superficies 𝑧=𝑧𝑒(𝑋,𝑌,𝑡) y 𝑧=𝑧𝑖(𝑋,𝑌,𝑡) para extradós e intradós respectivamente, con (𝑋,𝑌) siendo las coordenadas de los ejes inerciales. Como se concluyó en el apartado 2.3, en este ala por ser finita, aparece una estela de torbellinos donde existe un gradiente de velocidades muy grande causado por los efectos viscosos, esta zona será supuesta como una superficie de discontinuidad de espesor nulo. Además, se hará la hipótesis de que el ángulo entre la estela y el ala es muy pequeño. A la hora de imponer las condiciones de contorno, vamos a evitar las impuestas en temperatura, puesto que no es necesario resolver el problema térmico para calcular 2-15. De este modo, las condiciones de contorno impuestas son: 1. Infinito no perturbado: Como se ha comentado que las condiciones impuestas en temperatura no serán tenidas en cuenta, la condición de contorno 2-5, se reduce a imponer condiciones para presiones y velocidades para 𝒙→∞, sin embargo, la expresión 2-27 nos ofrece una relación entre la perturbación en presiones y la perturbación en velocidades, y de ella se extrae que si no hay perturbación en velocidad, tampoco existirá perturbación en presión, por lo que, la condición 2-5 puede reducir a una condición de contorno en velocidades para 𝒙→∞: 𝒙→∞, 𝒗→𝑼∞, 𝒗′→0 (2–28) 2. Condición para velocidades verticales en el ala: Al igual que ocurre con 2-5, la ecuación 2-6 se reduce también a imponer una condición en velocidades para 𝒙∈Σ𝑎𝑙𝑎. Como hemos despreciado los efectos viscosos, no habrá nada que frene a la velocidad hasta hacerla cero en la pared, esto ocurrirá dentro de la capa límite que no es objeto de estudio en nuestro problema. Por tanto, lo que se impondrá para 𝒙∈Σ𝑎𝑙𝑎 será que ya que el fluido se desplaza con la superficie del ala, imponemos que el ala sea una superficie fluida lo cual implica, como se verá a continuación, una condición para las velocidades verticales del fluido en el ala. Sea la función 𝐹𝑒,𝑖, que igualada a 0 nos da una expresión implícita de la superficie alar: 𝐹𝑒,𝑖(𝑥,𝑦,𝑧,𝑡)=𝑧−𝑧𝑒,𝑖(𝑥,𝑦,𝑡)=0 (2–29) siendo 𝐹𝑒,𝑖(𝑥,𝑦,𝑧,𝑡)=0 una superficie fluida, se cumple que los puntos (𝑥+𝑑𝑥,𝑦+𝑑𝑦,𝑧+𝑑𝑧) 17 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible pertenecen a la mima superficie en 𝑡+𝑑𝑡, es decir, 𝐹𝑒,𝑖(𝑥,𝑦,𝑧,𝑡)=𝐹𝑒,𝑖(𝑥+𝑑𝑥,𝑦+𝑑𝑦,𝑧+𝑑𝑧,𝑡). Por lo que: 𝑑𝐹𝑒,𝑖=𝐹𝑒,𝑖(𝑥+𝑑𝑥,𝑦+𝑑𝑦,𝑧+𝑑𝑧,𝑡)−𝐹𝑒,𝑖(𝑥,𝑦,𝑧,𝑡)=0 (2–30) Desarrollando 2-30 𝑑𝐹𝑒,𝑖=𝜕𝐹 𝜕𝑥 𝑑𝑥+𝜕𝐹 𝜕𝑦 𝑑𝑦+𝜕𝐹 𝜕𝑧 𝑑𝑧+𝜕𝐹 𝜕𝑡 𝑑𝑡=0 (2–31) La expresión 2-32 es análoga a 2-33: 𝜕𝐹 𝜕𝑥𝑑𝑥 𝑑𝑡+𝜕𝐹 𝜕𝑦𝑑𝑦 𝑑𝑡+𝜕𝐹 𝜕𝑧𝑑𝑧 𝑑𝑡+𝜕𝐹 𝜕𝑡 =0 (2–32) Sabiendo que: 𝒗=𝑑𝑥 𝑑𝑡𝒙 +𝑑𝑦 𝑑𝑡𝒚 +𝑑𝑧 𝑑𝑡𝒛=𝑣𝑥𝒙 +𝑣𝑦𝒚 +𝑣𝑧𝒛 (2–33) Llegamos a: 𝜕𝐹 𝜕𝑥𝑣𝑥+𝜕𝐹 𝜕𝑦𝑣𝑦+𝜕𝐹 𝜕𝑧𝑣𝑧+𝜕𝐹 𝜕𝑡 =0 (2–34) Incluyendo el cálculo de las derivadas para la superficie fluida −𝜕𝑧𝑒,𝑖 𝜕𝑥 𝑣𝑥−𝜕𝑧𝑒,𝑖 𝜕𝑦 𝑣𝑦+𝑣𝑧+𝜕𝑧𝑒,𝑖 𝜕𝑡 =0 (2–35) Y desarrollando las componentes de la velocidad 𝒗=(𝑈∞+𝑣𝑥′)𝒙 +𝑣𝑦′𝒚 +𝑣𝑧′𝒛: −𝜕𝑧𝑒,𝑖 𝜕𝑥 𝑈∞−𝜕𝑧𝑒,𝑖 𝜕𝑥 𝑣𝑥′−𝜕𝑧𝑒,𝑖 𝜕𝑦 𝑣𝑦′+𝑣𝑧′+𝜕𝑧𝑒,𝑖 𝜕𝑡 =0 (2–36) Finalmente, despreciando los términos de orden superior: 𝑣𝑧′(𝑥,𝑦,𝑧=𝑧𝑒,𝑖)=𝜕𝑧𝑒,𝑖(𝑥,𝑦,𝑡) 𝜕𝑡 +𝑈∞𝜕𝑧𝑒,𝑖(𝑥,𝑦,𝑡) 𝜕𝑥 (2–37) Como el grosor de los perfiles del ala es muy pequeño en relación a la longitud de cuerda que tienen, se puede evaluar la condición 2-38 en 𝑧=0± en vez de en la superficie real del ala, con lo que finalmente, la condición para las velocidades verticales nos queda: 𝑣𝑧′(𝑥,𝑦,𝑧=0±,𝑡)=𝜕𝑧𝑒,𝑖(𝑥,𝑦,𝑡) 𝜕𝑡 +𝑈∞𝜕𝑧𝑒,𝑖(𝑥,𝑦,𝑡) 𝜕𝑥 (2–38) Hasta ahora, hemos simplificado las condiciones que se habían impuesto en el apartado 2.1, sin embargo, dado que se ha incluido una nueva región de discontinuidades –la estela–, ahora hay que imponer condiciones de contorno en esta nueva región: 3. Continuidad en la estela: La conservación de la masa (ecuación de continuidad) debe cumplirse a través de la estela, esta condición Fundamentos teóricos 18 queda impuesta tomando un volumen de control en torno a la estela e igualando el gasto que sale al que entra. Recordemos que el gasto se definía como 𝐺=𝜌𝑣𝐴, siendo 𝜌 la densidad del fluido, 𝑣 la velocidad y 𝐴 la superficie a través de la cual quiere calcularse dicho gasto. La aplicación de esta condición lleva a: 𝜌∞𝑣𝑧′(𝑥,𝑦,0−)𝐴=𝜌∞𝑣𝑧′(𝑥,𝑦,0+)𝐴 𝑣𝑧′(𝑥,𝑦,𝑧=0−)=𝑣𝑧′(𝑥,𝑦,𝑧=0+) (2–39) Conviene notar que la condición 2-40 se impone en 𝑧=0 ya que hemos supuesto que el ángulo entre el ala y la estela era prácticamente despreciable. 4. Igualdad de presiones en la estela: Dado que la estela es una superficie de discontinuidad que carece de masa, hay que imponer la igualdad de presiones entre extradós e intradós de la estela 3 . 𝑝𝑒(𝑥,𝑦,𝑧=0+)=𝑝𝑖(𝑥,𝑦,𝑧=0−) A partir de 2-27, se llega a: 𝜌∞(𝜕 𝜕𝑡+𝑈∞𝜕 𝜕𝑥)(𝜙𝑒′−𝜙𝑖′)=0 (2–40) Esta expresión se conoce como la condición de Kutta generalizada, puesto que para el caso estacionario (𝜕/𝜕𝑡=0) puede observarse que queda su expresión más conocida: 𝑣𝑥𝑒 ′=𝑣𝑥𝑖′. Para ver las implicaciones que tiene esta expresión, conviene ponerla en términos de la densidad de circulación. Para ello vamos a realizar un pequeño desarrollo matemático: Sea la curva R formada por los puntos ABCD de la figura 2-4 lo suficientemente cercana a uno de los perfiles del ala, de forma que 𝜙(𝐵)=𝜙𝑒 y 𝜙(𝐶)=𝜙𝑖. Se cumple que si calculamos la circulación Γ(𝑡,𝑥,𝑦=𝑦0) para 𝑦=𝑦0 𝑐𝑡𝑒 análogamente a como se hizo en la expresión 2-24: Γ(𝑡,𝑥,𝑦=𝑦0)=∮𝒗⋅𝑑𝒍=∮∇𝜙⋅𝑑𝒍=∮𝑑𝜙 𝑅= 𝑅𝑅 𝜙(𝐵)−𝜙(𝐶)=𝜙𝑒−𝜙𝑖=𝜙𝑒′−𝜙𝑖′ (2–41) 3 Prestar especial atención a que esta condición se impone en la estela, no en el ala. En el ala, esta diferencia de presiones es distinta de cero y es precisamente la que genera sustentación. 19 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Figura 2-4. Curva ABCD sobre la que calculamos circulación Γ(𝑡,𝑥,𝑦) Si introducimos 2-42 en 2-41: 𝜌∞(𝜕 𝜕𝑡+𝑈∞𝜕 𝜕𝑥)Γ(𝑡,𝑥,𝑦)=0 (2–42) La expresión 2-41 implica que si un observador viaja con la corriente incidente 𝑈∞, siempre ve la misma circulación Γ, es decir, el valor de la circulación Γ en 𝑥 para un instante t se traslada hasta el punto 𝑥+𝑈∞𝑑𝑡 en un tiempo 𝑑𝑡. Lo cual puede comprobarse sin más que observar que la variación que se experimentaría en la circulación en ese incremento de tiempo 𝑑𝑡 es cero: Γ(𝑡+𝑑𝑡,𝑥+𝑈∞𝑑𝑡,𝑦)−Γ(𝑡,𝑥,𝑦) 𝑑𝑡 =𝜕Γ 𝜕𝑡+𝑈∞𝜕Γ 𝜕𝑥=0 (2–43) Este hecho será de suma importancia a la hora de implementar en el método Vortex-Lattice la evolución de la circulación en la estela, puesto que hemos comprobado que la circulación calculada en 𝑥 para un instante t será la misma que tendremos en 𝑥+𝑈∞𝑑𝑡 para el instante 𝑡+𝑑𝑡. 5. Teorema de Bjerkness-Kelvin Una vez hemos impuesto todas las condiciones de contorno vistas, sólo queda garantizar el cumplimiento del teorema de Bjerkness-Kelvin por parte de nuestra solución. 2.6 Descomposición del problema Si se observan la ecuación del potencial 2-25 y las condiciones de contorno, puede comprobarse que son lineales, por lo que es posible aplicar el principio de superposición y separar el problema en dos: uno estacionario y otro no estacionario, entre los que la única diferencia está en la forma de imponer la condición de impenetrabilidad. El hecho de que la única diferencia esté en la forma de imponer la condición de impenetrabilidad, se debe a que las ecuaciones que describen la geometría del ala pueden descomponerse en dos componentes, una estacionaria (E) y otra no estacionaria (N): 𝑧𝑒,𝑖(𝑡,𝑥,𝑦)=𝑧𝑒,𝑖 𝐸(𝑥,𝑦)+𝑧𝑒,𝑖 𝑁(𝑡,𝑥,𝑦) Fundamentos teóricos 20 El problema estacionario puede descomponerse a su vez en problema de espesor y curvatura: 𝑧𝑒𝐸(𝑥,𝑦)=𝑧𝑚 𝐸(𝑥,𝑦)+𝑧𝑠𝐸(𝑥,𝑦) 𝑧𝑖𝐸(𝑥,𝑦)=𝑧𝑚 𝐸(𝑥,𝑦)−𝑧𝑠𝐸(𝑥,𝑦) donde el problema de espesor está denotado por 𝑧𝑠𝐸 y no contribuye a la sustentación del perfil, la cual es aportada en su totalidad por el problema de curvatura, que está denotado por 𝑧𝑚 𝐸. Una descomposición similar puede hacerse para el problema no estacionario. Sin embargo, se va a suponer que para nuestro problema, el espesor de los perfiles del ala no varía con el tiempo, es decir, se mantienen constantes; esto implica que el problema de espesor no incluye términos no estacionarios. Por tanto, para nuestro problema no estacionario sólo queda el problema de curvatura que sí puede presentar variaciones en el tiempo: 𝑧𝑒𝑁(𝑡,𝑥,𝑦)=𝑧𝑝𝑁(𝑡,𝑥,𝑦)=𝑧𝑖𝑁(𝑡,𝑥,𝑦) donde 𝑧𝑝 denota el problema de curvatura para el caso no estacionario A partir de las expresiones para describir la geometría del ala, podemos formular las siguientes condiciones de impenetrabilidad según estemos en el problema estacionario o en el no estacionario: Estacionario: 𝑣𝑧′𝐸(𝑥,𝑦,𝑧=0+)=𝑈∞𝜕𝑧𝑒𝐸(𝑥,𝑦) 𝜕𝑥 𝑣𝑧′𝐸(𝑥,𝑦,𝑧=0−)=𝑈∞𝜕𝑧𝑖𝐸(𝑥,𝑦) 𝜕𝑥 No estacionario: 𝑣𝑧′𝑁(𝑡,𝑥,𝑦,𝑧=0+)=𝑣𝑧′𝑁(𝑡,𝑥,𝑦,𝑧=0−)=𝜕𝑧𝑝𝑁(𝑡,𝑥,𝑦) 𝜕𝑡 +𝑈∞𝜕𝑧𝑝𝑁(𝑡,𝑥,𝑦) 𝜕𝑥 El problema estacionario es ampliamente abordado en las clases de Aerodinámica II de la Escuela Superior de Ingeniería de Sevilla, por lo que en esta memoria no se tratará. El objetivo de este trabajo es centrarse en la resolución del problema no estacionario, para lo cual se extrapolará el método Vortex-Lattice estudiado en el en clase para el problema estacionario al problema no estacionario. 3 MÉTODO VORTEX-LATTICE na vez que hemos enunciado el problema objeto de estudio de esta memoria, se va a proceder a explicar el método numérico que se va a implementar para la resolución del mismo. El objetivo de este capítulo es precisamente desarrollar y justificar el método numérico Vortex-Lattice. Este método es, como ya se ha comentado, una adaptación del método utilizado en las clases de Aerodinámica II (problema 3D estacionario), tomando de las referencias [3] y [9] los elementos necesarios para su desarrollo. Para llevar a cabo la explicación del método, el capítulo se estructura en:  Un primer apartado donde se enuncia el problema no estacionario planteando todas las ecuaciones y condiciones de contorno a las que se llegó en el capítulo 2.  Un segundo apartado en el que se explica la analogía de nuestro problema con el electromagnetismo, de donde se tomarán resultados que serán de vital importancia para la implementación del método Vortex-Lattice.  En el tercer apartado, se desarrollará en vista de los resultados alcanzados en el apartado 3.2 que la solución pasa por llevar a cabo una combinación de hilos de torbellinos. En este apartado se detallará también qué combinación es la más adecuada para esos hilos de torbellinos presentando un mallado para el ala y la estela.  Posteriormente, en el apartado 3.4 se ilustrará el procedimiento de resolución a seguir a partir de los resultados alcanzados en el apartado 3.3.  Finalmente, se incluirá el código implementado en MatLab para la ejecución del método VortexLattice. U “No basta tener buen ingenio, lo principal es aplicarlo bien”. - René Descartes - Método Vortex-Lattice 22 3.1 Ecuaciones del problema incomprensible no estacionario Tal y como se ha visto a lo largo de todo el capítulo 2, el problema no estacionario para el caso incompresible está definido por las siguientes ecuaciones y condiciones de contorno:  Las ecuaciones del problema en términos del campo de velocidades de perturbación 4 se expresan: ∇⋅𝒗 󰇍 󰇍 =0 ∇×𝒗 󰇍 󰇍 =0 (3–1) También puede expresarse en términos del potencial de perturbación, donde es la ecuación de Laplace la que rige nuestro problema: ∇2𝜙=0 (3–2) siendo 𝒗 󰇍 󰇍 =∇𝜙 puesto que el campo de velocidades es conservativo.  Condición de infinito no perturbado: 𝒙→∞, 𝒗→𝟎 (3–3)  Condición para velocidades verticales: 𝑤(𝑡,𝑥,𝑦,𝑧=0±)=𝜕𝑧𝑝(𝑡,𝑥,𝑦) 𝜕𝑡 +𝑈∞𝜕𝑧𝑝(𝑡,𝑥,𝑦) 𝜕𝑥 ,∀(𝑥,𝑦)∈Σ𝑎𝑙𝑎 (3–4)  Continuidad en la estela: 𝑤(𝑡,𝑥,𝑦,𝑧=0+)=𝑤(𝑡,𝑥,𝑦,𝑧=0−); ∀(𝑥,𝑦)∈Σ𝑒𝑠𝑡𝑒𝑙𝑎 (3–5)  Condición de Kutta generalizada: (𝜕 𝜕𝑡+𝑈∞𝜕 𝜕𝑥)Γ(𝑡,𝑥,𝑦)=0 (3–6)  Teorema Bjerkness-Kelvin: La circulación alrededor de cualquier línea fluida cerrada se mantiene constante en el curso del movimiento. 4 Con el objetivo de facilitar la notación, se denota ahora la velocidad de perturbación como 𝒗 en lugar de como 𝒗′, así mismo, las componentes (𝑣𝑥′,𝑣𝑦′,𝑣𝑧′) pasan a denotarse por (𝑢,𝑣,𝑤). 23 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 3.2 Analogía con Electromagnetismo Para hallar la solución a las ecuaciones 3-1, podemos acudir a soluciones ya halladas en problemas gobernados por la misma ecuación diferencial, como por ejemplo el campo magnético generado por un cable por el que circula una corriente de intensidad 𝐼. A partir de las ecuaciones de Maxwell, en concreto, de la ley de Gauss para el campo magnético, sabemos que: ∇⋅𝑩 󰇍 󰇍 =0, siendo 𝑩 el campo magnético. Por otro lado, la ley de Ampère establece, tal y como se enuncia en [10], que la circulación del campo magnético para una curva alrededor de dicho cable cumple: ∮𝑩 󰇍 󰇍 ⋅𝑑𝒍 =𝜇0𝐼 (3–7) Figura 3-1. Ley de Ampère I. Imagen obtenida de [10] Por otro lado, para una curva que no encierre al cable por el que circula la intensidad I, se cumple: ∮𝑩 󰇍 󰇍 ⋅𝑑𝒍 =0 (3–8) Método Vortex-Lattice 24 Figura 3-2. Ley de Ampère II. Imagen obtenida de [10] A partir de este resultado, se puede deducir que el campo magnético cumple ∇×𝑩 󰇍 󰇍 =0 para todo punto exterior al cable por el que circula la corriente de intensidad I: Partimos de: ∮𝑩 󰇍 󰇍 ⋅𝑑𝒍 =0 Por el Teorema de Stokes, se cumple que para la superficie encerrada por la curva cerrada anterior: ∬∇×𝑩 󰇍 󰇍 =0 Si extendemos a un diferencial de superficie esta integral, llegamos a: ∇×𝑩 󰇍 󰇍 =0 (3–9) Por tanto, obtenemos que el campo magnético exterior al cable cumple las siguientes condiciones, que pueden ser relacionadas con el campo de velocidades de perturbación de nuestro problema: ∇⋅𝑩 󰇍 󰇍 =0 → ∇⋅𝒗 󰇍 󰇍 =0 ∇×𝑩 󰇍 󰇍 =0 → ∇×𝒗 󰇍 󰇍 =0 (3–10) Por esta razón, se puede definir análogamente un hilo de torbellinos con intensidad 𝛤 que sustituya al hilo de corriente con intensidad I, de forma que puede establecerse también la siguiente analogía: ∮𝑩 󰇍 󰇍 ⋅𝑑𝒍 =𝜇0𝐼 → Γ=∮𝒗 󰇍 󰇍 ⋅𝑑𝒍 (3–11) A partir de la ley de Biot y Savart, conocemos el campo magnético generado por un hilo de corriente eléctrica y utilizando 3-12, podemos establecer una analogía que nos permita calcular el campo de velocidades que generaría un hilo de torbellinos: 25 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible  Ley de Biot y Savart obtenida de [10]: 𝑩 󰇍 󰇍 (𝒓 󰇍 )=∫𝜇0𝐼 4𝜋 𝐶𝑑𝒍𝟎 󰇍 󰇍 󰇍 ×𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |3,𝑐𝑢𝑚𝑝𝑙𝑒 {∇⋅𝑩 󰇍 󰇍 =0 ∇×𝑩 󰇍 󰇍 =0 (3–12)  Analogía para campo de velocidades de la ley de Biot y Savart: 𝒗 󰇍 󰇍 (𝒓 󰇍 )=∫ Γ 4𝜋 𝐶𝑑𝒍𝟎 󰇍 󰇍 󰇍 ×𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |3,𝑐𝑢𝑚𝑝𝑙𝑒 {∇⋅𝒗 󰇍 󰇍 =0 ∇×𝒗 󰇍 󰇍 =0 (3–13) 3.3 Solución por combinación de hilos de torbellinos De este modo, hemos llegado a una expresión que cumple 3-1, solo nos queda cumplir con las condiciones de contorno, lo cual podrá llevarse a cabo mediante una combinación de diferentes hilos de torbellinos repartidos por la superficie del ala y de la estela, por tanto, el campo de velocidades solución de nuestro problema cumple: 𝒗 󰇍 󰇍 (𝒓 󰇍 )=∑𝒗 󰇍 󰇍 𝑖(𝒓 󰇍 ) 𝑖=∑∮Γi 4𝜋𝑑𝒍𝟎 󰇍 󰇍 󰇍 ×𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |3 𝑖 (3–14) Nótese que la integral es sobre una línea cerrada, esto se debe a que para poder cumplir con la condición de Bjerkness-Kelvin, los hilos de torbellinos no deben comenzar ni terminar en el dominio fluido. Tomando hilos de torbellinos cerrados y de intensidad Γi constante, la condición de Bjerkness-Kelvin queda automáticamente satisfecha según [9]. La solución 3-15 satisface además la condición de contorno establecida en 3-3, tal y como se puede comprobar: 𝒗 󰇍 󰇍 →0,|𝒓 󰇍 |→∞. Si los hilos de torbellinos se colocan en el plano 𝑧=0, la condición 3-5 queda satisfecha también. Por otro lado, si definimos 𝑤𝑝 como: 𝑤𝑝(𝑡,𝑥,𝑦)=(𝜕 𝜕𝑡+𝑈∞𝜕 𝜕𝑥)𝑧𝑝(𝑡,𝑥,𝑦) (3–15) entonces la condición para las velocidades verticales queda: ∑∮Γi 4𝜋𝑑𝒍𝟎 󰇍 󰇍 󰇍 ×𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |3 𝑖⌋𝑧=𝑤𝑝(𝑡,𝑥,𝑦); ∀(𝑥,𝑦)∈Σ𝑎𝑙𝑎 (3–16) donde ∑∮Γi 4𝜋𝑑𝒍𝟎 󰇍 󰇍 󰇍 ×𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |𝒓 󰇍 −𝒓𝟎 󰇍 󰇍 󰇍 󰇍 |3 𝑖⌋𝑧denota la componente ‘z’ del campo de velocidades generado por la distribución de hilos de torbellinos. Por tanto, nuestro objetivo ahora es hallar una distribución de hilos de torbellinos tal que verifique la ecuación 3-17 y la condición de Kutta generalizada 3-6. 3.3.1 Mallado del ala e introducción de hilos de torbellinos El primer paso a la hora de obtener la distribución de hilos de torbellinos que cumpla 3-17 y 3-6, es dividir en paneles el ala y la estela de forma adecuada para posteriormente introducir dichos hilos. Así mismo, el tiempo de simulación será dividido también en 𝑇 instantes de tiempo separados 𝛥t, de forma que los instantes de tiempo 𝑘=1,…,𝑇 suceden cuando 𝑡𝑘 = (𝑘−1)𝛥𝑡. En un instante 𝑘, hay 𝑀×𝑁 paneles en el ala y 𝑇×𝑁 paneles en la estela. Como veremos más adelante, los M Método Vortex-Lattice 32 𝑃=𝑀×𝑁 Π=𝑃+1 𝜏=(𝑀+𝑇)×𝑁 (3–28) El sistema queda introduciendo 3-29: [𝑤11 ⋯ 𝑤𝑃1 ⋮ ⋱ ⋮ 𝑤1𝑃 ⋯ 𝑤𝑃𝑃][Γ1 ⋮ Γ𝑃]=[𝑤𝑝1 ⋮ 𝑤𝑝𝑃]−[𝑤Π1 ⋯ 𝑤𝜏1 ⋮ ⋱ ⋮ 𝑤Π𝑃 ⋯ 𝑤𝜏𝑃][ΓΠ ⋮ Γ𝜏] (3–29) Definiendo las siguientes matrices: 𝑨=[𝑤11 ⋯ 𝑤𝑃1 ⋮ ⋱ ⋮ 𝑤1𝑃 ⋯ 𝑤𝑃𝑃] 𝒙=[Γ1 ⋮ Γ𝑃] 𝑾𝒑=[𝑤𝑝1 ⋮ 𝑤𝑝𝑃] 𝑩=[𝑤Π1 ⋯ 𝑤𝜏1 ⋮ ⋱ ⋮ 𝑤Π𝑃 ⋯ 𝑤𝜏𝑃] 𝒛=[ΓΠ ⋮ Γ𝜏] (3–30) Finalmente, para resolver el problema no estacionario para el caso incompresible no tenemos más que despejar la matriz x para los diferentes instantes de simulación que tomemos: 𝑨𝒙=𝑾𝒑−𝑩𝒛 (3–31) 𝒙=𝑨−𝟏(𝑾𝒑−𝑩𝒛) (3–32) 3.4.2 Cálculo de 𝑪𝑳 Una vez que se ha obtenido la solución al sistema 3-32, se tiene el valor de la circulación en todo el ala, para el cálculo de 𝐶𝐿 utilizando 2-15 necesitamos el valor de la presión en todo el ala. 𝐶𝐿≅1 1 2𝜌∞𝑈∞ 2𝐴𝒛⋅∫∫ (𝑝−𝑝∞)(−𝒏𝑠)𝑑𝜎 Σ𝑎𝑙𝑎 (2-15) La expresión 2-15 puede simplificarse algo más teniendo en cuenta las simplificaciones que si hicieron en los apartados 2.4, 2.5 y 2.6, puesto que el ala ha dejado de tener grosor es una superficie que se extiende en el plano 𝑧=0, de esta forma, podemos suponer que 𝒛≅𝒏𝒔, por lo que: 33 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 𝐶𝐿≅1 1 2𝜌∞𝑈∞ 2𝐴∫∫ (𝑝−𝑝∞)𝑑𝜎 Σ𝑎𝑙𝑎 (3–33) Si se integra por la superficie de extradós e intradós y además se incluye la discretización hecha en el apartado 3.3.1, el cálculo del coeficiente de sustentación se reduce a 3-35: 𝐶𝐿≅1 1 2𝜌∞𝑈∞ 2𝐴∑∑Δ𝑝𝑖𝑗Δ𝑆𝑖𝑗 𝑗𝑖 (3–34) donde Δ𝑝𝑖𝑗 es el salto de presiones entre extradós e intradós que hay en cada panel y Δ𝑆𝑖𝑗 es la superficie de cada panel. Por tanto, se requiere el cálculo del salto de presiones entre extradós e intradós a partir del valor obtenido de la circulación para cada panel del ala. Para ello, utilizamos 2-27 y 2-42 para llegar a la siguiente expresión: Δ𝑝𝑖𝑗=𝜌∞(𝜕 𝜕𝑡+𝑈∞𝜕 𝜕𝑥)Γ𝑖𝑗 (3–35) Para el cálculo numérico de las derivadas, se utilizan las siguientes aproximaciones: 𝜕Γ𝑖𝑗 𝜕𝑡 =Γ𝑖𝑗 𝑘−Γ𝑖𝑗 𝑘−1 Δ𝑡 𝜕Γ𝑖𝑗 𝜕𝑥 =Γij−Γi−1,j ℎ𝑖𝑗 (3–36) donde:  Γ𝑖𝑗 𝑘 es la circulación en el panel (𝑖,𝑗) para el instante de tiempo 𝑘.  Γ𝑖𝑗 𝑘−1 es la circulación en el panel (𝑖,𝑗) para el instante de tiempo 𝑘−1.  Γ𝑖−1,𝑗 𝑘 es la circulación en el panel (𝑖−1,𝑗) para el instante de tiempo 𝑘.  ℎ𝑖𝑗 es la cuerda media del panel (𝑖,𝑗). De esta forma se implementa el cálculo del coeficiente de sustentación a partir de la circulación obtenida para cada panel del ala. 3.5 Implementación en MatLab Para la implementación en MatLab creamos dos módulos principales que se ejecutarán dentro del programa ‘VLattice.m’. El primero de los módulos se encargará de generar las matrices del sistema que se mantienen constantes a lo largo de todo el tiempo de simulación, este módulo denominado ‘gen_mat.m’ se ejecutará una vez al principio de ‘VLattice.m’ para obtener las matrices que se utilizarán durante toda la simulación. El segundo módulo calculará para cada instante de tiempo 𝑘 los valores de la circulación, presión, sustentación y coeficiente de sustentación. Se denominará ‘circ.m’ y será ejecutado en cada instante de simulación. El código de cada módulo así como el del programa ‘VLattice’ se encuentra desarrollado y explicado en los siguientes apartados. Método Vortex-Lattice 34 3.5.1 ‘gen_mat.m’ function [A,B,S,b,xg,yg,xcmat,hmat,bmat,Sp,Mn,Tau,inc_t]=gen_mat(M,N,T,t_final,U_inf,A R,cr,E,psi) %========================================================================== %Función que genera parámetros invariables en el tiempo requeridos por el %método Vortex-Lattice. %========================================================================== %Variables de entrada: %M: Número de paneles en x %N: Número de paneles en y %T: Número de paneles en estela = Número de instantes de tiempo a simular %t_final: Instante de tiempo final de la simulación %U_inf: Velocidad de la corriente incidente %AR: Alargamiento del ala %cr: Cuerda en la raíz del ala %E: Estrechamiento del ala %psi: Ángulo de flecha %========================================================================== %Variables de salida: %A: Matriz A del sistema de ecuaciones 3-32 %B: Matriz B del sistema de ecuaciones 3-32 %S: Superficie del ala %b: Envergadura del ala %xg: Matriz con coordenada x punto centrado de cada panel %yg: Matriz con coordenada y punto centrado de cada panel %hmat: Matriz con cuerda media de cada panel %bmat: Matriz con envergadura de cada panel %Sp: Matriz con superficie de cada panel %P: Parámetro P definido en 3-29 %Tau: Parámetro Tau definido en 3-29 %inc_t: Incremento de tiempo en cada instante de simulación %========================================================================== %INICIO: %Geometría del ala b = AR*cr; S = b*cr; ct = E*cr; %Mallado en t: inc_t = t_final/(T-1); %Mallado en y y(1:(N+1)) = -b/2+((1:(N+1))-1)*b/N; %Mallado en x xa(1:(N+1)) = tan(psi)*abs(y); xs(1:(N+1)) = cr+((b/2*tan(psi)+ct-cr)/(b/2))*abs(y); %Mallado en x e y: xv = zeros(M+1,N+1); yv = zeros(M+1,N+1); for j=1:(N+1) xv(1:(M+1),j) = xa(j)+((1:(M+1))-1)*(xs(j)-xa(j))/M; xv((M+2):(M+T+1),j) = xs(j)+(1:T)*U_inf*inc_t; yv(1:(M+T+1),j) = y(j); end 35 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible %Vértices 1, 2, 3, 4: x1 = zeros(1,(M+T)*N); y1 = zeros(1,(M+T)*N); x2 = zeros(1,(M+T)*N); y2 = zeros(1,(M+T)*N); x3 = zeros(1,(M+T)*N); y3 = zeros(1,(M+T)*N); x4 = zeros(1,(M+T)*N); y4 = zeros(1,(M+T)*N); for i = 1:M+T for j = 1:N I = (i-1)*N+j; x1(I) = xv(i,j); y1(I)=yv(i,j); x2(I) = xv(i,j+1); y2(I)=yv(i,j+1); x3(I) = xv(i+1,j+1); y3(I)=yv(i+1,j+1); x4(I) = xv(i+1,j); y4(I)=yv(i+1,j); end end %Vértices A, B, C, D: xA(1:(M+T)*N) = x1+(x4-x1)/4; yA(1:(M+T)*N)=y1; xB(1:(M+T)*N) = x2+(x3-x2)/4; yB(1:(M+T)*N)=y2; xC(1:(M+T-1)*N) = xB((N+1):(M+T)*N); yC(1:(M+T-1)*N)=yB((N+1):(M+T)*N); xD(1:(M+T-1)*N) = xA((N+1):(M+T)*N); yD(1:(M+T-1)*N)=yA((N+1):(M+T)*N); %Introducimos puntos C y D de la última hilera: xC(((M+T-1)*N+1):(M+T)*N) = xB(((M+T-1)*N+1):(M+T)*N)+U_inf*inc_t; yC(((M+T1)*N+1):(M+T)*N) = yB(((M+T-1)*N+1):(M+T)*N); xD(((M+T-1)*N+1):(M+T)*N) = xA(((M+T-1)*N+1):(M+T)*N)+U_inf*inc_t; yD(((M+T1)*N+1):(M+T)*N) = yA(((M+T-1)*N+1):(M+T)*N); %Punto de colocación: %Cuerda media de cada panel: h(1:(M+T)*N)=(x4+x3-x2-x1)/2; %Coordenadas del punto de colocación: xc(1:(M+T)*N)=(x1+x2)/2+3/4*h; yc(1:(M+T)*N)=(y2+y1)/2; %Introducimos parámetros para notación en matrices: P=M*N; Pi=P+1; Tau=(M+T)*N; %MATRIZ A %Calculamos factores de influencia para matriz A wIC=zeros(P,P); for c=1:P %Recorremos en I el wIk xr=xc(c); yr=yc(c); [wAB_I(1:P)]=BiotSavart 5 (xr,yr,xA(1:P),yA(1:P),xB(1:P),yB(1:P)); %Recorremos en k el wIk [wBC_I(1:P)]=BiotSavart(xr,yr,xB(1:P),yB(1:P),xC(1:P),yC(1:P)); [wCD_I(1:P)]=BiotSavart(xr,yr,xC(1:P),yC(1:P),xD(1:P),yD(1:P)); [wDA_I(1:P)]=BiotSavart(xr,yr,xD(1:P),yD(1:P),xA(1:P),yA(1:P)); wIC(c,1:P)=wAB_I(1:P)+wBC_I(1:P)+wCD_I(1:P)+wDA_I(1:P); end A=wIC; %Matriz B %Calculamos factores de influencia para matriz A wIC=zeros(P,Tau-P); for c=1:P %Recorremos en I el wIk 5 La función ‘BiotSavart.m’ se encuentra explicada en Anexo A. Método Vortex-Lattice 36 xr=xc(c); yr=yc(c); [wAB_I(1:TauP)]=BiotSavart(xr,yr,xA(Pi:Tau),yA(Pi:Tau),xB(Pi:Tau),yB(Pi:Tau)); %Recorremos en k el wIk [wBC_I(1:TauP)]=BiotSavart(xr,yr,xB(Pi:Tau),yB(Pi:Tau),xC(Pi:Tau),yC(Pi:Tau)); [wCD_I(1:TauP)]=BiotSavart(xr,yr,xC(Pi:Tau),yC(Pi:Tau),xD(Pi:Tau),yD(Pi:Tau)); [wDA_I(1:TauP)]=BiotSavart(xr,yr,xD(Pi:Tau),yD(Pi:Tau),xA(Pi:Tau),yA(Pi:Tau)); wIC(c,1:Tau-P)=wAB_I(1:Tau-P)+wBC_I(1:Tau-P)+wCD_I(1:Tau-P)+wDA_I(1:TauP); end B=wIC; %Pasamos de I a (i,j) xg=zeros(M,N); yg=zeros(M,N); Sp=zeros(M,N); hmat=zeros(M,N); bmat=zeros(M,N); for i=1:M for j=1:N I=(i-1)*N+j; %Implementamos coordenadas de los puntos de centrado: xg(i,j)=(x1(I)+x2(I))/2+h(I)/4; yg(i,j)=(y1(I)+y2(I))/2; %Calculamos áreas de los paneles: Sp(i,j)=((x4(I)-x1(I))+(x3(I)-x2(I))).*(y2(I)-y1(I))/2; %Matriz de cuerda media de cada panel: hmat(i,j)=h(I); %Matriz de envergadura de cada panel: bmat(i,j)=(y2(I)+y3(I)-y4(I)-y1(I))/2; end end 37 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 3.5.2 ‘circ.m’ function [x,Gamma,p,Lnum,cLnum]=circ(Gamma_0,A,Wp,B,z,M,N,rho,inc_t,U_inf,hmat,Sp) %========================================================================== %Función que calcula presiones y fuerzas aerodinámicas sobre el perfil a %partir de parámetros específicos de cada instante de tiempo. %========================================================================== %Variables de entrada: %Gamma_0: Matriz de circulaciones en ala en instante k-1 %A: Matriz A del sistema de ecuaciones 3-32 %Wp: Matriz Wp del sistema de ecuaciones 3-32 %B: Matriz B del sistema de ecuaciones 3-32 %z: Vector de circulaciones en estela para instante k %M: Paneles en dirección x %N: Paneles en dirección y %rho: Densidad del aire %inc_t: Incremento del tiempo en cada instante de simulación %U_inf: Velocidad de la corriente incidente %hmat: Matriz de cuerda media por panel %Sp: Matriz de superficie por panel %========================================================================== %Variables de salida: %x: Vector de circulaciones por panel %Gamma: Matriz de circulaciones en cada panel del ala %p: Matriz de diferencias de presión entre extradós e intradós por panel %Lnum: Sustentación en instante de tiempo k %cLnum: Coeficiente de sustentación en instante de tiempo k x=A\(Wp-B*z); %Pasamos de I a (i,j) for i=1:M for j=1:N I=(i-1)*N+j; Gamma(i,j)=x(I,1); end end %==================================================================== %CÁLCULO DE DENSIDADES DE CIRCULACIÓN %==================================================================== %Empezamos calculando la densidad en todos los paneles excepto en los del %borde de ataque: for i=2:M for j=1:N gamma(i,j)=(Gamma(i,j)-Gamma(i-1,j)); end end %Para los paneles del borde de ataque: for j=1:N gamma(1,j)=Gamma(1,j); end gamma=gamma./hmat; %==================================================================== %CÁLCULO SALTO PRESIONES %==================================================================== p=rho*(Gamma-Gamma_0)/inc_t+rho*U_inf*gamma; %==================================================================== %CÁLCULO DE SUSTENTACIÓN Y cL Método Vortex-Lattice 38 %==================================================================== Lnum=sum(sum(Sp.*p)); S=sum(sum(Sp)); %Superficie del ala completa cLnum=Lnum/(1/2*rho*U_inf^2*S); 39 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 3.5.3 ‘VLattice.m’ %PARÁMETROS MALLADO M=50; N=50; Uinf_inct_c=1/16; %Parámetros externos %Caracteristicas aire U_inf=1; rho=1; %Parámetros ala alpha=5*pi/180; AR=10; E=1; psi=0; cr=1; %Partición en el tiempo t_final=20*cr/U_inf; inc_t=Uinf_inct_c*cr/U_inf; T=floor(t_final/inc_t+1); t_final=inc_t*(T-1); [A,B,S,b,xg,yg,xcmat,hmat,bmat,Sp,Mn,Tau,inc_t]=gen_mat(M,N,T,t_final,U_inf,A R,cr,E,psi); %Matriz Wp Wp=[-U_inf*alpha*ones(Mn,1);zeros(Mn-Mn,1)]; %Matriz z y Gamma_0 iniciales z=zeros(Tau-Mn,1); Gamma_0=zeros(M,N); for k=1:T t(k)=(k-1)*inc_t; %Cálculo de Sustentación y cL para instante k [x,Gamma,p,Lnum,cLnum]=circ(Gamma_0,A,Wp,Wg,B,z,M,N,rho,inc_t,U_inf,hmat,Sp); L(k)=Lnum; cL(k)=cLnum; %Actualización de variables para instante k+1 Gamma_0=Gamma; z(((Tau-Mn+1-N):(Tau-Mn)))=[]; z=[x((Mn+1-N):Mn,1);z]; end 4 RESULTADOS ras el desarrollo e implementación del método Vortex-Lattice, llegó el momento de comprobar que los resultados que ofrece son correctos y pueden encontrar una interpretación física. Ese es el objetivo del presente capítulo. La estructura del capítulo consta de:  Un primer apartado donde se exponen los parámetros necesarios para definir la geometría del ala.  Un segundo apartado enunciando y resolviendo el problema de Wagner. En este apartado se incluirán además dos subapartados donde se estudiarán los efectos del estrechamiento y del ángulo de flecha en el ala.  Un tercer apartado donde se tratará el problema de Theodorsen para 3D y se compararán las soluciones obtenidas con el caso 2D y la solución analítica. T Resultados 48 Figura 4-7. Coeficiente de sustentación para diferentes estrechamientos En la figura 4-7, podemos observar que existe un valor óptimo del estrechamiento para el cual se da un 𝑐𝐿 máximo cuando se acerca al estacionario. Si buscamos ese valor óptimo, encontramos que se encuentra en torno a 𝐸=0.3. Conviene destacar además, que a la vista de los resultados expuestos en la Figura 4-7, el peor estrechamiento posible es el del ala rectangular con 𝐸=1, siendo su coeficiente de sustentación el más bajo, esto se debe a que debido al intenso rebordeo que existe de la corriente en las puntas del ala, un ala rectangular aprovecha de forma poco eficiente el flujo incidente. Se puede demostrar (referencia [8]) que es el ala elíptica la que mejor aprovecha el flujo, debido a que minimiza el rebordeo existente en las puntas del ala. Cuanto más parecida sea la geometría de la forma en planta del ala a una elipse, mejor aprovechamiento se hará de la corriente incidente. Efectivamente, si comparamos la forma en planta del ala con 𝐸=0.34, podemos observar en la figura 4-8 que son razonablemente parecidas. 49 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Figura 4-8. Comparación geométrica ala elíptica y hexagonal Por tanto, se puede concluir que el efecto del estrechamiento es beneficioso para el aprovechamiento del flujo incidente siempre que se intente asemejar la forma en planta del ala a la forma elíptica. 4.2.2 Efecto del ángulo de flecha Una vez hemos discutido el efecto del estrechamiento, le toca el turno al ángulo de flecha. Para ello, realizaremos nuevamente diversos análisis sobre las formas en planta de la figura 4-9. Los resultados obtenidos se exponen en la figura 4-10 Resultados 50 Figura 4-9. Formas en planta para distintos estrechamientos incluyendo ángulo de flecha Figura 4-10. Coeficiente de sustentación incluyendo ángulo de flecha En la figura 4-10 se puede observar como el coeficiente de sustentación sigue teniendo una tendencia similar a la observada en la figura 4-7 con la diferencia de que el estacionario al que aspiran las curvas es sustancialmente inferior. La causa de este fenómeno guarda una estrecha relación con el hecho de que la corriente incidente normal que ve cada perfil del ala ahora ha cambiado debido a la introducción de un ángulo de flecha. Para analizar el efecto que tiene ese cambio en la corriente incidente normal a los perfiles del ala, vamos a centrarnos en el caso bidimensional tomando valores del alargamiento convenientemente altos. Teniendo esto en cuenta, se realizan nuevas simulaciones para diferentes valores del ángulo de flecha en estas alas de gran alargamiento. Los resultados obtenidos se reflejan en la figura 4-11 y muestran una tendencia decreciente conforme aumentamos el ángulo de flecha. 51 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Figura 4-11. Coeficiente de sustentación para diferentes ángulos de flecha Curiosamente, en la figura 4-11 se observa que la dependencia de 𝐶𝐿 con 𝜓 puede aproximarse como: 𝐶𝐿(𝜓)=𝐶𝐿(𝜓=0)cos(𝜓). Esta aproximación puede justificarse por el hecho de que la corriente incidente normal que ve cada perfil del ala ahora ha cambiado debido a la introducción de un ángulo de flecha. Efectivamente, la componente normal de la velocidad que ve cada perfil está contenida en un plano que forma vertical que forma un ángulo 𝜓 con la dirección del eje ‘x’ tal y como podemos observar en la figura 2b.3: Figura 4-12. Plano en el que está contenida la componente normal de la velocidad que ve cada perfil Por esta razón, el ángulo de ataque con la corriente incidente normal es algo menor que 𝛼, si hallamos el valor del ángulo de ataque que ve, tendremos que para el caso bidimensional, 𝑐𝑙=2𝜋𝛼𝑒𝑓𝑓. 𝑥 𝜓 𝜓 Plano de ecuación 𝑦=−𝑥 𝑡𝑎𝑛𝜓 𝑦 Resultados 52 Para encontrar 𝛼𝑒𝑓𝑓 se puede hacer un análisis geométrico, del que se extrae: sen𝛼𝑒𝑓𝑓=tan𝛼 √1+tan2𝜓+tan2𝛼 De la expresión anterior, puede hacerse la siguiente aproximación: sen𝛼𝑒𝑓𝑓≅𝛼𝑒𝑓𝑓 𝑝𝑎𝑟𝑎 𝛼𝑒𝑓𝑓≪1 tan𝛼≅𝛼 𝑝𝑎𝑟𝑎 𝛼≪1 √1+tan2𝜓+tan2𝛼≅√1+tan2𝜓 𝑝𝑎𝑟𝑎 𝛼≪1 Por tanto, la expresión para 𝛼𝑒𝑓𝑓 nos queda: 𝛼𝑒𝑓𝑓=𝛼 √1+tan2𝜓=𝛼cos𝜓 Es decir, para el 𝑐𝑙 para el caso bidimensional que tiene cada perfil para el ala en flecha es: 𝑐𝑙=2𝜋𝛼𝑒𝑓𝑓=2𝜋cos𝜓 𝛼 Por lo que, el 𝑐𝐿 al que tiende el ala en flecha para alargamientos muy altos es 𝐶𝐿(𝜓)=𝐶𝐿(𝜓=0)cos(𝜓), lo cual efectivamente se comprueba al observar la figura 4-11. De este modo, se concluye que la inclusión de un ángulo de flecha en el ala resulta en un decremento del coeficiente de sustentación al que se puede aspirar. Sin embargo, y para concluir, se sabe que el ángulo de flecha tiene otros efectos beneficiosos como:  Favorecer la reducción de efectos transónicos en el ala, ya que la componente normal de la velocidad es la que más afecta a los perfiles. Los efectos transónicos son muy perjudiciales para la sustentación, puesto que tienden a igualar presiones entre extradós e intradós del ala, lo cual nos reduce la sustentación y además, ocasiona la aparición de momentos de cabeceo. Por tanto, es conveniente evitarlos.  Ofrecer una mayor estabilidad a balanceo. Al incluir dicho ángulo, aparece una velocidad de resbalamiento al aplicar un momento de balanceo que tiende a aumentar la sustentación de un ala y a reducir la de la otra de forma que genera un momento recuperador oponiendose al momento ejercido.  Garantizar una mejor estabilidad a cabeceo, ya que tiende a retrasar el centro aerodinámico con respecto al centro de gravedad, lo cual nos interesa para generar un momento recuperador gracias al peso cuando aumentamos el ángulo de ataque. Estas son las razones por las que se incluye un ángulo de flecha en alas, a pesar de la reducción que puede ocasionar en el coeficiente de sustentación aerodinámica. 4.3 Problema de Theodorsen El problema de Theodorsen para la Aerodinámica no estacionaria, está caracterizado por el movimiento armónico de una placa teniendo dos grados de libertad: ℎ, que es el desplazamiento vertical (positivo hacia abajo) del eje de giro (también denominado eje elástico) y 𝛼, que es el giro alrededor de dicho eje. La posición del eje de giro se denota por 𝑥𝑒. Como el movimiento es armónico, cualquier magnitud genérica 𝜓 puede expresarse mediante un fasor 𝜓 tal que: 𝜓=ℜ( 𝜓𝑒𝑗𝜔𝑡). La sustentación de un perfil bidimensional puede ser determinada gracias a la fórmula de Theodorsen, desarrollada en [6], a través de los fasores ℎ,𝛼 y las variables 𝜔 y 𝑥𝑒. Al igual que en el apartado 4.2, aquí tomaremos distintos valores del alargamiento y comprobaremos que para alargamientos relativamente altos, los resultados de la solución numérica dada por el método Vortex-Lattice se 53 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible aproximan razonablemente a la solución analítica de Theodorsen. Figura 4-13. Problema de Theodorsen La solución analítica puede obtenerse a partir del fasor de 𝑧𝑝: 𝑧𝑝 =−ℎ−𝛼(𝑥−𝑥𝑒) Utilizando 3-16, nos queda: 𝑤𝑝 =−𝑗𝜔𝑧𝑝 +𝑈∞𝜕𝑧𝑝  𝜕𝑥=−𝑗𝜔ℎ−[𝑈∞+𝑗𝜔(𝑥−𝑥𝑒)]𝛼 Introduciendo en el programa, 𝑤𝑝=ℜ(𝑤𝑝 𝑒𝑗𝜔𝑡), obtenemos la solución representada en la figura 4-14. Además, conforme se aumenta el valor del alargamiento, más se aproxima la solución obtenida al caso 2D. Para obtener la solución 2D que aparece, se ha tomado un valor de estrechamiento elevado (Λ=100) y como se puede comprobar prácticamente coincide con la solución analítica predicha por Theodorsen. 𝑦 𝑧𝑝 ℎ 𝑥𝑒 𝛼 Resultados 54 Figura 4-14. Coeficiente de sustentación para problema de Theodorsen con diferentes valores de Λ 5 FLAMEO AEROELÁSTICO al y como se adelantaba en el capítulo 1, el estudio de la Aerodinámica no estacionaria puede ayudar a predecir la aparición de fenómenos perjudiciales para la estructura del ala. El principal efecto perjudicial es la aparición del flameo, que como ya se comentó, supone la aparición de oscilaciones armónicas que no sufren amortiguamiento alguno con el tiempo. Este fenómeno se debe a que el ala comienza a extraer energía del fluido, extrae la misma cantidad que pierde a causa de las fuerzas de amortiguamiento. Por tanto, se va a modelar el ala de forma que se pueda predecir la aparición de este fenómeno, la principal diferencia con los resultados obtenidos hasta ahora es que la posición de los perfiles del ala es desconocida en este caso. En los resultados presentados en el capítulo 4, la función 𝑧=𝑧𝑝(𝑡,𝑥,𝑦) era conocida, sin embargo, para esta nueva aplicación, la función 𝑧=𝑧𝑝(𝑡,𝑥,𝑦) describirá un movimiento oscilatorio causado por el acoplamiento de fuerzas elásticas, de amortiguamiento y aerodinámicas. Para llevar a cabo este análisis, el capítulo se estructurará en tres apartados:  Un primer apartado donde se expondrá un modelo de dos grados de libertad para un ala rígida similar al problema de Theodorsen del apartado 4.3.  Un segundo apartado donde se expondrá la ecuación del movimiento para una placa elástica y la metodología seguida para resolverla.  Y finalmente un tercer apartado donde se discutirán los resultados obtenidos por los dos modelos. 5.1 Placa plana rígida con dos grados de libertad El primer modelo que se va a desarrollar para el flameo estará formado por una placa plana rígida que tiene dos grados de libertad: ℎ y 𝛼. La placa tiene una masa 𝑚 con el centro de gravedad en (𝑥=𝑥𝐺,𝑦=0) y con momento de inercia 𝐼𝐺 respecto del centro de gravedad. La placa estará sujeta con dos muelles de rigideces 𝑘ℎ y 𝑘𝛼 tal y como se muestra en la figura 5-1. T Flameo Aeroelástico 56 Figura 5-1. Placa plana con dos grados de libertad (ℎ,𝛼) Las ecuaciones que rigen el movimiento de la placa se derivan de aplicar las ecuaciones de Lagrange y son: 𝑚ℎ󰇘=−𝑘ℎℎ−𝐿−𝑚(𝑥𝐺−𝑥𝑒)𝛼󰇘 𝐼𝐺𝛼󰇘=𝑘ℎℎ(𝑥𝐺−𝑥𝑒)−𝑘𝛼𝛼+𝑀𝐿 (5–1) donde 𝐿 y 𝑀𝐿 son fuerza de sustentación total sobre la placa y el momento que genera esa fuerza con respecto al centro de gravedad. 5.1.1 Cálculo de fuerza y momento De las expresiones 5-1, se busca despejar ℎ󰇘 y 𝛼󰇘 para lo cual es necesario conocer previamente los valores de la sustentación y el momento para un instante de tiempo dado. Para el cálculo de 𝐿 se utilizará la expresión 5-2 que está basada en 3-34: 𝐿=1 2𝜌∞𝑈∞ 2𝐴𝐶𝐿=∫∫ (𝑝−𝑝∞)𝑑𝑥𝑑𝑦 Σ𝑝𝑙𝑎𝑐𝑎 (5–2) Para calcular 𝑀𝐿, hacemos uso de la definición de momento ejercido por una fuerza y realizamos la integración en la superficie de la placa: 𝑀𝐿=∫∫ |(𝒙 󰇍 󰇍 −𝒙 󰇍 󰇍 𝑮)×(−𝒏 󰇍 󰇍 𝒔)(𝑝−𝑝∞)|𝑑𝜎 Σ𝑝𝑙𝑎𝑐𝑎 (5–3) 𝑘𝛼 𝑘ℎ 57 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible donde 𝒙 󰇍 󰇍 y 𝒙 󰇍 󰇍 𝑮, representan los vectores posición de un punto genérico del ala (la variable de integración) y la del centro de gravedad de la placa respectivamente. Teniendo en cuenta que la distribución de presiones en la placa será simétrica con respecto al eje 𝑥, la expresión 5-3 puede reducirse a: 𝑀𝐿=∫∫ (𝑝−𝑝∞)(𝑥−𝑥𝐺)𝑑𝜎 Σ𝑝𝑙𝑎𝑐𝑎 (5–4) Para el cálculo numérico de estas expresiones, se va a utilizar la matriz de presiones Δ𝑝𝑖𝑗 obtenida de 3-36, quedando finalmente para un instante de tiempo 𝑡𝑘 dado: 𝐿(𝑡𝑘)≅∑∑Δ𝑝𝑖𝑗(𝑡𝑘)Δ𝑆𝑖𝑗 𝑗𝑖 𝑀𝐿(𝑡𝑘)≅∑∑Δ𝑝𝑖𝑗(𝑡𝑘)Δ𝑆𝑖𝑗(𝑥𝑔𝑖𝑗+𝑥𝑐𝑖𝑗 2−𝑥𝐺) 𝑗𝑖 (5–5) donde 𝑥𝑔𝑖𝑗 y 𝑥𝑐𝑖𝑗 representan las posiciones de los puntos de centrado y colocación de cada panel. Para el cálculo de Δ𝑝𝑖𝑗, es necesario imponer 𝑤𝑝, sabiendo que 𝑧𝑝(𝑡,𝑥,𝑦)=ℎ(𝑡)+𝛼(𝑡)(𝑥𝑒−𝑥) 𝑤𝑝(𝑡,𝑥,𝑦)=𝜕𝑧𝑝 𝜕𝑡+𝑈∞𝜕𝑧𝑝 𝜕𝑥=ℎ󰇗(𝑡)+𝛼󰇗(𝑡)(𝑥𝑒−𝑥)−𝑈∞𝛼(𝑡) 5.1.2 Integración en el tiempo A partir de las ecuaciones 5-1, se busca calcular ℎ(𝑡) y 𝛼(𝑡) a partir de las siguientes integrales: 𝛼󰇗(𝑡)=∫𝛼󰇘(𝑡′)𝑑𝑡′ 𝑡 0 𝛼(𝑡)=∫𝛼󰇗(𝑡′)𝑑𝑡′ 𝑡 0 (5–6) ℎ󰇗(𝑡)=∫ℎ󰇘(𝑡′)𝑑𝑡′ 𝑡 0 ℎ(𝑡)=∫ℎ󰇗(𝑡′)𝑑𝑡′ 𝑡 0 (5–7) Para el cálculo numérico de 5-6 y 5-7 utilizando 5-1, se implementarán las siguientes operaciones para el cálculo de ℎ(𝑡𝑘) y 𝛼(𝑡𝑘) en un instante de tiempo 𝑡𝑘 dado y siendo las magnitudes en 𝑡𝑘−1 conocidas: 𝛼󰇘(𝑡𝑘)=1 𝐼𝐺(𝑀𝐿(𝑡𝑘)+𝑘ℎℎ(𝑥𝐺−𝑥𝑒)−𝑘𝛼𝛼(𝑡𝑘−1)) 𝛼󰇗(𝑡𝑘)=𝛼󰇗(𝑡𝑘−1)+𝛼󰇘(𝑡𝑘)Δ𝑡 𝛼(𝑡𝑘)=𝛼(𝑡𝑘−1)+𝛼󰇗(𝑡𝑘−1)Δ𝑡 (5–8) Flameo Aeroelástico 64 Sin embargo, para poder aplicar 5-32, primero hay que calcular la media de presiones 〈Δ𝑝𝑖〉 en la envergadura para el perfil central, y ese es el objetivo del presente subapartado. Para hallar 〈Δ𝑝〉, como se trata de una media, se realizará el siguiente cálculo: 〈Δ𝑝(𝑥,𝑡)〉=1 𝐻∫ 𝑝(𝑥,𝑦,𝑡)𝑑𝑦 𝐻/2 −𝐻/2 (5–33) La expresión 5-33 puede aproximarse numéricamente: 〈Δ𝑝𝑖𝑘〉=〈Δ𝑝(𝑥𝑖,𝑡𝑘)〉=1 𝐻∑𝑝(𝑥𝑖,𝑦𝑗,𝑡𝑘)𝑏𝑖𝑗 𝑁 𝑗=1 =1 𝑁∑𝑝(𝑥𝑖,𝑦𝑗,𝑡𝑘) 𝑁 𝑗=1 =1 𝑁∑Δ𝑝𝑖𝑗 𝑘 𝑁 𝑗=1 (5–34) donde 𝑏𝑖𝑗 es la longitud en dirección de la envergadura de cada panel. Para hallar Δ𝑝𝑖𝑗 𝑘, se utiliza como siempre 3-36, que requiere el cálculo de Γ𝑖𝑗 en cada panel, resultado que obtenemos a partir de imponer 𝑤𝑝. El cálculo de 𝑤𝑝 queda: 𝑤𝑝=𝜕𝑧𝑝∗ 𝜕𝑡∗+𝑈∞𝜕𝑧𝑝∗ 𝜕𝑥∗=∑𝑞𝑖󰇗(𝑡∗)𝜓𝑖(𝑥∗) 𝑛 𝑖=1 +𝑈∞∑𝑞𝑖(𝑡∗)𝜕𝜓𝑖(𝑥∗) 𝜕𝑥∗ 𝑛 𝑖=1 (5–35) Poniendo 5-35 en forma matricial se obtiene finalmente: 𝑾𝒑 𝑘=𝑈∞ [ 𝜕𝜓1(𝑥𝑐1 ∗) 𝜕𝑥∗⋯𝜕𝜓𝑛(𝑥𝑐1 ∗) 𝜕𝑥∗ ⋮ ⋱ ⋮ 𝜕𝜓1(𝑥𝑐𝑀 ∗) 𝜕𝑥∗⋯𝜕𝜓𝑛(𝑥𝑐𝑀 ∗) 𝜕𝑥∗ ] 𝒒𝒌+[𝜓1(𝑥𝑐1 ∗)⋯ 𝜓𝑛(𝑥𝑐1 ∗) ⋮ ⋱ ⋮ 𝜓1(𝑥𝑐𝑀 ∗)⋯ 𝜓𝑛(𝑥𝑐𝑀 ∗)]𝒒󰇘𝒌 (5–36) 5.2.6 Implementación en MatLab Al igual que se hizo para la implementación del método Vortex-Lattice en el capítulo 3, para la implementación de este método se utilizarán diversos módulos para hacer más fácil tanto la programación como el entendimiento del código. Los nuevos módulos que se implementan son:  ‘mat_psi.m’ que calcula las matrices que tiene por elementos los modos de vibración evaluados en los puntos de centrado y colocación del perfil, para poder ser utilizadas en 5-32 y5-36.  ‘wp.m’ que aplica 5-36 a través del vector 𝒙=[𝒒 𝒒󰇗].  ‘media.m’ calcula la media de las presiones en la envergadura a través de 5-34.  ‘pesos.m’ es una función diseñada a partir del apéndice B de [3] donde se expone un método de integración realizando una aproximación a través de un determinado número de pasos de integración. El código final para la resolución de este problema queda: 65 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible %Parámetros externos: Mast=0.74; U=5.09; AR=10; E=1; psi=0; L=1; rho=1; Uinf=1; n=3; %Número de modos de vibración p=2; %Pasos de integración %Inicio: qk=[0.04 0.04 0.04]'; qdot=zeros(n,1); xdot=zeros(2*n,0); x=[qk;qdot]; N=40; M=50; %Paneles en ala T=2000; %Paneles en estela t_final=10; %Cálculo matrices [A,B,S,b,xg,yg,hmat,bmat,Sp,Mn,Tau,inc_t]=gen_mat(M,N,T,t_final,Uinf,AR,L,E,p si); %Inicialización simulación z=zeros(Tau-Mn,1); Gamma_0=zeros(M,N); %Matrices psi [Psic,Psicprima,Psig,beta]=mat_psi(M,n,L); %Vector tiempo t=linspace(0,t_final,T); %========================================================================== %SIMULACIÓN %========================================================================== for k=1:T %Matriz Wp Wp=wp(Uinf,Psicprima,Psic,x,M,N); %Cálculo de presiones [gamsol,Gamma,presion,Lnum,cLnum]=circ(Gamma_0,A,Wp,B,z,M,N,rho,inc_t,Uinf,hm at,bmat,Sp); %Cálculo de presiones medias pmed=media(M,N,presion,bmat,b); %Actualización de parámetros para siguiente iteración Gamma_0=Gamma; z(((Tau-Mn+1-N):(Tau-Mn)))=[]; z=[gamsol((Mn+1-N):Mn,1);z]; %Cálculo de aceleraciones Qaero=Psig*L/M*pmed; q2dot=Mast*Qaero-1/U^2*diag(beta)*qk; xdot=[xdot,[qdot;q2dot]]; Flameo Aeroelástico 66 %Cálculo pesos integración if k<=p pesos_integracion=pesos(t,k); else xdot(:,1)=[]; end %Integración en tiempo x=[qk;qdot]; x=x+xdot*pesos_integracion; qk=x(1:n); qdot=x(n+1:2*n); q(:,k)=qk; end 5.3 Resultados obtenidos Una vez se han desarrollado e implementado los dos modelos ideados para la caracterización del flameo en un ala, llegó el momento de analizar los resultados que se obtienen. Los resultados obtenidos para el primer modelo, muestran que un ala finita tiende a tener su velocidad de flameo más baja conforme aumenta la envergadura, es decir, las alas de menor envergadura son más estables, o dicho de otro modo, un ala finita presenta una mayor estabilidad que un ala infinita. Este resultado es mostrado por [5] y diversos autores también citados en dicha referencia; en ese aspecto, el modelo del ala rígida con dos grados de libertad nos da una idea cualitativa bastante buena del comportamiento aeroelástico de un ala. Además, si se representa la evolución de la velocidad de flameo con respecto al estrechamiento, se observan tendencias muy parecidas a las que aparecen en la bibliografía consultada. Por otro lado, un ala infinitamente rígida como la expuesta en este modelo no se da en la realidad, por lo que cabría esperar que los resultados obtenidos por el segundo modelo sean cuantitativamente mejores. Sin embargo, la inestabilidad que presenta el método de integración utilizado cuando la frecuencia de la respuesta del sistema es muy alta ha hecho imposible obtener soluciones satisfactorias por parte del segundo modelo. Además, no sólo la inestabilidad del método de integración ha sido el problema, la partición que se puede realizar en el tiempo con el método Vortex-Lattice diseñado está limitada por la memoria que puede utilizar MatLab, puesto que la partición en el tiempo está íntimamente relacionada con el mallado de la estela, a mayor partición en el tiempo, más paneles se incluyen en la estela y mayores son las matrices con las que tiene que trabajar el programa. En definitiva, la partición en el tiempo máxima que se ha podido realizar no ha sido capaz de hacer frente a la inestabilidad presente en el método de integración, dicha inestabilidad sólo puede solucionarse utilizando un método más robusto o consiguiendo ser capaces de incrementar aún más la partición en el tiempo mediante alguna modificación en el programa Vortex-Lattice. Todo esto ha provocado que en el segundo modelo sólo se puedan obtener resultados aislados para los parámetros 𝑀∗ y 𝑈∗, a partir de los cuales no se pueden establecer relaciones claras. Lo que sí se puede asegurar es que los modos de vibración del ala en esos casos particulares en los que el método ofrece resultados concuerdan con bastante precisión con las observaciones experimentales. A continuación, se desarrollan dos subapartados donde se exponen los resultados obtenidos. 5.3.1 Resultados de placa plana rígida con dos grados de libertad Tras ejecutar el programa mostrado en 5.1.3, se obtienen las siguientes gráficas: 67 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Figura 5-3. Respuesta para la velocidad de flameo 𝑈𝐹 del ala de alargamiento Λ=100 Figura 5-4. Respuesta para una velocidad superior a la de flameo 𝑈𝐹 del ala de alargamiento Λ=100 Flameo Aeroelástico 68 Figura 5-5. Respuesta para una velocidad inferior a la de flameo 𝑈𝐹 del ala de alargamiento Λ=100 A partir de las figuras 5-3, 5-4 y 5-5 se puede observar como existe una determinada velocidad 𝑈𝐹 para la que el ala comienza a tener oscilaciones armónicas que no se amortiguan con el tiempo (figura 5-3). Por encima de dicha velocidad, las oscilaciones divergen hacia al sistema inestable (figura 5-4), esto se debe a que la placa está extrayendo del flujo incidente más energía de la que pierde por amortiguamiento. Y finalmente, por debajo de la velocidad de flameo, las oscilaciones se amortiguan (figura 5-5) haciendo al sistema estable, la interpretación física de este fenómeno se basa en el hecho de que el perfil no extrae la suficiente energía del fluido como para compensar las pérdidas debidas al amortiguamiento. Por otro lado, si mantenemos la velocidad de 𝑈∞=0.55, que resulta ser la de flameo para Λ=100, y reducimos el alargamiento del ala, obtenemos los resultados mostrados en la figura 5-6. 69 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Figura 5-6. Respuesta para 𝑈∞=0.55 para un ala de alargamiento Λ=10 En la figura 5-6, se puede observar como la respuesta del sistema tiende a ser amortiguada para Λ=10, es decir, se encuentra por debajo de su velocidad de flameo, a pesar de estar utilizando la misma velocidad que resultaba ser de flameo para el ala de alargamiento Λ=100. Este resultado refleja lo que ya se comentaba al principio del apartado, cuanto menor sea la envergadura, más estable será el ala. Esto implica que un ala finita es más estable que un ala infinita tal y como se enuncia en [5]. Se puede representar la evolución de 𝑈𝐹 frente al alargamiento observando así la comentada tendencia. Esta evolución ha sido representada en la figura 5-7. Flameo Aeroelástico 70 Figura 5-7. Evolución de 𝑈𝐹 frente al alargamiento Λ Si se compara la figura 5-7 con la figura 5-8 extraída de [5], se puede comprobar que la tendencia obtenida es muy parecida a la que se observa en dicha referencia. Se puede ver como para alargamientos del orden de 0.1, la velocidad de flameo se incrementa considerablemente; mientras que para alargamientos en torno a 10-100, la variación no es tan significativamente notable, a pesar de que se sigue obteniendo la tendencia de que a menor envergadura, más estabilidad hay. Figura 5-8. Evolución de 𝑈𝐹 frente al alargamiento Λ según [5] 71 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible 5.3.2 Resultados de placa flexible semiempotrada La ejecución del programa implementado en 5.2.6 no ha tenido resultados satisfactorios a causa de la gran inestabilidad numérica que presenta el método explícito de integración por pesos utilizado y desarrollado en [3]. En la referencia [3] ya se comentaba dicho problema, que podía ser subsanado utilizando un método de integración implícito o aumentando la partición en el tiempo. El método Vortex-Lattice desarrollado para el caso tridimensional no permite la implementación del método de integración implícito y la partición en el tiempo se ha escogido lo más fina que permite la potencia de MatLab. A pesar de una utilizar una partición tan fina para el tiempo, los resultados obtenidos son muy inestables para la gran mayoría de combinaciones de los parámetros 𝑀∗ y 𝑈∗. Sin embargo, hay ciertas combinaciones de estos parámetros para los que la frecuencia de la solución no resulta excesivamente alta y el método consigue realizar la integración en el tiempo correctamente. Los resultados que se pueden obtener para los modos de vibración se pueden observar en la figura 5-9 y si se comparan con la figura 5-10 se puede comprobar que guardan mucho parecido con los resultados experimentales. Figura 5-9. Vista de perfil de la vibración de una placa flexible Flameo Aeroelástico 72 Figura 5-10. Vista experimental de perfil de la vibración de una placa flexible Este gran parecido se debe a que los modos de vibración están convenientemente bien elegidos, por lo que para unas determinadas oscilaciones de sus coordenadas generalizadas, los resultados obtenidos guardan una gran similitud con las observaciones experimentales. Sin embargo, el hecho de que estos resultados se hayan encontrado de forma tan puntual no permite establecer analogía alguna con otras referencias como [4] y [5] donde se analiza la dependencia de 𝑈∗ con 𝐻∗ y 𝑀∗. Para asegurar que el problema que aparecía se debía a la inestabilidad del método de integración, se resolvió con el mismo método un problema del que se conociera su solución analítica, como por ejemplo, la vibración libre frente a condiciones iniciales no nulas de un sistema mecánico. Los resultados resultan bastante aproximados cuando la frecuencia de la respuesta no es excesivamente alta como se muestra en la figura 5-11, sin embargo, cuando se impone una frecuencia mayor y no se aumenta lo suficiente la partición en el tiempo, la solución numérica diverge considerablemente de la solución analítica tal y como sucede para el caso del flameo, este fenómeno queda reflejado en la figura 5-12. Se remarca nuevamente que dicha divergencia de la solución numérica se produce para nuestro problema con una partición en el tiempo que no puede aumentarse más debido a la limitación de memoria impuesta por MatLab a la hora de realizar el mallado de la estela. 73 Implementación del método Vortex-Lattice para el cálculo de la Aerodinámica no Estacionaria de alas en régimen incompresible Figura 5-11. Resultados del método de integración para frecuencias bajas Figura 5-12. Resultados del método de integración para frecuencias altas Referencias 80 REFERENCIAS [1] S.S. Arora: Study of Vibration Characteristics of Cantilever Beams of Different Materials. Thapar University, Patiala, 2012. [2] A. Barrero Ripoll, M. Pérez-Saborid: Fundamentos y Aplicaciones de la Mecánica de Fluidos, McGrawHill.. [3] M. Colera Rico: Proyecto Fin de Carrera, Cálculo numérico del flujo potencial linealizado no estacionario sobre perfiles en los regímenes subsónico y supersónico , Escuela Superior de Ingeniería, Sevilla. 2015. [4] C. Eloy, R. Lagrange, C. Souillez, L. Schouvellier: Aeroelastic instability of cantilevered flexible plates in uniform flow. Journal of Fluid Mechanichs, Vol 611, pags 97.106, 2006. [5] C. Eloy, C. Souilliez, L. Schouveiler: Flutter of a rectangular plate. Volume 23, Issue 6, Pages 904-919, Journal of Fluids and Structures.2007. [6] Y.C. Fung: An Introduction to the Theory of Aeroelasticity. Dover, 2008. [7] I. E. Garrick: Propulsion of a Flapping and Oscillating Airfoil. NACA Tech. Rep. 567, 1936. [8] J. M. Gordillo Arias de Saavedra, G. Riboux Acher: Introducción a la Aerodinámica Potencial, 1ª ed, Paraninfo, 2012. [9] J. Katz, A. Plotkin: Low Speed Aerodynamics. Cambridge University Press. 2001. [10] Young, Hugh D. y Roger A. Freedman: Física Universitaria, con física moderna, volumen 2. Pearson Educación.