Diseño de perfiles aerodinámicos mediante metodología inversa

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado de Ingeniería Aeroespacial Diseño de perfiles aerodinámicos mediante metodología inversa Autor: Aitor Robles Corpa 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, 2015 Trabajo Fin de Grado Grado de Ingeniería Aeroespacial Diseño de perfiles aerodinámicos mediante metodología inversa Autor: Aitor Robles Corpa Tutor: Miguel Pérez-Saborid Sánchez-Pastor Profesor Titular Dep. Ingeniería Aeroespacial y Mecánica de Fluidos Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2015 Trabajo Fin de Grado: Diseño de perfiles aerodinámicos mediante metodología inversa Autor: Aitor Robles Corpa Tutor: Miguel Pérez-Saborid Sánchez-Pastor El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: A mi familia, por todo su apoyo, comprensión y cariño durante todos estos años. A mis amig@s, por esos largos días de biblioteca resolviendo problemas juntos. A Pilar y Saúl, por enseñarme que las matemáticas son algo más que números y operaciones. A mi tutor Miguel Pérez-Saborid, por su orientación, apoyo y dedicación. A todos, muchas gracias Resumen El objetivo de este proyecto es el desarrollo de un programa mediante Matlab que sea capaz de calcular las características geométricas de un perfil aerodinámico a partir de una distribución de presiones Cp especificada. Con el fin de conseguir la mejor precisión posible se utilizan tres metodologías diferentes: la primera basada en el método de paneles de Vortex-Lattice, la segunda con series de Fourier y la tercera adapta la fórmula de inversión de Söhngen. Los códigos han sido desarrollados considerando flujo potencial, es decir, no tienen en cuenta ni efectos viscosos ni compresibilidad del fluido. La comparación de los resultados se realiza viendo si los perfiles calculados reproducen las secciones buscadas. Para ello, se toman una serie de perfiles en distintas condiciones, y se introducen las presiones del flujo potencial en los programas. Las distribuciones de presiones son obtenidas bien mediante el método de Green implementado por Gordillo[11], o usando el software libre XFLR5. III XÍndice de Figuras 3.3. Cambio de variable x→θ38 3.4. Comparación entre mallas de integración 41 3.5. Integración de la pendiente dze dx 41 3.6. Esquema de programación 43 3.7. Comparación de distribuciones de presiones con flujo potencial a α=0° 44 3.8. Triángulos semejantes 45 3.9. Comparación de distribuciones de presiones tras la aproximación lineal en el borde de salida 46 4.1. Perfil lenticular 47 4.2. Distribución de presiones asociada a perfiles parabólicosCp(x,t)48 4.3. Perfiles parabólicos con t=0.149 4.4. Perfiles parabólicos con t=0.249 4.5. Perfil elíptico 50 4.6. Distribución de presiones asociada a perfiles elípticos Cp(x,t)51 4.7. Perfiles elípticos con t=0.251 4.8. Perfiles elípticos con t=0.652 4.9. NACA 0012 , α=0° 54 4.10. NACA 1412 , α=0° 54 4.11. NACA 2412 , α=0° 54 4.12. NACA 0012 , α=2° 55 4.13. NACA 1412 , α=2° 55 4.14. NACA 2412 , α=2° 55 4.15. NACA 0024 , α=0° 56 4.16. NACA 1424 , α=0° 56 4.17. NACA 4420 , α=0° 56 4.18. NACA 0024 , α=2° 57 4.19. NACA 1424 , α=2° 57 4.20. NACA 4420 , α=2° 57 4.21. Estudio del mallado para NACA 0012 , α=2° 58 4.22. Estudio de los parámetros de entrada para NACA 1412 , α=0° 60 4.23. Estudio del ángulo de ataque para NACA 0012 61 4.24. Estudio de la curvatura y el espesor con perfiles NACA 4 dígitos, α=0° 63 4.25. Comparación flujo real y potencial: NACA 4412 , α=0° 65 4.26. S. Fourier: Clark Y , α=0° 66 4.27. S. Fourier: Eppler 211 , α=0° 66 4.28. S. Fourier: Göttingen 394 , α=0° 66 4.29. Comparación de métodos de diseño inverso , Boeing 727, α=0° 67 A.1. Pantalla inicial del XFLR5 74 A.2. Pantalla inicial de diseño por método directo 74 A.3. Diseño de perfiles: NACA 1424 y NACA 0012 con flaps y slaps 75 A.4. Pantalla inicial del análisis de perfiles 76 A.5. Configuración del tipo de análisis 76 A.6. Polar view: características globales 77 A.7. Distribucion de presiones Cppara NACA 1424 a α=3° 78 A.8. Presentación de los datos exportados 78 B.1. Superficies para la integral de Green [11] 79 B.2. Discretización de la geometría en Npaneles. Definición de punto medio y vectores tangente y normal de un panel k[11] 80 Índice de Tablas 1.1. Serie NACA: ventajas, inconvenientes y aplicaciones [1] 6 4.1. Tiempos de ejecución al variar el mallado 59 XI Introducción Si por algo se distingue el ingeniero aeronáutico del resto de ingenieros es por sus conocimientos sobre Aerodinámica. Es la clave que permite explicar y comprender los diversos fenómenos que provoca una corriente fluida sobre un objeto, y cómo aprovecharlos en nuestro beneficio. Existen multitud de superficies que aprovechan estas fuerzas aerodinámicas para distintas aplicaciones como: álabes de compresor, alabes de turbinas de gas y vapor, hélices, rotores, alas, estabilizadores, alerones, spoilers, etc. No obstante, todas ellas parten de la misma base, los perfiles aerodinámicos. La motivación principal de este proyecto consiste en dar un uso práctico a los conocimientos y herramientas que se han adquirido durante el Grado de Ingeniería Aeroespacial en materia de perfiles. En distintas asignaturas se define la geometría y las propiedades de las secciones, se estudia el comportamiento del flujo a su alrededor, se explican los efectos de la capa límite, y muchas otras cosas; pero en ningún momento se dedica tiempo para hablar sobre cómo se diseñan. Por esta razón surge mi interés, y el origen de este proyecto. Por otro lado, el motivo de encaminar el trabajo hacia un programa de cálculo es demostrar que hoy en día, bajo la guía de un profesor, cualquier alumno puede crear sus propias rutinas para resolver problemas reales de ingeniería. Tradicionalmente la enseñanza en escuelas de ingeniería se ha basado en clases teóricas en las que se deducen las ecuaciones y se trata de alcanzar soluciones analíticas, que a menudo conllevan largas y tediosas manipulaciones donde el alumno pierde la noción de los conceptos más importantes. Es lógico que durante la formación del ingeniero este deba conocer el origen de las expresiones, pero aún más importante es que sepa aplicarlas adecuadamente y sea capaz de interpretar los resultados. La resolución de las ecuaciones no debe suponer un impedimento. Actualmente, gracias a programas de cálculo numérico como Matlab u Octave, es posible conseguir soluciones de forma numérica con unas cuantas lineas de código, permitiendo al alumno centrarse en las variables del problema. De esta manera el aprendizaje se vuelve más eficiente y ameno. A esto hay que añadirle que, la implementación de programas para resolver problemas prácticos hace que el alumno aplique y desarrolle los conocimientos adquiridos sobre cálculo numérico, imprescindibles en ingeniería. Los objetivos principales del proyecto son dos. En primer lugar investigar acerca de los posibles métodos que existen para el diseño de perfiles y justificar la necesidad de los mismos. Y segundo, desarrollar un programa de cálculo inspirado en la metodología inversa, bajo las simplificaciones típicas del flujo potencial incompresible. Se evaluarán los resultados y se comparará sus posibilidades respecto a códigos actuales. El documento está estructurado en cuatro capítulos. El primer capitulo introduce los conceptos básicos relacionados con los perfiles aerodinámicos para justificar la importancia de los diseños, y explica los métodos que existen para abordar el problema de diseño. En el segundo se plantean las ecuaciones del problema que relacionan la forma del perfil con la distribución de presiones, dividiendo el problema en espesor y curvatura. El tercero explica la forma de implementar las soluciones anteriormente encontradas mediante tres métodos numéricos distintos. También se especifica la manera de obtener los datos de entrada necesarios. Por último, en el cuarto capítulo se exponen los resultados, se analizan las características de los tres métodos empleados y se determina cual de ellos es el mejor. XIII 1 Fundamentos del diseño inverso Co mo es natural, antes de abordar cualquier proyecto en ingeniería, es necesario realizar una investigación previa sobre el tema del que se va a tratar, familiarizándose con los conceptos y terminología empleada en dicho campo. El propósito de este primer capitulo es hacer comprender al lector sobre la importancia del diseño de perfiles y cuales son sus fundamentos. Se comenzará por exponer una introducción sobre perfiles aerodinámicos en la que se explica: qué son, cómo surgieron y cuáles son sus características que los hace ser elementos imprescindibles. Además se realizará un estudio sobre los efectos y fuerzas presentes en los perfiles que cualquier ingeniero tiene que tener en cuenta a la hora del diseño. Por último, se hablará sobre los distintos métodos que han surgido a lo largo de la historia con el objetivo de conseguir geometrías cada vez más eficientes. Cabe destacar, que será en este apartado donde el lector podrá comprender el concepto de diseño inverso, que a priori puede resultar extraño o confuso. 1.1 Perfiles aerodinámicos 1.1.1 Reseña histórica El primer trabajo serio sobre el desarrollo de las secciones de perfil aerodinámico se inició en la década de 1800. Aunque por aquel entonces ya se sabía que las placas planas producían sustentación al establecer un ángulo de incidencia, algunos intuían que las secciones con curvatura, más parecidas a las alas de las aves, podrían conseguir mayor sustentación y eficiencia. Bajo esta premisa, HF Phillips patentó una serie de formas aerodinámicas (fig. 1.1.a) en 1884 después de probarlos en uno de los primeros túneles de viento de la historia, en el cual la "corriente artificial de aire" se generaba a partir de la inducción por un chorro de vapor en un tronco de madera a modo de conducto. Casi al mismo tiempo, Otto Lilienthal tenía ideas similares, puesto que creía que la clave del éxito del vuelo residía en la curvatura. Después de estudiar cuidadosamente las formas de las alas de las aves, diseñó y experimentó con una serie de perfiles (fig. 1.1.b) en los que varía tanto el radio de borde de ataque, como el espesor y que están recogidos en su libro Bird Flight as the Basis of Aviation, 1894. En 1903, los hermanos Wright emplearon unas secciones muy parecidas a las propuestas por Lilienthal: delgadas y muy curvadas. Posiblemente la razón de esta selección se debiera a que las primeras pruebas sobre perfiles aerodinámicos se realizaron a muy bajo número de Reynolds, donde dichas secciones se comportan mucho mejor que las más gruesas. Esta creencia errónea de que las superficies aerodinámicas eficientes tenían que ser finas y altamente curvadas fue una de las razones por las que algunos de los primeros aviones eran biplanos. Pasaron los años y se desarrollaron muchas otra superficies sustentadoras basándose principalmente en ensayo y error. Algunas de las secciones más exitosas de esa época fueron la Clark Y y Göttingen 398, presentes en la fig. 1.1.c, que fueron tomadas como base para una familia de secciones analizadas por la NACA (Comité Consejero Nacional de Aeronáutica) en la década de 1920. Es a partir de la Segunda Guerra Mundial, cuando en los Estados Unidos la NACA, desarrolla la mayor parte de los perfiles empleados en la actualidad. Todos ellos quedaron recogidos en catálogos, como Summary of Airfoil Data de Abbott y Von 1 21. Fundamentos del diseño inverso (a) H.F. Phillips, 1884 (b) O.Lilienthal, 1894 (c) Desarrollos s.XX Figura 1.1 Evolución histórica de los perfiles aerodinámicos [3]. Doenhoff [ 12 ], donde se define su geometría y sus curvas aerodinámicas. Entre estas secciones destacaban las clasificadas como series de 4 y 5 dígitos. En los años 50, la NACA dejó de lado el diseño de perfiles para centrarse en el estudio de la aerodinámica a gran velocidad, cuyo objetivo era comprender los fenómenos que ocurren durante un vuelo en régimen transónico y supersónico. Entrados ya en los 60, Somers [ 22 ] explica que la NASA (Administración Nacional de la Aeronáutica y del Espacio), sucesora de la NACA, volvió a entrar en el campo de los perfiles con el diseño de perfiles supercríticos por parte de Richard T. Whitcomb. Añade que el conocimiento adquirido durante esos años sirvió para intentar conseguir una nueva serie de flujo turbulento y baja velocidad, mucho más eficiente que los anteriores NACAs. Pese a lograr valores mayores en los coeficiente de sustentación máxima, la resistencia en crucero no fue tan baja como las series de 4 y 5 dígitos. Lejos de rendirse, la NASA siguió investigando, lo que la llevó hacia los perfiles de flujo laminar natural (NLF). Estas secciones se caracterizan por conseguir mantener el flujo laminar durante gran parte de la cuerda ( ≥30% ), tanto en extradós como en intradós, simplemente a base de gradientes favorables de presión. Como se explicará en los siguientes apartados, el mantener flujo laminar sobre el perfil beneficia a la reducción de la resistencia. Por otro lado, Somers [ 22 ] apunta que las investigación sobre NLF fueron fomentadas por la aparición de las estructuras compuestas. Hoy en día los perfiles son diseñados específicamente para cada aplicación, donde encontramos desde hélices, alas y estabilizadores hasta alabes y carenados. De forma similar a los catálogos existentes en aquella época, actualmente se encuentran en Internet multitud de bases de datos donde pueden ser consultadas las características de muchos perfiles utilizados en aviación civil. Un ejemplo de ello está en www.airfoiltools.com, que no sólo permite consultar las gráficas de un perfil sino que además puede mostrar comparativas entre varias secciones. Sin embargo, no toda la información relacionada con los perfiles es pública, como por ejemplo aquellas geometrías asociadas a programas militares o empresas privadas. 1.1 Perfiles aerodinámicos 3 1.1.2 Definición y descripción de un perfil aerodinámico En aeronáutica se denomina perfil aerodinámico, perfil alar, o simplemente perfil, a la forma de la sección transversal de un elemento, que al moverse en el seno de un fluido es capaz de crear una distribución de presiones a su alrededor de manera que genere sustentación. Tal y como explica Anderson [ 5 ], la configuración de esta sección, es una de las consideraciones más importantes a la hora de diseñar superficies sustentadoras como alas, o de otros elementos similares como las palas de hélices, los rotores en helicópteros, los álabes de una turbina y/o compresor y estabilizadores. Según el propósito que se quiera conseguir en el diseño de estos cuerpos, los perfiles pueden ser más finos o gruesos, simétricos o no, curvos o poligonales, e incluso poseer una sección variable a lo largo de la envergadura. Sin embargo, este concepto no se limita solamente a las aeronaves, ya que todo objeto posee un perfil característico. Dependiendo de la relación de aspecto entre la longitud transversal ( `T ) y la longitudinal ( `L ), dichas secciones se dividen entre cuerpos romos y cuerpos fuselados. Los primeros están caracterizados por una relación de aspecto de orden unidad ( `T/`L∼O(1) ), mientras que los fuselados son más estrechos ( `T/`LO(1) ). No obstante, no todos ellos son actos para generar sustentación ( L ) de forma eficiente, puesto que ademas de esta fuerza también se tiene en cuenta la resistencia ( D ) que ofrecen al movimiento. Por este motivo, en Aerodinámica es muy común evaluar la eficiencia de los diseños a través del cociente L/D o de su equivalente con coeficientes adimensionales CL/CD , definiéndose como eficiencia aerodinámica. Gordillo [ 11 ] señala que cuanto mayor sea este valor, menor será el empuje necesario para lograr una fuerza de sustentación igual a la del peso total de la aeronave. Añade también, que para la aviación comercial el valor de CL/CDsuele estar entre 10 y 20, pudiendo llegar a ser superior a 50 en el caso de ultraligeros. En conclusión, un diseño aerodinámico será tanto mejor cuanto mayor sea su eficiencia y esto se logra diseñando superficies que consigan CL∼O(1) y CDO(1) . Esta condición impone de manera indirecta restricciones sobre los tipos de geometría que se deben usar para las superficies sustentadoras. Observando la fig. 1.2.a , se ve que el coeficiente de resistencia CD es de orden unidad para los cuerpos romos y que cae rápidamente al disminuir la relación de aspecto. Además, experimentalmente se demuestra que CL∼O(1) sólo se alcanza para ciertos perfiles fuselados a un ángulo de ataque pequeño como el expuesto en la fig. 1.2.b . Por lo tanto, suponiendo posible conseguir valores de CL∼O(1) , el cumplimiento de la condición CL/CD1 requiere que CD1 , y como ya se ha visto esto sólo se logra mediante cuerpos fuselados. Esa es la razón por la cual en el campo de la aeronáutica se emplean cuerpos esbeltos para los perfiles de las superficies sustentadoras y de control. Por este motivo el trabajo se centrará en el diseño de este tipo de geometrías. (a) Coef. de resistencia en función del número de Reynolds para distintas geometrías (b) Coef. de sustentación y de momento para un NACA 0012 Figura 1.2 Justificación del uso de perfiles fuselados para superficies sustentadoras [11]. 41. Fundamentos del diseño inverso No obstante, lo anterior no excluye la presencia de cuerpos romos en la vida real. Aunque no son de interés para el trabajo desarrollado en este documento, la aerodinámica también los estudia con otros propósitos. El ejemplo más cercano lo encontramos en los automóviles (fig. 1.3), cuyas formas han cambiado para conseguir vehículos más eficientes al disminuir su resistencia al avance. (a) Renault 8 (1962-1971) (b) Seat 1500 (1963–1972) (c) Renault Clio (2015) (d) Seat Leon (2015) Figura 1.3 Evolución de perfiles en automoción. Para comenzar a trabajar con perfiles es preciso definir las características geométricas que los describen. La nomenclatura fue establecida por la NACA en los años 30 cuando comenzó a experimentar con distintas secciones creadas de forma sistemática. La terminología se ilustran en la fig. 1.4, y sus definiciones pueden encontrarse en Anderson [5] y Carmona [13]. Figura 1.4 Nomenclatura de un perfil alar. •Extradós (upper surface) : parte superior del perfil, generalmente asociada a la cara de succión. •Intradós (lower surface) : parte inferior del perfil, generalmente asociada a la cara de presión. •Borde de ataque (leading edge) : punto más adelantado de la sección (B.A). •Borde de salida (trailing edge) : punto más retrasado de la sección (B.S). • Cuerda (chord) : línea recta que une el borde de ataque con el borde de salida. Es una de las dimensiones principales del perfil. • Linea de curvatura media o linea media (mean camber line) : línea formada por los puntos equidistantes entre extradós e intradós. La forma de esta línea es muy importante en las características aerodinámicas del perfil, puesto que fija la curvatura de la sección. Si la linea de curvatura media está por encima de la cuerda como en la fig. 1.4 , se dice que la curvatura es positiva, negativa si cae por debajo y de doble curvatura si en un tramos va por arriba y en otro por debajo. • Radio del borde de ataque (leading-edge radius) : define la forma del borde de ataque como un círculo tangente al extradós e intradós, y cuyo centro se sitúa en la linea tangente en el origen a la linea de curvatura media. Su magnitud define la agudeza del B.A, lo que tiene efectos sobre la entrada en pérdida. • Punto de máximo espesor (maximum thickness) : es el punto sobre la cuerda, cuya distancia entre extradós e intradós medida de forma perpendicular a la cuerda, sea máxima. El espesor máximo es una característica importante, que se suele expresar en porcentaje de la cuerda. El valor varía desde un 3% en los más esbeltos hasta un 18% en los muy gruesos. 1.1 Perfiles aerodinámicos 5 • Curvatura máxima (maximum camber) : es la máxima distancia entre la línea de curvatura media y la cuerda del perfil. Este valor y su posición a lo largo de la cuerda ayudan a definir la forma de la línea de curvatura media. El valor de la ordenada máxima y su posición suelen darse en forma de porcentaje de la cuerda. • Ángulo de ataque (angle of attack) : es el que forma la dirección de la corriente incidente respecto a la cuerda del perfil. (fig. 1.5) Figura 1.5 Definición del ángulo de ataque. 1.1.3 Clasificación de perfiles. Familias Con el fin de estandarizar la nomenclatura de los perfiles y poder identificarlos de forma inequívoca, se desarrollaron las diferentes series. Sin embargo, cada familia identifica sus perfiles de manera particular con un código numérico en el que cada dígito especifica una característica del perfil siguiendo unas reglas de diseño propias de cada serie. Entre las diferentes familias de perfiles se pueden destacar las series NACA, TsAGI, Joukowsky, Clark Y, RAF, Gottingen, Eppler, NLF, GA, Wortmann, etc. Las series NACA son las más empleadas en la actualidad, no limitándose unicamente al diseño de alas ni solamente a la aviación. Las series TsAGI son el homólogo ruso de los perfiles NACA, y por ello son muy utilizados en la industria aeronáutica rusa. Por su parte, las secciones Joukowsky fueron las primeras que se definieron de forma teórica, no obstante actualmente no se usan. Los perfiles Clark Y (EE.UU) y RAF (Reino Unido), son los primeros intentos de definir una familia de perfiles, actualmente se emplean para hélices. Por otra parte, las series Gottingen, Eppler y Wortmann son secciones alemanas; donde la primera familia es de uso general, mientras que la segunda y tercera están especializadas en planeadores. Para terminar, los perfiles NLF y GA son los últimos diseños desarrollados por la NASA, empleándose los primeros en planeadores y los segundos en aviones ligeros. [20] Figura 1.6 Familias de perfiles aerodinámicos. A pesar de la gran variedad de perfiles que existen, no siempre son escogidos. Hoy en día, gracias a los recursos computacionales de los que se dispone, es más fácil y efectivo para los diseñadores crear y optimizar un perfil específico adaptado a una aplicación particular, en lugar de elegir uno ya existente. A continuación, puesto que son los más extendidos en la industria aeronáutica y serán los que se usen en este proyecto, se presenta la nomenclatura de las series NACA. 12 1. Fundamentos del diseño inverso Desde el punto de vista de un vuelo real, estas gráficas son muy importantes puesto que indican a los pilotos el margen de maniobra del que disponen. Durante el vuelo de una aeronave, la sustentación está determinada en cada momento por la condición de vuelo. Suponiendo un valor fijo de la sustentación l , y teniendo en mente su relación con el coeficiente Cl , ec. (1.5) , es posible ver que los pilotos pueden realizar dos acciones para alcanzar dicho valor. Bien pueden variar el ángulo de ataque α para conseguir el valor Cl(α) necesario con una cierta velocidad de vuelo U∞ , o bien variar la velocidad U∞ sin cambiar el ángulo de ataque α . Luego a altas velocidades, el ángulo de ataque α será pequeño (poco Cl ), mientras que a bajas velocidades ocurre lo contrario, se incrementará α para poder mantener el mismo nivel de sustentación. No obstante, se ha visto que existe un Cl,max asociado a un αmax , lo que implica que existe un valor mínimo requerido de la velocidad para que el vuelo sea posible. Dicho valor recibe el nombre de velocidad de entrada en perdida, ya que si la velocidad es inferior a este valor, la aeronave no es capaz de generar la sustentación necesaria y caería. l=1 2ρ∞U2 ∞c Cl(α)(1.5) Centrándonos ahora en la resistencia de los perfiles, se ha dicho que tiene dos componentes: una debida a presión y la otra a fricción. Ambas están presentes en el coeficiente Cd . Otra de las curvas más habituales que se usan para estudiar y comparar los perfiles entre sí es la polar del perfil, donde se representa Cl frente a Cd . Es preciso señalar que la fig. 1.13 muestra que los valores de la resistencia siempre son positivos, es decir, aunque el perfil sustente de forma negativa, esta fuerza siempre se opone al movimiento. Figura 1.13 Polar en función del espesor t [11]. Respecto a la resistencia de presión, también es conocida como resistencia de forma, ya que depende en gran medida de la forma de la sección del cuerpo. Carmona [ 13 ] señala que esta resistencia está relacionada con la forma de la estela (más o menos gruesa) que deja el objeto tras su paso, puesto que está asociada al punto de desprendimiento de la capa límite y por ende a los esfuerzos viscosos. Cuanto antes se desprende la corriente, mayor es la estela y mayor es la resistencia que se obtiene. Por lo anteriormente expuesto, también puede denominarse resistencia de estela. Por otra parte, la resistencia de fricción está vinculada al rozamiento que existe entre el fluido y el sólido en la capa límite, debido a los esfuerzos viscosos. Depende en gran medida del tipo predominante de capa límite que haya sobre el perfil. Normalmente su valor suele ser inferior a la debida por la presión, pero no tanto como para ser despreciada en los diseños. Visto lo visto, el comportamiento de la resistencia pasa inevitablemente por el estudio de la capa límite y la viscosidad. 1.2 Fuerzas aerodinámicas 13 1.2.2 Capa límite e influencia de la viscosidad Llegados a este punto es necesario explicar que la corriente alrededor de un perfil, es habitualmente descompuesta en dos regiones: exterior e interior. Idea planteada por Prandtl en 1904. La diferencia entre estas dos zonas está en la influencia de los esfuerzos viscosos sobre la corriente, siendo despreciables en el caso exterior donde puede considerarse que el fluido es ideal, pero no en el interior por muy pequeño que sea el valor de la viscosidad del fluido. Por este motivo a la región exterior también se la conoce como región no viscosa, mientras que a la interior se la denomina región viscosa o capa límite. Perez-Saborid [ 6 ] aclara que la importancia de los efectos viscosos dentro de la capa límite radica en que los gradientes de velocidad en la dirección normal a la superficie son muy acusados al darse en una capa muy delgada, y no por el valor del coeficiente de viscosidad. La diferencia entre estas dos zonas se comprueba experimentalmente observando detenidamente los perfiles de velocidad sobre la sección. Figura 1.14 Perfiles de velocidad dentro de la capa límite . (a) Placa plana (b) Borde de ataque Figura 1.15 Capa límite sobre secciones . Estudiando el perfil de velocidades de la fig. 1.14, se aprecia que para un punto x sobre la sección, lejos de la pared la velocidad es prácticamente constante, lo que se corresponde con la corriente libre. Sin embargo, al acercarnos a la pared, la velocidad comienza a disminuir rápidamente hasta hacerse cero en la pared. Esto es así puesto que, tal y como indica Gordillo [ 11 ], los esfuerzos viscosos que actúan sobre las partículas fluidas en las inmediaciones de la superficie del perfil, son los únicos capaces de decelerar el fluido en la dirección tangencial a la superficie del sólido para así hacer cumplir la condición de no deslizamiento del fluido sobre la pared del mismo. Tras esta explicación, se deduce que el espesor δ de la capa límite es la distancia entre la pared del perfil y el punto donde el valor de la velocidad sea igual al 99 % de la corriente libre. Se trata de una región en la que se disipa energía, y por tanto cuanto mayor sea su extensión sobre el perfil, afectará a más capas del fluido, aumentando su espesor y disipando más energía, lo que implica una mayor resistencia de fricción. 14 1. Fundamentos del diseño inverso Al igual que ocurre en cualquier flujo, dentro de la capa límite la corriente puede ser laminar oturbulenta. En el primer caso, las lineas de corriente fluyen de forma paralela, mientras que en el caso turbulento estas se entrelazan y la corriente se vuelve desordenada. Carmona [ 13 ] explica que generalmente el desarrollo de la capa límite empieza por ser laminar cerca del borde de ataque. A medida que esta avanza por la superficie, disipa cada vez más energía haciendo que su espesor crezca lentamente, hasta que llegada una cierta distancia del borde de ataque, la capa límite comienza a verse afectada por unas perturbaciones de tipo ondulatorio que provocan un aumento repentino de su espesor y la destrucción del flujo laminar, pasando a ser turbulento. Este desarrollo queda ilustrado en la fig. 1.16. La región en la que ocurre el paso de laminar a turbulenta se denomina región de transición, aunque al ser una zona tan pequeña suele hablarse de punto de transición. Figura 1.16 Evolución de la capa límite [13]. Figura 1.17 Capa límite: laminar y turbulenta [8] . En una capa límite turbulenta las partículas se entremezclan. Las que están alejadas de la pared pasan por zonas cercanas a esta, y viceversa, las próximas a la pared fluyen alejándose de ella. El resultado global es un aumento de la velocidad dentro de la capa límite respecto a la laminar, y esto es debido a que las partículas exteriores de la capa límite absorben energía de la corriente libre y la comunican a las zonas próximas a la pared, consiguiendo una mayor velocidad que para el caso laminar. Comparando los perfiles de velocidad característicos, fig. 1.18, se observa que cerca de la pared se cumple dV dy lam <dV dy tur , lo que implica que las fuerzas de fricción son mayores para una capa turbulenta que para una laminar. Así que de cara al diseño para reducir la resistencia de fricción, conviene tener capa laminar sobre perfil o en su defecto, intentar retrasar la transición lo máximo posible. Para ello el diseñador debe tener cuidado con la rugosidad de la superficie, las turbulencias de la corriente exterior, y sobretodo controlar los gradientes adversos de presiones sobre la sección. 1.2 Fuerzas aerodinámicas 15 (a) Comparación de perfiles [13] (b) Perfil laminar [8] Figura 1.18 Perfiles de velocidad en la capa límite laminar y turbulenta . Tras lo explicado, uno puede deducir sin dudar, que existe una relación entre la velocidad de la corriente libre y la de la capa límite. No obstante, este razonamiento puede llegar más lejos, puesto que como ya se ha expuesto, la velocidad exterior depende de la distribución de presiones sobre el perfil. Por lo tanto, la capa límite también se ve influenciada por las presiones. En efecto, si el gradiente de presiones es favorable (se fluye hacia zonas de menor presión) la corriente exterior se acelera, y tanto la velocidad exterior como la de la capa límite tienen el mismo sentido. En esta circunstancia se dice que la capa límite permanece adherida, porque la corriente exterior fluye alrededor de una región viscosa que coincide prácticamente con el sólido[ 11 ]. Por el contrario, si el gradiente de presiones es adverso (se fluye hacia zonas con mayor presión), el flujo libre se ve frenado, siendo mayor esta deceleración en la capa límite debido a los efectos viscosos, lo que puede provocar la inversión del flujo en la región viscosa. Este fenómeno se conoce como desprendimiento de la capa límite, pues los efectos viscosos ya no quedan restringidos a una pequeña zona, provocando una estela turbulenta tras el objeto. Figura 1.19 Capa límite sobre perfiles [11] . 16 1. Fundamentos del diseño inverso Trasladados estos conceptos a los perfiles alares, el comportamiento de la capa límite queda ilustrado en la fig. 1.19. Se muestra que de forma genérica la capa permanece adherida en todo el intradós y parte del extradós donde el gradiente es favorable y la corriente se acelera. Además, se ve como tras el pico de succión en el extradós cerca del borde de ataque, el gradiente se vuelve adverso decelerando la corriente hasta llegar al desprendimiento. Esto genera una estela turbulenta cuya presión permanece constante y próxima al valor de la presión en el punto de desprendimiento, pd<p∞ , lo que origina la fuerza de resistencia de forma o estela causada por la diferencia de presiones. Debido a la existencia de viscosidad siempre habrá desprendimiento de la capa límite sobre los perfiles, aunque este tenderá a suceder cerca del borde de salida para ángulos pequeños. Si por el contrario, el punto de desprendimiento se adelanta hacia el borde de ataque, ocurre el fenómeno de entrada en pérdida, pues se rompe el flujo laminar sobre el perfil y no es capaz de generar sustentación, lo que ademas conlleva un gran incremento de la resistencia de presión. El punto de desprendimiento es un punto fijo de cada sección y se caracteriza por sus valor máximo Cl,max asociado a un αmax. Carmona [ 13 ] muestra que existen tres tipos de entrada en pérdida y que estas depende de la forma del perfil. El primer tipo consiste en una entrada en pérdida por el borde de salida, típica de las secciones gruesa ( espesor >14% ) con borde de ataque muy redondeado. La capa límite es turbulenta, aumentando su espesor poco a poco con el ángulo de ataque. Cuando este vale aproximadamente 10 ° , la capa empieza a desprenderse del B.S., adelantándose el punto de separación hacia el B.A. conforme el ángulo crece. La pérdida de sustentación es gradual y suave incluso una vez alcanzado el valor máximo. El segundo tipo es una entrada en pérdida por el borde de ataque, característica de perfiles de espesor medio ( ≈10% ). A un ángulo pequeño la capa límite se desprende cerca del B.A., por lo que la capa todavía es laminar. Tras separarse, se vuelve inestable y se convierte en turbulenta, adhiriéndose de nuevo al perfil y dejando encerrada una pequeña burbuja de recirculación. Al seguir aumentando el ángulo, la burbuja se hace más pequeña hasta que se desprende cerca del B.A. provocando una pérdida brusca. Por último, el tercer tipo es propio de secciones muy finas ( ≈5% ). A ciertos ángulos de ataque se forma una burbuja que provoca la disminución de la pendiente de sustentación, pero no la entrada en pérdida. Conforme crece el ángulo, la burbuja se extiende hasta el B.S., momento en el que se alcanza el Cl,max. La pérdida se produce de forma gradual. Figura 1.20 Formas de entrada en pérdida [13] . Figura 1.21 Perfil en pérdida [8] . 1.3 Métodos de diseño 17 Al diseñar habrá que tener en mente la entrada en pérdida, por ser un fenómeno no deseado, ya no desde el punto de vista del incremento de la resistencia de forma, si no por seguridad, porque podría ser la causa de un accidente. Para retrasar el desprendimiento de la capa límite, lo más fácil e inmediato es aumentar la velocidad de la capa límite, dado que se origina por la pérdida de esta. Una de las técnicas para conseguir este aumento consiste en posicionar generadores de vórtices en la superficie del perfil, de tal manera que la capa límite pase a ser turbulenta. Se ha comentado que este tipo de capa posee mayor energía, por lo que es capaz de permanecer adherida durante más tiempo (fig. 1.17), aunque también supondría una mayor resistencia de fricción. Otra forma empleada en los diseños actuales, es la inyección de aire proveniente desde el intradós hasta el extradós, lo que supone un aporte de energía a la capa límite para seguir adherida. Esto último suele combinarse con los dispositivos hipersustentadores, slaps y flaps, que amplían el rango operativo de los perfiles a bajas velocidades. Como se ha podido comprobar, desde el punto de vista del diseño es fundamental conocer el comportamiento de la capa límite, ya que pese a ser una región muy estrecha cerca del perfil, influye en gran medida en las características aerodinámicas de los perfiles: •Resistencia de fricción •Resistencia de presión •Valor máximo del coeficiente de sustentación Cl,max •Forma de entrada en pérdida 1.3 Métodos de diseño Al comienzo de este documento, se explicó que los primeros diseños fueron fruto de la observación de la naturaleza, generalmente de las alas de las aves, y que se perfeccionaron a base de ensayo y error en rudimentarios túneles de viento. Pronto se dieron cuenta de la importancia de la forma de las secciones para generar las fuerzas aerodinámicas, siendo la piedra angular de la investigación aeronáutica [22]. El objetivo que persigue todo diseñador es conseguir que su diseño cumpla una serie de condiciones y/o que posea una serie de características según la aplicación para la que esté destinada la sección. La idea básica que está detrás de todos los diseños aerodinámicos modernos es: "las fuerzas aerodinámicas y las características de la capa límite son resultado de la distribución de presiones que se crea alrededor del perfil, y que a su vez depende de la forma del propio perfil". Esto implica que el diseñador debe ser capaz de trasladar los requisitos de su diseño a una distribución de presiones que le permita alcanzar sus objetivos. El problema de diseño se reduce entonces a encontrar la forma del perfil asociada a dicha distribución. Durante los primeros años se empleó de forma exitosa el uso de los catálogos. Estos disponían de información acerca de los datos de las pruebas, lo que hacia posible seleccionar aquel perfil que más se acercara a las necesidades según la aplicación. Algunos de estos catálogos son: Summary of Airfoil Data de Abbott y Von Doenhoff [ 12 ], Profilkatalog Stuttgarter de Althaus y Wortmann, Low Reynolds Number Airfoil catalog de Althaus y Airfoils at Low Speeds de Selig. [2] Con el paso del tiempo las aplicaciones se han diversificado, y las propiedades requeridas se han vuelto cada vez más sofisticadas, quedándose fuera del rango cubierto por los catálogos. Aquí es donde surgen multitud de métodos para tratar de determinar los perfiles más adecuados. No obstante, todos los procesos pueden ser agrupados en dos filosofías: métodos directos o métodos inversos. (fig. 1.22) La necesidad del desarrollo de estos métodos de diseño, está justificada por tres motivos: •Diseño de perfiles fuera del alcance de los catálogos. •Diseño específico para cada aplicación. •Ahorro económico en concepto de pruebas y ensayos. 18 1. Fundamentos del diseño inverso Figura 1.22 Clasificación de los procesos de diseño. Hoy en día, gracias al desarrollo de métodos numéricos y su aplicación al campo de la mecánica de fluidos, existen en el mercado los códigos CFD (Computational Fluid Dynamics). Estos programas permiten simular y estudiar la evolución del campo fluido en infinidad de situaciones. Dentro de las diferentes aplicaciones de los CFD, se encuentra una rama especializada en aerodinámica y diseño de perfiles. La gran mayoría utilizan el método directo simulando ser túneles de viento virtuales donde probar los diseños, dejando los métodos inversos sólo para los códigos más sofisticados. No hay que olvidar que los CFD son una herramienta más para el diseñador, siendo este el que debe tomar las decisiones sobre el diseño. Además, a pesar de las buenas aproximaciones ofrecidas por estos programas, siempre será necesario realizar pruebas en túneles de viento para corroborar con los resultados experimentales. En el caso de este estudio, se ha escogido el programa XFLR5 del que se hablará más adelante. 1.3.1 Diseño directo En un proceso directo, se toma como punto de partida la forma del perfil aerodinámico. Especificada la geometría de la sección, se procede al cálculo de la distribución de presiones a su alrededor, lo que permitirá conocer sus características aerodinámicas tales como: desarrollo de la capa límite, pendiente de sustentación, polar del perfil, ángulo de ataque máximo, posición del centro de presiones, etc. Este proceso queda reflejado en la fig. 1.23. Con toda esta información, se evalúa la forma dada y a continuación se modifica la geometría para tratar de mejorar el o los rendimientos aerodinámicos según los objetivos marcados. Se trata por tanto de un proceso iterativo de optimización hasta conseguir el perfil con el rendimiento deseado. La optimización siempre está ligada a la minimización de una función objetivo relacionada con las propiedades del perfil controladas por estos rendimientos.[16] Como puede intuirse, es un proceso largo y costoso computacionalmente, pues son necesarias muchas iteraciones para obtener las características deseadas. Además, requiere de experiencia previa, ya que el diseñador tiene que intuir como cambiará la distribución de presiones según las modificaciones que se realicen en la geometría. Uno de los primeros en trabajar sobre este tipo de metodología fue Theodorsen [ 23 ]. En los años 30, presentó un método teórico para calcular la distribución de velocidades a partir de una geometría arbitraria, sin considerar efectos viscosos ni compresibilidad. La idea era emplear la transformación conforme, convirtiendo la geometría del perfil en un círculo, y determinar sobre este tanto la velocidad como la presión. Para finalizar, invertía el proceso y obtenía estas distribuciones asociadas al perfil. Es un método válido para cualquier perfil aerodinámico, incluso para los de gran espesor. La validez de los resultados fruto de este desarrollo analítico quedan patentes en la fig. 1.24. 1.3 Métodos de diseño 19 Figura 1.23 Proceso de diseño directo [16]. (a) Clark Y α=9°(b) Clark Y α=−1° Figura 1.24 Theodorsen, Comparación distribución de presiones experimental y teórica [23] . 20 1. Fundamentos del diseño inverso 1.3.2 Diseño inverso Al contrario que en el diseño directo, los métodos inversos comienzan con la definición de las especificaciones aerodinámicas y tratan de determinar la forma de la sección que reproduzca dichas características. El problema se ha invertido, pues las propiedades del perfil que antes eran calculadas como salida del problema, ahora son la entrada para calcular la geometría. Lo más común en todos estos métodos inversos, es que el diseñador defina como función objetivo la distribución de presiones o velocidades sobre el perfil que generen las características aerodinámicas deseadas. Está filosofía de diseño queda reflejada en la fig. 1.25. Poseen una base teórica sólida, pero es necesario que la distribución de presiones definida, tenga sentido físico. Jepson[ 16 ] apunta que la motivación detrás de estos métodos era aprovechar las relaciones conocidas entre la distribución de presiones o velocidades, y otras propiedades aerodinámicas como sustentación, resistencia, momento de cabeceo y capa límite. Como ya se explicó, de la diferencia de presiones entre extradós e intradós se obtiene la sustentación a cualquier ángulo de ataque. Además, los gradientes de presiones determinan el desarrollo de la capa límite que influye tanto en la resistencia como en la entrada en pérdida que define el αmax y el Cl,max . Se puede concluir, que el rendimiento total del perfil está indirectamente controlado mediante la presión. Esto sugiere que, simplemente ajustando las distribuciones de presiones según los requisitos para el diseño, se pueden conseguir perfiles óptimos para cada aplicación. Figura 1.25 Proceso de diseño inverso [16]. Por ejemplo, suponiendo el caso en que el diseñador esté interesado en perfiles para planeadores, deberá buscar secciones capaces de generar gran sustentación y poca resistencia, a pequeños ángulos de ataque, para que la aeronave se mantenga en vuelo el mayor tiempo posible. Atendiendo a estas necesidades, para lograr alta sustentación el perfil deberá generar un pico de succión importante que provoque una gran diferencia de presiones entre extradós e intradós, teniendo cuidado con el gradiente adverso de presiones para no entrar en pérdida; mientras que si se quiere disminuir la resistencia, es fundamental controlar que la capa límite sea laminar. Bajo estas condiciones, el diseñador podría plantear las curvas Cp mostradas en la fig. 1.26.a como las deseadas para un perfil posicionado a α=3° . Cada una de las curvas representa una de las características buscadas: perfiles de alta sustentación en rojo, y perfiles de flujo laminar en azul. Se observa claramente como ambas condiciones parecen contrarias entre sí, ya que la diferencia de presiones en los perfiles laminares es inferior a los de alta sustentación, no consiguen mucha sustentación a pequeños ángulos de ataque puesto que se centran en controlar el gradiente adverso. Invirtiendo el problema, los perfiles obtenidos a partir de esas curvas serían NLF1015 y FX 60-100, representados en fig. 1.26.b . Así que, en definitiva, el diseñador debería llegar a una solución de compromiso entre las dos. 1.3 Métodos de diseño 21 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 Cp x Flujo Laminar Alta sustentación (a) Curvas Cpobjetivo 0 0.2 0.4 0.6 0.8 1 −0.25 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 0.25 z x NLF 1015 FX 60−100 (b) Perfiles asociados Figura 1.26 Ejemplo de diseño . Ahora bien, lo realmente complicado de esta metodología está en el hecho de encontrar el procedimiento adecuado para relacionar las presiones con la forma del perfil. Atendiendo a la forma de resolución del problema inverso, se diferencian dos tipos de procesos: directos e iterativos. Acoplados o Directos: emplean una formulación alternativa del problema directo. Hacen que las variables geométricas aparezcan como variables dependientes (explícita o implícitamente). Teóricamente son capaces de encontrar la geometría como función del campo fluido en un único paso [4]. Desacoplados o Iterativos: basados en la actualización iterativa de la forma del perfil hasta conseguir la distribución de presiones deseada. Resuelven el problema directo sobre el perfil y comparan la distribución de presiones o los rendimientos obtenidos con los deseados, para a continuación, conforme a los criterios implementados, actualizar la forma del perfil hasta lograr los objetivos marcados [4]. Los métodos inversos son una herramienta que mejoran significativamente la capacidad de diseño de superficies sustentadoras. Sin embargo, las primeras teorías llevadas acabo por Mangler y Lighthill, tenía un punto débil, ya que sólo eran capaces de calcular los perfiles para un único punto de operación definido por la distribución de velocidades a un ángulo de ataque para flujo no viscoso. Para evaluar el rendimiento de la sección fuera de este punto era necesario recurrir a un análisis pos-diseño. Podría pensarse que es suficiente con un único punto y optimizar el rendimiento del perfil para dicho punto, haciendo que el perfil trabajase siempre en ese punto. Sin embargo, desde el punto de vista de las superficies sustentadoras, los perfiles están sometidos a diferentes condiciones determinadas por los segmentos de vuelo y es preciso que tengan un buen rendimiento dentro de un rango operativo, no exclusivamente en un punto. Nace así la imperiosa necesidad de perfeccionar la metodología hacia procesos de diseño inverso multipuntuales. Código de Eppler En 1957, Eppler publicó uno de los primeros trabajos prácticos realizados sobre esta materia. Somers [ 22 ] explica que el código desarrollado por Eppler contenía tres funcionalidades: un método de trazado conformal para el diseño de perfiles a partir de una distribución de velocidades especificada, un método de paneles para analizar el flujo potencial sobre los perfiles, y un método integral para la capa límite. Con este programa era posible tanto diseñar perfiles aerodinámicos con las características de capa límite especificadas, como analizar los perfiles cuya geometría fuera conocida. Representó un gran salto en el campo de diseño de perfiles, y por ello comenzó a ser utilizado por la NASA a partir de 1975. El método de transformación conforme usado por Eppler es único porque permitía definir la distribución de velocidades a lo largo de la sección para diferentes ángulos de ataque. Es una característica muy importante, pues gracias a ella era posible diseñar secciones cuyo comportamiento se ajusta al definido en distintos puntos de operación (proceso multipuntual). 28 2. Planteamiento del problema Se concluye que en primera aproximación se puede despreciar el término de esfuerzos viscosos para el flujo exterior, no siendo así dentro de la capa límite. Se trata de un flujo potencial. • N º Froude : surge al comparar el término de fuerzas másicas ρ−→ fm con la aceleración convectiva, dando lugar a la relación Fr =U2 ∞/(gc) . Aplicando los ordenes de magnitud vistos anteriormente se deduce que Fr 1 , lo que implica de nuevo , que en primera aproximación, se pueden despreciar las fuerzas másicas. • N º Strouhal : contrasta la variación de la velocidad debido a dos fenómenos: cambio de posición ρ−→ v·∇−→ v y cambio de las condiciones de contorno ρ∂−→ v ∂t . Se trata de estimar si las partículas del campo fluido permanecen el tiempo suficiente como para experimentar una variación de la velocidad debida al tiempo. En el caso abordado St =c/(U∞tc)1 , lo que indica que el flujo es estacionario como ya se declaró. Aplicando estas deducciones, se reduce la complejidad de las ecuaciones de continuidad y cantidad de movimiento, ya que se desprecian los términos dependientes del tiempo, las fuerzas másicas y los esfuerzos viscosos. ∇·(ρ−→ v) = 0(2.4) ρ−→ v·∇−→ v=−∇p(2.5) Por otra parte, el flujo que se forma alrededor del perfil puede entenderse como el resultado de las perturbaciones sobre la corriente incidente debido a la presencia del perfil en el campo fluido. De esta manera, las magnitudes fluidas pueden ser modeladas mediante los valores de la corriente incidente ( −→ x→∞ ) más las pequeñas perturbaciones. −→ v(x,z) = −→ U∞+−→ v0=U∞−→ i+u0−→ i+w0−→ k p(x,z) = p∞+p0 ρ(x,z) = ρ∞+ρ0 donde se cumple que: |v0|U∞;|p0| p∞;|ρ0|ρ∞ |−→ x|→∞:−→ v0=0 ; p0=0 ; ρ0=0 Teniendo en cuenta esta consideración de pequeñas perturbaciones provocadas por perfiles esbeltos, se linealiza el sistema formado por las ecuaciones (2.3) - (2.5) . El resultado obtenido es un sistema de ecuaciones (2.6) - (2.8) cuyas incógnitas pasan a ser las perturbaciones en velocidad, presión y densidad. Además, tras esta manipulación, se ha conseguido una relación más directa entre las magnitudes, lo que permite sustituir unas expresiones en otras para llegar a la conexión buscada en este capítulo. ∇·−→ v0=−U∞ ρ∞ ∂ρ0 ∂x(2.6) ρ∞U∞ ∂−→ v0 ∂x=−∇p0(2.7) ρ0=p0 a2 ∞ (2.8) 2.1 Formulación de las ecuaciones 29 Antes de sustituir y obtener la expresión final, se debe proyectar la ecuación de cantidad de movimiento linealizada sobre la componente −→ i . Al hacer esto se consigue la relación directa entre la velocidad y la presión (2.9) necesaria para poder sustituir en el sistema. Al mismo tiempo, permite conocer una expresión del coeficiente de presiones en función de la velocidad de perturbación (2.10) que será muy práctica en este estudio. p0=−ρ∞U∞u0(2.9) Cp=p−p∞ 1 2ρ∞U2 ∞ =p0 1 2ρ∞U2 ∞ =−2u0 U∞ (2.10) Usando la definición de la divergencia ( ∇·−→ v0=∂u0 ∂x+∂w0 ∂z ) para el caso bidimensional y tras sustituir las expresiones anteriores, se llega a la ecuación buscada (2.11) . Se trata de una igualdad cuya solución proporciona el campo de velocidades alrededor de un perfil. Según el desarrollo seguido, la expresión es válida tanto para fluidos compresibles como incompresibles, ya sea en régimen subsónico o supersónico; el único requisito es que las perturbaciones sean pequeñas. (1−M2 ∞)∂u0 ∂x+∂w0 ∂z=0(2.11) Para terminar de cerrar el problema queda definir las condiciones de contorno que debe cumplir el campo de velocidades alrededor de la sección. La primera ya fue incluida con la consideración de pequeñas perturbaciones, y es que lejos del objeto las perturbaciones son nulas (2.12) . La siguiente es conocida como condición de impenetrabilidad, según la cual la velocidad del flujo sobre el perfil es tangente a la superficie (2.13) . De esta forma se consigue relacionar la distribución de velocidades con la forma del perfil. Al tratarse de secciones esbeltas, esta condición puede evaluarse en primera aproximación en z→0 . Por último, debe cumplirse la condición de Kutta-Joukowski según la cual la solución real del problema es aquella en la que la corriente no rebordea el borde de salida. Esta condición junto con el hecho de que las presiones en el borde de salida son iguales tanto en intradós como en extradós, implica que las velocidades en este punto sean las mismas por ambas caras (2.14) y (2.15) . Tanto la demostración de estas condiciones como una explicación más extensa de las mismas son ofrecidas por Katz [15] y Gordillo [11]. −→ v0→0|−→ x|→∞(2.12) −→ v·−→ ns=0 ; w0(x,z→0±) = U∞ dz± p dx ∀x∈Perfil (2.13) u0(x,z→0+) = u0(x,z→0−)x∈BS Cond. de KJ (2.14) w0(x,z→0+) = w0(x,z→0−)x∈BS Cond. de KJ (2.15) Llegados a este punto, el problema puede ser expresado de forma más sencilla si se expresa en función de un potencial. En este caso puesto que se está considerando el flujo alejado de la superficie del perfil, las partículas fluidas no rotan al estar sometidas a un campo de presiones uniforme. Esto indica que la vorticidad es cero ( ∇×−→ v=−→ ω=0 ), lo que implica que la velocidad proviene de un potencial ( −→ v=∇Φ ). Por otra parte, si además se incluye la simplificación de considerar régimen subsónico incompresible ( M2 ∞1 ), entonces el problema deriva en la resolución de la ecuación de Laplace. Gordillo [ 11 ] explica que la aproximación casi incompresible da buenos resultados cuando el número de Mach satisface M∞<0.3 . También afirma que a diferencia de la ecuación (2.1) la de Laplace es lineal, lo que simplifica enormemente la resolución del flujo exterior a la capa límite. 30 2. Planteamiento del problema ∇2Φ0=0(2.16) −→ v0=∇Φ0→0|−→ x|→∞(2.17) −→ v·−→ ns=0 ; ∂Φ0 ∂z(x,z→0±) = U∞ dz± p dx ∀x∈Perfil (2.18) ∂Φ0 ∂x(x,z→0+) = ∂Φ0 ∂x(x,z→0−)x∈BS Cond. de KJ (2.19) ∂Φ0 ∂z(x,z→0+) = ∂Φ0 ∂z(x,z→0−)x∈BS Cond. de KJ (2.20) 2.2 Resolución teórica Lo primero a tener en cuenta para hallar la solución del problema planteado es que gracias a que la ecuación de Laplace es lineal, la solución puede construirse mediante superposición de soluciones elementales Φi . Históricamente el problema se ha estudiado separándolo en un problema de espesor o simétrico, y un problema de curvatura o antisimétrico (2.22) . Esta metodología se apoya en el hecho de que el contorno de los perfiles zp(x) está determinado a través de la función de espesor ze(x) , y la línea de curvatura zc(x) . En el caso de los perfiles NACA de 4 dígitos, estas expresiones están basadas en el valor del máximo espesor t , máxima curvatura hy posición de la máxima curvatura p. (a) P. Completo (b) P. Espesor (c) P. Curvatura Figura 2.1 División del problema completo [11]. zp(x,t,h,p) = zc(x,h,p)±ze(x,t)(2.21) Desde otro punto de vista este fenómeno físico se interpreta a través de la condición de impenetrabilidad. Si no existiera el perfil, las lineas de corriente serían uniformes. Al introducir la sección, las lineas de corriente se perturban para cumplir la condición de impenetrabilidad, apareciendo una velocidad vertical w0 . En conclusión, el origen de esta velocidad vertical tiene dos componentes, una debida al espesor del perfil y otra debido a su curvatura (2.23). Φ0=Φ0 espesor +Φ0 curvatura (2.22) ∂Φ0 espesor ∂z(x,z→0±)+ ∂Φ0 curvatura ∂z(x,z→0±) = U∞(dzc dx ±dze dx )(2.23) 2.2 Resolución teórica 31 2.2.1 Problema de espesor o simétrico Dada las características del problema de espesor, Kantz [ 15 ] y Gordillo [ 11 ] usan una solución elemental de la ecuación de Laplace cuyo campo de velocidades de perturbación generado sea simétrico. De está manera la solución para este problema consiste en una distribución continua de fuentes situadas a lo largo del eje x (fig. 2.2). La intensidad de estas fuentes vo(xo) debe ser calculada para cumplir con la condición de impenetrabilidad, obteniendo así la relación buscada entre la velocidad y el espesor del perfil. Figura 2.2 Distribución continua de fuentes [15]. La velocidad en cualquier punto del dominio fluido viene determinada por la integración de la distribución de las fuentes (2.24). Se distinguen las dos componentes: velocidad horizontal u0 sim−→ iy vertical w0 sim−→ k. −→ v0sim(x,z) = 1 2πZc 0 vo(xo)(x−xo) (x−xo)2+z2dxo−→ i+Zc 0 vo(xo)z (x−xo)2+z2dxo−→ k(2.24) En primer lugar, se calcula el valor de la velocidad vertical w0 sim(x,z→0) , quedando en función de la intensidad de la fuente vo(x) . Al relacionar este hecho con la condición de impenetrabilidad, da lugar a la expresión (2.25) que proporciona el valor de la intensidad de cada fuente [ 15 ]. Al sustituir esta condición en la integral que proporciona la componente horizontal de velocidad, se obtiene una expresión que la relaciona de forma directa con el espesor del perfil (2.26). w0 sim(x,z→0±) = ±vo(x) 2=±U∞ dze dx (x);vo(x) = 2U∞ dze dx (x)(2.25) u0 sim(x,z→0) = U∞ πZc 0 dze dx (xo) (x−xo)dxo(2.26) A pesar de encontrar la función buscada, esta está escrita en términos de velocidades. Sin embargo el interés de este estudio reside en la distribución de presiones que se forma alrededor de las secciones. Como ya se ha visto anteriormente ambas magnitudes están vinculadas entre sí a través de la igualdad (2.9) o (2.10) . Reescribiendo la ecuación anterior, se llega a la expresión buscada para el problema simétrico (2.27). Cp,sim(x) = −2 πZc 0 dze dx (xo) (x−xo)dxo(2.27) Se trata de una ecuación integral cuyo dato de entrada es la distribución de presiones simétrica Cp,sim . Su resolución proporciona los valores de la pendiente de la función de espesor dze dx (x) que serán integrados a lo largo del perfil para calcular el espesor ze(x). 32 2. Planteamiento del problema 2.2.2 Problema de curvatura o antisimétrico Siguiendo la misma estructura, para el problema de curvatura se debe utilizar una solución elemental Φi cuyo campo de velocidades de perturbación sea antisimétrico. En esta ocasión la distribucion continua sobre el eje x es de torbellinos (fig. 2.3). De igual manera que en el anterior apartado, la intensidad de estos torbellinos γo(xo)debe ser la adecuada para cumplir la condición de impenetrabilidad. Figura 2.3 Distribución continua de torbellinos [15]. La velocidad en cualquier punto del dominio fluido viene determinada por la integración de la distribución de los torbellinos (2.28) . También se distinguen las dos componentes: velocidad horizontal u0 anti−→ i y vertical w0 anti−→ k. −→ v0anti(x,z) = 1 2πZc 0 γo(xo)z (x−xo)2+z2dxo−→ i+Zc 0−γo(xo)(x−xo) (x−xo)2+z2dxo−→ k(2.28) Esta vez se comienza por calcular el valor de la velocidad horizontal u0 anti(x,z→0) , obteniendo una relación para el valor de la intensidad de cada torbellino (2.29) . Sustituyendo esta expresión en la integral que proporciona la componente vertical w0 anti de la velocidad, e imponiendo la condición de impenetrabilidad, se consigue análogamente una expresión que relaciona directamente la curvatura del perfil con la velocidad de perturbación horizontal (2.30). u0 anti(x,z→0±) = ±γo(x) 2;γo(x) = 2u0 anti(x)(2.29) w0 anti(x,z→0) = U∞ dzc dx (x) = −1 πZc 0 u0 anti(xo) (x−xo)dxo(2.30) Al igual que antes la ecuación está escrita en términos de velocidad, pero utilizando las mismas expresiones anteriores puede ser reescrita en función de la distribución de presiones. dzc dx (x) = 1 2πZc 0 Cp,anti(xo) (x−xo)dxo(2.31) Se trata de una ecuación cuyo dato de entrada es la distribución de presiones antisimétrica Cp,anti . A diferencia del problema simétrico, ahora no se trata de resolver una ecuación, sino en calcular una integral. Esta integral proporciona directamente los valores de la pendiente de la función de curvatura dzc dx (x) que serán integrados a lo largo del perfil para calcular la curvatura zc(x). 2.2 Resolución teórica 33 2.2.3 Distribución de presiones Por último, es necesario dedicar un pequeño apartado a la distribución de presiones, ya que las expresiones encontradas dependen de estos valores. Habitualmente al estudiar el campo de presiones alrededor de un perfil se distingue la presión del extradós Cp,E y la del intradós Cp,I , ya que la diferencia entra ambas proporciona la sustentación. Por lo tanto, lo normal sería emplear estos datos como entrada del problema de diseño. Sin embargo, las expresiones desarrolladas están en función de la distribución simétrica Cp,sim y antisimétrica Cp,anti . Basta con aplicar la definición de ambos coeficientes para poder utilizar los datos de entrada habituales. Cp,sim =Cp,E+Cp,I 2(2.32) Cp,anti =Cp,E−Cp,sim ;Cp,anti =Cp,I−Cp,E 2(2.33) A modo representativo se presentan una serie de gráficas en las que se puede apreciar la diferencia entre las distribuciones de presiones en el extradós e intradós y las distribuciones simétricas y antisimétricas empleadas en las ecuaciones. 0 0.2 0.4 0.6 0.8 1 −0.5 0 0.5 1 Cp x Cp Cp simétrico Cp antisimétrico (a) NACA 0012 , α=0° 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 Cp x Cp Cp simétrico Cp antisimétrico (b) NACA 0012 , α=3° 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 Cp x Cp Cp simétrico Cp antisimétrico (c) NACA 2424 , α=0° 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 Cp x Cp Cp simétrico Cp antisimétrico (d) NACA 2424 , α=3° Figura 2.4 Distribuciones de presión Cp. 3 Resolución computacional Fo rmuladas las ecuaciones que proporcionarán los valores de espesor ze(x)(2.27) y curvatura zc(x)(2.31) , se procede a su implementación numérica usando el programa de cálculo Matlab. En este tercer capítulo se van a exponer todos los elementos usados para la resolución computacional de las ecuaciones, desde los diversos métodos numéricos hasta la obtención de los parámetros de entrada Cp,EyCp,I. 3.1 Problema de espesor A priori, el problema simétrico es el más complejo de resolver. El motivo de esta afirmación reside en el hecho de que se trata de una ecuación integral cuya incógnita ze no sólo se encuentra en el integrando, sino que además está derivada, aunque esto último merece menos importancia. Para intentar obtener la solución más exacta posible de este problema se han planteado tres métodos: método de Vortex-Lattice, series de Fourier y fórmula de inversión de Söhngen. Se hace hincapié en resolver este problema porque supone la parte central del programa que se quiere implementar para diseñar perfiles de forma inversa, puesto que el cálculo de la curvatura no presenta ninguna dificultad añadida. Un detalle a tener en cuenta es que si la distribución de presiones introducida fuera simétrica ( Cp,E=Cp,I ), el perfil buscado sería simétrico y estaría posicionado a ángulo de ataque nulo. Esto implica que bastaría con resolver el problema de espesor para conocer la geometría. No obstante esa situación será poco frecuente porque la utilidad de los perfiles se centra en su capacidad de generar sustentación a partir de la diferencia de presiones entre extradós e intradós. Sin embargo, será un planteamiento útil para comprobar el correcto funcionamiento de los programas implementados en este apartado. 3.1.1 Método de Vortex-Lattice En primer lugar se decidió adaptar el método de Vortex-Lattice por tres motivos: la formulación del problema que resuelve es semejante al aquí tratado, es fácil de programar, y por estar familiarizado con él tras la asignatura de Aerodinámica II y conocer su efectividad. Anderson [ 5 ] y Katz [ 15 ] explican que se trata de una metodología ampliamente empleada en aerodinámica desde los años 70. Es usado para la estimación de las características aerodinámicas como la sustentación y la resistencia inducida sobre cuerpos de geometría, espesor y orientación arbitraria. Para el caso de alas en 3D, está técnica numérica consiste en mallar el ala dividiéndola en pequeños paneles i j e imponer sobre cada uno de ellos una circulación Γi j tal que se cumpla la condición de impenetrabilidad (fig. 3.1.a) . Esta filosofía puede ser adaptada para el estudio de perfiles, que es el caso de este documento. En la literatura existen referencias a esta adaptación a 2D como Método de paneles. Al tratarse de una técnica para la estimación de la sustentación y la resistencia inducida, no tiene en cuenta el espesor ni la viscosidad porque estas fuerzas tienen su origen en el problema de curvatura. Se parte de la misma expresión deducida en el problema de curvatura (3.1) , con el objetivo de determinar la circulación Γ del perfil. En primer lugar, siguiendo la analogía 3D, se aproxima la curvatura del perfil por 35 36 3. Resolución computacional (a) Mallado en 3D (b) Adaptación a 2D (c) Panel elemental Figura 3.1 Vortex-Lattice. N elementos rectos de longitud h que simulan ser placas planas. Adoptando la simplificación de pequeñas perturbaciones, el mallado es trasladado al eje x (fig. 3.1.b). Sobre cada uno de estos paneles j se impone un torbellino de circulación Γj=γj·h , situado en la posición xg j . Adicionalmente se supone que la intensidad de circulación γj es uniforme en cada placa. Este mallado permite que la integral sea aproximada como un sumatorio, derivando en una expresión donde las incógnitas son las N densidades de circulación γj . Para terminar de cerrar el problema es necesario evaluar la igualdad en N puntos diferentes. Con este propósito se define sobre cada panel otro punto xci , conocidos como puntos de colocación. De esta manera se obtiene un sistema de ecuaciones N×N(3.2). La clave del éxito de este método está en la posición escogida para cada uno de los puntos xg j y xci . Tomando un elemento genérico j de longitud h , el primer punto xg j se sitúa a una distancia h/4 del comienzo del elemento j , mientras que el punto de colocación xci se posiciona a h/4 del final del elemento (fig. 3.1.c) . Esta elección no es ni mucho menos al azar. Al situar el torbellino a esa distancia y evaluarla en el otro punto, se calcula de forma exacta la circulación sobre una placa plana ( Γ=αcπU∞ ). En conclusión, las buenas aproximaciones conseguidas por este método están fundamentadas en la aproximación del perfil mediante placas planas, ya que es capaz de reproducir de forma exacta la circulación sobre ellas. U∞ dzc dx (x) = Zc 0 γo(xo) 2π(xo−x)dxo(3.1) U∞ dzc dx (xci) = N ∑ j=1 γjh 2π(xg j −xci)−→ bi=Ai j ·γj(3.2) 3.1 Problema de espesor 37 Para demostrar la eficacia de este método, se va a considerar una placa plana posicionada a un ángulo de ataque α cualquiera, cuya curvatura es zc=−αx . Al resolver el sistema numéricamente, se obtiene la distribución de intensidades γj a lo largo del eje x . Al compararlas con la distribución teórica ec. (3.3) , como se muestra en la fig. 3.2, se observa que los resultados numéricos casan perfectamente con los teóricos, es decir, esta metodología de integración es capaz de reproducir los resultados teóricos para las intensidades γ(x) . Se puede concluir entonces que Vortex-Lattice es un método de discretización efectivo para resolver problemas integrales. 0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 2.5 3 3.5 γ/U∞ x γ numérico γ teórico Figura 3.2 Intensidad de torbellinos γ(x) para placa plana. γ U∞ =2α1+cosθ sinθdonde θ=arccos1−2x c(3.3) Volviendo al problema de espesor que se busca resolver, si se hecha un vistazo a la ecuación planteada (3.4) y a la que soluciona el método de Vortex-Lattice (3.1) , puede verse que son problemas integrales del mismo tipo. Es por ello que se procede a su resolución empleando este proceso (3.5) . Sin embargo, hay que tener en cuenta que los valores conseguidos se corresponderán con la pendiente en los puntos xg j y no con el espesor del perfil. No obstante estos últimos son calculados al integrar de forma numérica la pendiente (3.6) . Por último, para que estos resultados se puedan combinar con los resueltos por el problema de curvatura, se decide interpolar los valores a los puntos de colocación xci. Cp,sim(x) = 2 πZc 0 dze dx (xo) (xo−x)dxo(3.4) Cp,sim(xci) = N ∑ j=1 2 π dze dx (xg j)h (xg j −xci)−→ bi=Ai j ·dze dx j (3.5) ze(x) = Zx 0 dze dx (xo)dxo(3.6) 3.1.2 Series de Fourier Las series de Fourier constituyen una herramienta matemática básica del análisis de Fourier empleado para estudiar funciones periódicas a través de la descomposición de dicha función en una suma infinita de funciones sinusoidales mucho más simples. Es una aplicación usada en muchas ramas de la ingeniería, que incluyen desde análisis vibratorio, acústica, óptica, hasta procesamiento de imágenes y señales. Dada una función periódica f(t)de periodo T, la serie de Fourier que aproxima dicha función es: f(t)≈a0 2+ ∞ ∑ n=1ancos2nπ Tt+bnsin2nπ Tt (3.7) donde a0,anybnson los coeficientes de Fourier. Mediante esta herramienta y con ayuda de otros resultados que ya se verán, es posible transformar la ecuación integral planteada para el problema de espesor, en otra expresión más fácil de resolver. Con esta manipulación lo que se pretende es que en lugar de tener que calcular los valores del espesor ze(x) , las incógnitas pasen a ser los coeficientes de Fourier que modelan la función ze. 44 3. Resolución computacional Llegados a este punto, para garantizar la correcta resolución del flujo potencial alrededor de una geometría arbitraria, se decidió aprovechar el programa proporcionado por Gordillo [ 11 ], que fue utilizado en el grado durante las prácticas de Aerodinámica I. El código es capaz de calcular los coeficientes de presiones del extradós e intradós sobre un perfil NACA 4 dígitos genérico para corriente potencial. El programa está basado en la implementación numérica del método de Green [Anexo B]. Otra posible opción que se barajó fue usar el programa XFLR5 [Anexo A] . Se trata de un software CFD libre que cuenta con la opción de calcular la distribución de presiones para flujo no viscoso. A priori la intuición puede decantarse por pensar que los resultados proporcionados por el XFLR5 serán más precisos que los conseguidos por el método de Green, dados los precedentes del XFLR5. No obstante al tratarse de un programa CFD, no deja de ser una "caja negra" en la que se introduce una geometría y sin saber cómo opera, devuelve los coeficientes de presión. Se hace necesario contrastar los resultados entre ambos programas para tomar una decisión sobre cual emplear. A primera vista los datos proporcionados por ambos son muy similares entre sí. La mayor discrepancia se encuentra en el borde de salida. En estos últimos puntos se aprecia que el coeficiente de presiones tiende a la unidad, lo que se corresponde con un punto de remanso. Se observa en la fig. 3.7 que esta tendencia es más acusada con el método de Green. Sin embargo hay que aclarar que este fenómeno no se da en la realidad. El borde de salida no es un punto de remanso porque el gradiente de presiones necesario para anular la velocidad en dicho punto es muy desfavorable, lo que conduciría a un desprendimiento de la capa límite. 0 0.2 0.4 0.6 0.8 1 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Cp x Cp pot. Met. Green Cp pot. XFLR5 (a) NACA 0012 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 Cp x Cp pot. Met. Green Cp pot. XFLR5 (b) NACA 1424 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 Cp x Cp pot. Met. Green Cp pot. XFLR5 (c) NACA 2424 0 0.2 0.4 0.6 0.8 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Cp x Cp pot. Met. Green Cp pot. XFLR5 (d) NACA 4412 Figura 3.7 Comparación de distribuciones de presiones con flujo potencial a α=0°. 3.4 Datos de entrada: Coeficiente de presiones Cp45 Para tratar de aliviar este gradiente adverso predicho por el método de Green y hacer que los resultados sean aun más parecidos, se opta por realizar una aproximación lineal en el último punto. Se pretende mantener la misma pendiente en los últimos dos paneles. Para ello el método elegido es el teorema de Thales o más conocido como triángulos semejantes (fig. 3.8). Figura 3.8 Triángulos semejantes. Tras aplicar esta aproximación al método de Green, conviene volver a comparar las distribuciones de presiones de ambos programas. Como se muestra en la fig. 3.9, la diferencia pasa a ser prácticamente nula. Además se ha hecho una ampliación del borde de salida para confirmar la validez de la aproximación lineal. Puesto que ambos métodos proporcionan los mismos datos, se opta por emplear el método de Green implementado por Gordillo [ 11 ]. El motivo principal de esta elección se debe a que el programa está implementado en Matlab y no es una "caja negra". Puede ser adaptado para cualquier geometría e incorporado a las rutinas de calculo para demostrar la efectividad de las funciones desarrolladas en este proyecto. Si se quisiera utilizar el XFLR5, sería necesario generar una base de datos que recogiese las distribuciones de presiones para los diferentes NACA con distintas condiciones, para posteriormente ser cargadas Matlab. 46 3. Resolución computacional 0 0.2 0.4 0.6 0.8 1 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Cp x Cp pot. Met. Green Cp pot. XFLR5 (a) NACA 0012 0.975 0.98 0.985 0.99 0.995 1 1.005 1.01 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Cp x Cp pot. Met. Green Cp pot. XFLR5 (b) BS NACA 0012 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 Cp x Cp pot. Met. Green Cp pot. XFLR5 (c) NACA 1424 0.97 0.98 0.99 1 1.01 1.02 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Cp x Cp pot. Met. Green Cp pot. XFLR5 (d) BS NACA 1424 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 Cp x Cp pot. Met. Green Cp pot. XFLR5 (e) NACA 2424 0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Cp x Cp pot. Met. Green Cp pot. XFLR5 (f) BS NACA 2424 0 0.2 0.4 0.6 0.8 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Cp x Cp pot. Met. Green Cp pot. XFLR5 (g) NACA 4412 0.975 0.98 0.985 0.99 0.995 1 1.005 1.01 1.015 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Cp x Cp pot. Met. Green Cp pot. XFLR5 (h) BS NACA 4412 Figura 3.9 Comparación de distribuciones de presiones tras la aproximación lineal en el borde de salida. 4 Análisis de los resultados Fi nalmente tras el estudio de las ecuaciones y la posibilidad de resolverlas numéricamente con tres métodos diferentes, es hora de comprobar cómo de eficaces son y cuales son sus límites. Se mostrarán diferentes ejemplos de perfiles aerodinámicos, desde los más sencillos hasta los más complejos, con o sin ángulo de ataque, que se han usado en este proyecto para evaluar y verificar la validez de los códigos implementados. 4.1 Primeros pasos: perfiles parabólicos y elípticos Para comenzar, se escogieron los perfiles más sencillos presentes en la literatura. Generalmente en teoría potencial las primeras secciones que se estudian son los perfiles parabólicos y elípticos, debido a que su geometría está bien definida y es posible la resolución del campo fluido analíticamente. Estas propiedades suponen una gran ventaja para nuestros intereses por dos motivos. En primer lugar, al ser secciones simétricas, no es necesario resolver el problema de curvatura, por lo tanto los resultados sólo dependerán de las tres metodologías programadas que resuelven el problema de espesor, permitiendo verificar su buen funcionamiento. En segundo lugar, al poder calcular de forma analítica el campo fluido alrededor del perfil, será posible compararlo con la distribución de presiones ofrecida por Met. Green y ver como afecta a la geometría resultante. 4.1.1 Perfiles parabólicos Estos perfiles también son conocidos como lenticulares. Su geometría se define en función del parámetro t , que indica el valor del semiespesor máximo que se alcanza en la parte media de la sección (fig. 4.1). Bajo estas características y teniendo en mente que son simétricos ( zp,E=−zp,I ), la geometría viene dada por la siguiente ecuación: Figura 4.1 Perfil lenticular. zp,E(x,t) = ze(x,t) = 4tx(1−x)(4.1) donde: dze dx =−8t(x−1 2)(4.2) Empleando la ec. (2.27) vista en el capítulo 2, se procede a conseguir la expresión analítica que determina la distribución de presiones Cpsobre ambas caras del perfil. Cp(x,t) = −2 π"l´ ım ε→0Zx−ε 0 −8txo−1 2 x−xo dxo+l´ ım ε→0Z1 x+ε−8txo−1 2 x−xo dxo#(4.3) 47 48 4. Análisis de los resultados Cp(x,t) = −16t πl´ ım ε→0"Zx−ε 0 1 2−xo+x−x x−xo dxo+Z1 x+ε 1 2−xo+x−x x−xo dxo#= =−16t πl´ ım ε→0Zx−ε 0 x−xo x−xo dxo+1 2−xZx−ε 0 dxo x−xo+Z1 x+ε x−xo x−xo dxo+1 2−xZ1 x+ε dxo x−xo= =−16t πl´ ım ε→0x−ε−1 2−xln(x−xo)|x−ε 0+1−x−ε−1 2−xln(x−xo)|1 x+ε= =−16t πl´ ım ε→0x−ε−1 2−xlnε x+1−x−ε−1 2−xlnx−1 −ε= =−16t πl´ ım ε→01−2ε−1 2−xlnε x+ln1−x ε, obteniendo, Cp(x,t) = 16t π1 2−xln1−x x−1(4.4) Al compararla gráficamente con el que proporciona el Met. Green se observa cierta semejanza entre ellas, pero no llegan a ser iguales (fig. 4.2). En general, tomando como referencia los valores de la presión en los extremos y en el punto medio, Green proporciona unos valores más bajos, lo que se traduce en una mayor velocidad de la corriente sobre el perfil, provocando posiblemente que el espesor de la sección buscada sea mayor. Además, se hace notable que la diferencia entre ellas es mayor si aumenta el espesor t del perfil, ya que modelan el flujo de distinta forma. No obstante, lo importante es comprobar que los códigos funcionan bien y son capaces de reproducir la geometría deseada a partir de las distribuciones de presiones dadas. 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 Cp x Cp Analitico Cp Met. Green (a) t=0.1 0 0.2 0.4 0.6 0.8 1 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 Cp x Cp Analitico Cp Met. Green (b) t=0.2 Figura 4.2 Distribución de presiones asociada a perfiles parabólicos Cp(x,t). 4.1 Primeros pasos: perfiles parabólicos y elípticos 49 0 0.2 0.4 0.6 0.8 1 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (a) CpAnalítico 0 0.2 0.4 0.6 0.8 1 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (b) CpMet. Green Figura 4.3 Perfiles parabólicos con t=0.1. 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (a) CpAnalítico 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (b) CpMet. Green Figura 4.4 Perfiles parabólicos con t=0.2. A través de la fig. 4.3 y fig. 4.4 se muestra que los códigos reproducen casi con exactitud el perfil para el caso de proporcionar como entrada la distribución de presiones exacta ( Cp Analítica), con independencia del espesor t . Sin embargo, en el caso de introducir los valores determinados por Green, los perfiles conseguidos son menos aproximados al objetivo, aunque tienen cierto parecido en la forma. Efectivamente, como se adelantó para este segundo caso, las secciones calculadas tienen un mayor espesor debido a que la los gradientes de velocidades impuestos son mayores que los exactos. Entrando dentro de los resultados proporcionados por cada método, con independencia tanto del espesor t como de los parámetros de entrada Cp , destacan las aproximaciones logradas por las series de Fourier. Respecto a Vortex-Lattice y Söhngen, sorprende que a pesar de ser metodologías tan diferentes, los resultados logrados por ambos sean prácticamente idénticos en cualquiera de los casos aquí expuestos. Si echamos un vistazo a las geometrías, se observa que tanto para Vortex-Lattice como para Söhngen, el borde de salida es abierto, mientras que Fourier consigue geometrías cerradas. Esto se debe a que no se impuso ninguna condición de cierre en estos dos programas. Dependerá de que la presión tenga sentido físico en el borde de salida y tienda a recuperarse sin llegar a ser un punto de remanso. Todos estos hechos confirman aun más si cabe la gran importancia que tiene la distribución de presiones como parámetro de entrada y su influencia en los resultados. Será necesario realizar más pruebas con los programas y verificar si las tendencias aquí observadas persisten con otras secciones. 50 4. Análisis de los resultados 4.1.2 Perfiles elípticos El estudio de las geometrías elípticas supone un pequeño paso adelante en busca de los límites que presentan los códigos. Son secciones simétricas cuya geometría es conocida en función del espesor t (ec. (4.5) ). No obstante, a diferencia de los perfiles parabólicos, la pendiente de los elípticos (ec. (4.6) ) cuenta con dos singularidades que se corresponden con el borde de ataque ( x=0 ) y el de salida ( x=1 ). Efectivamente, en dichos puntos la pendiente se vuelve infinita, lo que puede traer problemas tanto a la hora de caracterizar el flujo en esas zonas como en reproducir la geometría. Figura 4.5 Perfil elíptico. (x−1)2 (1 2)2+z2 p (t 2)2=1 ; zp,E(x,t) = ze(x,t) = tpx(1−x) (4.5) donde: dze dx =t1−2x 2px(1−x)(4.6) Para tratar de modelar el campo de presiones sobre el perfil se emplean tres formas diferentes. La primera se corresponde con el resultado analítico obtenido a través de la ec. (2.27) . Se caracteriza por ser constante (ec. (4.7) ), pero para emular la caída de presión y la recuperación, se han supuesto como puntos de remanso tanto el B.A como el B.S. Su demostración puede ser consultada en Katz [ 15 ]. En segundo lugar se ha aprovechado una expresión (ec. (4.8) ) dada en los trabajos de Mecánica de Fluidos II, que caracteriza la velocidad del fluido sobre un cuerpo elíptico con el origen situado en el centro de la elipse y cuyos extremos se sitúan en −1 y 1 . Para su uso es preciso determinar los valores de la presión dentro del rango [−1,1] y posteriormente trasladarlos a [0,1] con el cambio de variable x=X+1 2 . Por último, el flujo se aproxima numéricamente con el Met. Green. Cp(t) = −2t(4.7) u0 e(X,t) = 1+t q1+t2X2 1−X2−1(4.8) De esta manera se intenta caracterizar el flujo en los extremos de la mejor forma posible, y evitar así que los resultados proporcionados por los programas no sean los adecuados debido a una mala aproximación de la distribución de presiones. En el caso de no coincidir con la potencial, será de utilidad para ver las desviaciones que se producen respecto a la geometría buscada. También es conveniente analizar el efecto que tiene el espesor t , pues como se ha justificado a lo largo del documento se está interesado en cuerpos esbeltos. Representando la presión con los tres modelos explicados para los espesores t=0.2 y t=0.6 se obtienen la fig. 4.6.a y fig. 4.6.b, respectivamente. Al compararlas, se ve como los modelos de presión implementados coinciden para espesores esbeltos, mientras que para cuerpos romos aparecen desviaciones respecto al potencial Cp=−2t . No obstante, aunque las presiones sean similares para el caso de t=0.2 , esto no garantiza que los perfiles calculados de forma inversa reproduzcan la forma deseada, y menos sin saber cómo son afectados por las singularidades del B.A y B.S. Analizando más detenidamente las desviaciones que aparecen con t=0.6 , la forma de las curvas recuerda a las vistas para perfiles parabólicos, donde la caída de presión y su recuperación es más gradual. Al igual que en el apartado anterior, estos modelos también presentan un valor mínimo de presión menor que la analítica. Así que, es de esperar que la forma que proporcionen los códigos para este caso sea parecida a un perfil lenticular cuyo espesor máximo sea mayor que el espesor de la elipse buscada. Lo único que se puede asegurar tras estudiar estas gráficas es que el Met. Green con el que se introducen los datos, sólo es válido para perfiles esbeltos. Al emplearse fuera de ese rango las distribuciones de presiones que proporciona no se corresponden con el flujo potencial que utilizan los programas. Por lo tanto, cuando se estudien perfiles alares reales deberá tenerse presente esta limitación y ceñirse a ese rango. 4.1 Primeros pasos: perfiles parabólicos y elípticos 51 Tras emplear los programas con los tres modelos de presión para las dos geometrías elípticas propuestas, se muestran los resultados en la fig. 4.7 y fig. 4.8. 0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 1 1.5 2 2.5 3 x Cp Cp=−2t ue Analitico Cp Met. Green (a) t=0.2 0 0.2 0.4 0.6 0.8 1 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 3 x Cp Cp=−2t ue Analitico Cp Met. Green (b) t=0.6 Figura 4.6 Distribución de presiones asociada a perfiles elípticos Cp(x,t). 0 0.2 0.4 0.6 0.8 1 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (a) Cp=−2tAnalítico 0 0.2 0.4 0.6 0.8 1 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (b) Cp(ue)Analítico 0 0.2 0.4 0.6 0.8 1 −0.2 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 0.2 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (c) CpMet. Green Figura 4.7 Perfiles elípticos con t=0.2. 52 4. Análisis de los resultados 0 0.2 0.4 0.6 0.8 1 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (a) Cp=−2tAnalítico 0 0.2 0.4 0.6 0.8 1 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (b) Cp(ue)Analítico 0 0.2 0.4 0.6 0.8 1 −0.5 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 0.5 x z Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (c) CpMet. Green Figura 4.8 Perfiles elípticos con t=0.6. En esta ocasión a pesar de contar con la distribución de presiones analítica Cp=−2t , los resultados logrados no son tan buenos como los obtenidos en perfiles lenticulares. La razón está en la pendiente de la curva de presiones. Al ser gradual, los métodos basados en la pendiente del espesor dze dx evolucionan suavemente, haciendo que el espesor crezca o disminuya poco a poco, como en el caso de las secciones parabólicas. Sin embargo en perfiles elípticos, lo que ocurre es que la presión cambia bruscamente en el B.A, provocando una pendiente de espesor que tiende a infinito y que descontrola el crecimiento del espesor. Al llegar al B.S. la presión vuelve a cambiar, haciendo que los valores de la pendiente de espesor sean negativos, pero no llegan a ser lo suficientemente negativos como para invertir el crecimiento descontrolado del espesor. El fenómeno anterior se aprecia muy bien cuando se comparan los tres modelos de distribución de presiones con los perfiles conseguidos. Con la analítica, al ser tan brusco el cambio de presión, el espesor no disminuye. Para los otros dos ejemplos, Cp(ue) y Met. Green, si se aprecia la disminución del espesor hacia el B.S al contar con un cambio de presión menos brusco. En general se mantienen las tendencias observadas para secciones lenticulares. Las series de Fourier consiguen muy buenas aproximaciones a la geometría buscada, a pesar de que los datos de entrada no se correspondan con los potenciales. El motivo es que de los tres métodos es el único que no se basa en el cálculo de la pendiente de espesor, si no en los coeficientes de Fourier para modelar la función de espesor. Por su parte, Vortex-Lattice y Söhngen siguen proporcionando resultados muy similares entre sí. Además ambos se ven afectados por las singularidades en las pendientes del B.A. y B.S., descontrolando los valores de espesor calculados en estos métodos. 4.2 Perfiles NACA 4 dígitos 53 En definitiva con estos primeros perfiles se ha podido verificar el correcto funcionamiento de los métodos implementados para la resolución del problema de espesor, y su comportamiento frente a los datos de entrada. Aunque puedan parecer que los resultados logrados no son los esperados, han ayudado a comprender que la inversión del problema no es tan sencilla como aparenta ser. Aun así, viendo los perfiles de presiones típicos de las secciones alares esbeltas como los mostrados en la fig. 3.9, se espera que no existan problemas con las pendientes al tratarse de curvas más suaves. 4.2 Perfiles NACA 4 dígitos Para continuar con la experimentación de los programas, se decide pasar al siguiente nivel, reproducir perfiles reales que se utilicen en la práctica. Con este propósito se escogen los perfiles NACA de 4 dígitos, secciones cuya geometría es sistemática, conocida y fácilmente reconocible con la nomenclatura de los perfiles. Se podría haber escogido cualquier otra familia, pero como ya se comentó en el primer capítulo, aun hoy en día los NACA son los más empleados para el diseño de alas. Con el fin de probar la consistencia de los programas se han seleccionado una serie de perfiles que abarcan todos los rangos posibles, desde los más esbeltos hasta los más grueso, con o sin curvatura. La misión de los códigos será intentar reproducirlos. Posteriormente esos mismos perfiles se han vuelto a evaluar pero con un ángulo de ataque α=2° . Como único parámetro de entrada se cuenta con las distribuciones de presiones proporcionadas por el Met. Green definidas con 500 nodos. Respecto a los métodos implementados: Vortex-Lattice y Söhngen se fijan con 200 paneles, y las series de Fourier con 25 coeficientes an para el problema de espesor y 200 paneles para el de curvatura. Los resultados obtenidos se muestran en las páginas siguientes desde la fig. 4.9 hasta la fig. 4.20. Cada uno de los perfiles se presenta junto a su distribución de presiones. A primera vista los resultados logrados son muy satisfactorios. Las aproximaciones que se consiguen respecto a las geometrías buscadas, son muy buenas, llegando a ser casi exactas en algunos casos. Se observan pequeñas desviaciones cuando el perfil comienza a tener una curvatura y/o espesor considerables, viendose incrementadas con el ángulo de ataque. Efectivamente, los primeros perfiles al ser tan esbeltos quedan bajo las simplificaciones realizadas considerando pequeñas perturbaciones, mientras que los últimos quedan algo más alejados, de ahí las desviaciones. Observando estas desviaciones, según las imágenes, parecen que quedan localizadas en el borde de ataque y sobre la parte media del extradós. La explicación de este fenómeno está en el pico de succión del borde de ataque y los valores que toma la pendiente de presión. Al igual que ocurría para los perfiles elípticos, si el cambio de presión es muy brusco, provoca que la pendiente de espesor alcance valores superiores a los reales. Esto se transmite desde el borde de ataque hasta la parte media del extradós, recuperando la normalidad poco a poco al llegar al borde de salida. Además teniendo en cuenta que el pico de succión aumenta con el ángulo de ataque, es lógico que las desviaciones se vean incrementadas con α=2°. Al estudiar los datos proporcionados por cada uno de los métodos y compararlos, destaca la serie de Fourier como método de inversión sobre los otros dos. No sólo es más preciso, si no que además es más rápido y cuenta con la condición de cierre, lo que le proporciona ventaja frente a los otros. Sin embargo, tampoco hay que despreciar los perfiles conseguidos por Vortex-Lattice y Söhngen. En especial sorprende Vortex-Lattice, ya que un método tan sencillo de implementar como es la conversión de una integral numérica en un sumatorio pueda dar resultados tan próximos. En comparación, Söhngen que logra casi los mismos perfiles que el método anterior, es más complicado de programar y mucho más lento con los cálculos. Por los motivos expuestos, parece que existe un claro vencedor como método inverso. Llegados a este punto es conveniente realizar un análisis más exhaustivo donde se controlen cada uno de los parámetros que influyen sobre los tres métodos. De esta manera se podrá determinar su sensibilidad ante el cambio de las condiciones, necesario para seguir evaluando las capacidades de los códigos. 60 4. Análisis de los resultados 0 0.2 0.4 0.6 0.8 1 −0.5 0 0.5 1 Cp x N = 25 N = 50 N = 100 N = 200 N = 300 (a) Datos Cpde entrada con Met. Green: N nodos 0 0.2 0.4 0.6 0.8 1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 z x Perfil N = 25 N = 50 N = 100 N = 200 N = 300 (b) Met. Vortex-Lattice 0 0.2 0.4 0.6 0.8 1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 z x Perfil N = 25 N = 50 N = 100 N = 200 N = 300 (c) S. Fourier 0 0.2 0.4 0.6 0.8 1 −0.15 −0.1 −0.05 0 0.05 0.1 0.15 z x Perfil N = 25 N = 50 N = 100 N = 200 N = 300 (d) I. Söhngen Figura 4.22 Estudio de los parámetros de entrada para NACA 1412 , α=0°. 4.3.3 Ángulo de ataque Llega el turno de analizar la posición de los perfiles respecto a la corriente libre, definido como ángulo de ataque. En particular, durante la formulación fue considerado junto a la curvatura para el problema antisimétrico sin dividir el problema entre ellos. Por lo tanto, aunque se trate de un perfil simétrico, si este está posicionado a un ángulo de ataque no nulo, el bloque que resuelve el problema de curvatura debe actuar para corregir los valores de espesor y colocarlos al mismo ángulo de ataque. Lo importante de este asunto es que se ha visto en resultados previos, que el ángulo de ataque incrementa la diferencia de presiones entre ambas caras, generando más sustentación. Fenómeno que ocurre debido al aumento del pico de succión en el extradós. Llegará un momento en que la aceleración del flujo sea tal, que las simplificaciones hechas basándose en pequeñas perturbaciones no sean válidas. Así que, el objetivo de este apartado es determinar el valor máximo del ángulo de ataque ( αmax ) para el cual los resultados calculados con los códigos siguen dando buenas aproximaciones. Para este estudio es preciso partir de un perfil simétrico que los programas reproduzcan perfectamente a ángulo de ataque nulo. Es necesario que sea simétrico para eliminar la influencia de la curvatura en el problema antisimétrico, dejando solo el ángulo de ataque. El proceso a seguir es simple, incrementar el ángulo de ataque de la sección deseada y compararla con los perfiles calculados por los programas. Bajo estas restricciones se decide escoger un perfil NACA 0012 y someterlo a ángulos comprendidos entre 0 ° y 7 ° . Respecto al resto de parámetros sobre los que el usuario tiene control, se mantiene lo establecido en los anteriores apartados: Green con 500 nodos y los códigos con 200 paneles y 25 coeficientes. A fin de tener más datos para el análisis, los perfiles resultantes se acompañan de una gráfica que muestra como varía la distribución de presiones con el aumento del ángulo de ataque. 4.3 Análisis de sensibilidad 61 0 0.2 0.4 0.6 0.8 1 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 Cp x α = 0º α = 1º α = 3º α = 5º α = 7º (a) Distribucion de presiones NACA 0012 para varios α 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (b) Perfiles α=0° 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (c) Perfiles α=1° 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (d) Perfiles α=3° 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (e) Perfiles α=5° 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (f) Perfiles α=7° Figura 4.23 Estudio del ángulo de ataque para NACA 0012. Analizando las secciones de las imágenes, se confirma lo que se ya se sospechaba. A medida que aumenta el ángulo de ataque, los códigos no son capaces de reproducir la forma deseada. En particular el error se encuentra en el cálculo del espesor y no en la posición del perfil, ya que si se observan los puntos inicial y final, estos coinciden con la sección buscada, y por tanto, también el ángulo que forman. Una vez más nos encontramos con que la singularidad en el borde de ataque provoca errores en el cálculo del espesor, pero esta vez afecta a los tres métodos por igual. La razón está en que al calcular Cp,sim necesario para el problema simétrico, en le B.A no se corresponde con la distribución para α=0° . Tras estudiar los datos, se puede establecer el límite de ángulo de ataque entre 2°y 3°. 62 4. Análisis de los resultados Realmente el usuario no tiene control sobre este parámetro, pero según las conclusiones sacadas, lo que el usuario debe ver es que si la posición del perfil resultante tiene un alto ángulo de ataque, probablemente los resultado no sean del todo válidos. 4.3.4 Curvatura y espesor Por último, pero no menos importante, hay que determinar donde se encuentran los límites de curvatura y espesor. Según lo visto en apartados anteriores, es lógico pensar que también existan unos valores de curvatura y de espesor a partir de los cuales las perturbaciones no sean pequeñas. De igual forma que para el ángulo de ataque, se quiere ver como actúa el bloque de resolución del problema antisimétrico cuando sólo existe curvatura. No obstante, también hay que tener presente el problema simétrico, ya que parece que es el que limita la capacidad de los códigos en el diseño inverso. Así que además de los límites, se estudiará qué influye más en las desviaciones de los resultados, si la curvatura o el espesor. Para distinguir los efectos producidos por la curvatura y el espesor sobre los resultados, se decide escoger como punto de partida dos perfiles simétricos con distinto espesor a α=0° . En este caso concreto se han seleccionando el NACA 0012 y NACA 0024 (fig. 4.24.a y fig. 4.24.b ). Estas secciones representan los dos extremos que nos preocupan, perfiles esbeltos y gruesos. Tomándolos como referencia, se procede a incrementar el valor de la curvatura poco a poco, dando lugar a los resultados mostrados en las fig. 4.24.c - fig. 4.24.h . La primera columna muestra los perfiles esbeltos (NACA XX12) y la segunda los gruesos (NACA XX24). Los parámetros que controla el usuario y con los que han sido obtenidos los resultados, siguen fijados en: 500 nodos, 200 paneles y 25 coeficientes. Observando los perfiles resultantes de los tres códigos, queda patenten que aunque no sean iguales a los buscados, las aproximaciones conseguidas en todos los casos son muy buenas. Son tanto mejores cuanto más simples son los perfiles, debido a lo que ya se ha explicado en varias ocasiones, el pico de succión en el extradós. Respecto a los métodos empleados, las series de Fourier continua siendo el mejor de los tres. Sin embargo, todo esto sólo reafirma las observaciones realizadas en los apartados anteriores, sin aportar nada nuevo. Lo interesante es fijarse en las desviaciones. Analizando los perfiles simétricos de referencia, se ve como para el más esbelto (fig. 4.24.a) los método reproducen casi con exactitud la sección, mientras que para el más grueso (fig. 4.24.b) aparecen las primeras desviaciones. Son muy discretas, pero existen. A medida que aumenta la curvatura, se hacen cada vez más notables. En el caso del perfil esbelto, las desviaciones comienzan a notarse cuando la curvatura alcanza un valor considerable (fig. 4.24.g). Con estos datos y después de experimentar con otros valores, es conveniente establecer un límite máximo de espesor del 18%, valor también sugerido por Carmona [ 13 ]. Por otro lado, con la curvatura se podría ser más permisivo dependiendo del espesor de la sección, pero sería recomendable que no excediera el 4%. Ante estos hechos se puede afirmar que la variable que más influye en los resultados proporcionados por los códigos, es el espesor de las secciones que se intenta imitar. Dejando en un segundo plano a la curvatura. Como ya ocurrió con la posición del perfil, ni la curvatura ni el espesor son parámetros que el diseñador controle, son fruto de la pendiente de la distribución de presiones y la diferencia entre el extradós e intradós. Deberá actuarse de la misma forma. Si durante el proceso de diseño los perfiles obtenidos cuentan con un espesor excesivo y algo de curvatura, el usuario debe ser consciente de que los resultados tienen ciertas desviaciones. Para terminar, es interesante fijarse y comparar las imágenes del apartado anterior (fig. 4.23) con las aquí presentes. Partiendo de la misma sección simétrica de control (NACA 0012), se ve como las desviaciones debidas al ángulo de ataque son más importantes que las provocadas por la curvatura y el espesor. Era de esperar, puesto que el pico de succión se ven influenciado en mayor medida por la posición del perfil que por la geometría. Por lo tanto, a la hora de evaluar la precisión de los resultados será necesario valorar primero la posición, y después su espesor. 4.3 Análisis de sensibilidad 63 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (a) NACA 0012 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (b) NACA 0024 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (c) NACA 1412 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (d) NACA 1424 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (e) NACA 2412 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (f) NACA 2424 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (g) NACA 4412 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Met. Vortex−Lattice S. Fourier I. Söhngen Perfil (h) NACA 4424 Figura 4.24 Estudio de la curvatura y el espesor con perfiles NACA 4 dígitos, α=0°. 64 4. Análisis de los resultados 4.4 Flujo real vs Flujo potencial Todo el proyecto ha sido enfocado al estudio de la corriente potencial alrededor de perfiles, flujo que se define de forma teórica bajo ciertas simplificaciones y que permite explicar muchos de los fenómenos que ocurren en la Aerodinámica. No obstante, su carácter teórico es el que puede hacer dudar de su utilidad en casos reales. Como defiende Gordillo [ 11 ], la Aerodinámica potencial es aquella que se dedica a estudiar el flujo alrededor de objetos fuselados en los casos en los que los esfuerzos viscosos se encuentren confinados en regiones muy estrechas del dominio fluido. Podría parecer que esta situación es muy particular y queda lejos de la realidad, sin embargo no es así, el flujo alrededor tanto de las superficies sustentadoras como de las estabilizadoras y de control, puede ser calculado usando técnicas de la Aerodinámica potencial en todo el dominio excepto en la capa límite y la estela. Atendiendo a este razonamiento, se entiende que los valores de presión empleados para el diseño inverso de perfiles, calculados del flujo potencial, son válidos siempre y cuando la capa límite sea una región muy estrecha y no interfiera en exceso con el flujo exterior. Dicho de otro modo, puede que la diferencia entre las distribuciones reales y potenciales de presión sea mínima. Así que, llegados a este punto resulta interesante plantearse la siguiente situación ¿qué ocurre si a los códigos se les introducen datos reales en lugar de los potenciales? A priori, según la sensibilidad vista que presenta los programas ante ligeros cambios en los datos de entrada, se espera que los perfiles resultantes no sean iguales, aunque esto dependerá de lo diferentes que sean las dos distribuciones entre sí. En especial, Vortex-Lattice y Söhngen, son los más sensibles ante está circunstancia, mientras que Fourier parece tolerar mejor las discrepancias. Para contrastar estas ideas basta con tomar ambas distribuciones y comparar los resultados. En esta ocasión no se usará el Met. Green implementado en Matlab, si no el programa XFLR5, ya que es capaz de calcular las presiones en el caso real con viscosidad y en el potencial. Los datos son obtenidos mediante un mallado uniforme de 300 paneles y almacenados en un archivo .txt para su importación en Matlab, conforme a la explicación dada en el Anexo A. Una vez introducidas las presiones, se aplica la misma rutina de cálculo bajo los parámetros fijados. El perfil escogido es un NACA 4412 con α=0° . Al tratarse de un perfil esbelto con una curvatura moderada y sin ángulo de ataque, se confía en que los resultados asociados a la distribución potencial sean buenos, haciendo más fácil la comparación respecto a los reales. Realizados los cálculos, los resultados son mostrados en la fig. 4.25, donde el color azul se asocia con el flujo real viscoso y el rojo con la corriente potencial. Al comparar las presiones, se ve como efectivamente ambas son muy parecidas excepto en las inmediaciones del borde de salida. Según los datos mostrados en XFLR5, la diferencia está causada por el crecimiento de la capa límite, la cual pasa de ser laminar a turbulenta en ese entorno. El flujo se frena y la presión se recupera. Como se intuyó, debido a estás discrepancias entre los efectos reales y el comportamiento teórico, las secciones asociadas a la corriente viscosa no consiguen reproducir la forma deseada cuando se emplea Vortex-Lattice o Söhngen (fig. 4.25.b y fig. 4.25.d). Sin embargo, la cosa cambia cuando se usan las series de Fourier, obteniendo unos perfiles muy próximos para ambos flujos. Puede verse en la fig. 4.25.c que en el entorno de la transición de la capa límite, la sección "real" tiene un poco más de espesor y termina con una mayor pendiente, lo que concuerda con las presiones. Aunque ya se tenia nociones sobre la tolerancia de este método ante desviaciones en los datos de entrada, sorprende ver lo flexible que es. Tras los datos obtenidos, se concluye entonces que por muy pequeños que puedan parecer los efectos viscosos al estar confinados en una región muy estrecha del espacio, siempre afectarán al flujo exterior en mayor o menor medida. Además habrá que tener especial cuidado con el crecimiento de la capa límite. Si la capa límite alcanza un determinado espesor, se da la transición a turbulento, lo que nos aleja del rango de las pequeñas perturbaciones, e invalida las ecuaciones empleadas. Todo esto deberá ser tenido en cuenta por el diseñador a la hora de definir la distribución de presiones según las características deseadas. Por otra parte se ha visto que las series de Fourier aceptan como datos de entrada las presiones reales. Para tratar de demostrar que no es un caso aislado, se quiere seguir experimentando con otros ejemplos. 4.4 Flujo real vs Flujo potencial 65 0 0.2 0.4 0.6 0.8 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Cp x Cp experimental Cp potencial (a) Distribución de presiones XFLR5 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Cp exp Cp pot (b) Met. Vortex-Lattice 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Cp exp Cp pot (c) S. Fourier 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Cp exp Cp pot (d) I. Söhngen Figura 4.25 Comparación flujo real y potencial: NACA 4412 , α=0°. Siguiendo el mismo proceso anterior que admite los datos exportados de XFLR5, se opta por probar perfiles de otras familias. Se busca que sean diferentes, pero al mismo tiempo deben respetar las restricciones fijadas durante el análisis de sensibilidad. En particular se pone especial atención a que sean esbeltos. Adicionalmente se van a suponer posicionados a ángulo de ataque nulo. Tras consultar la base de datos UIUC [ 19 ], se han escogido los siguientes: Clark Y , Eppler 221 y Göttingen 394. Se observa que para el perfil Clark Y (fig. 4.26), las presiones son muy parecidas salvo por el efecto de la transición en el borde de salida, igual que ocurría antes. En cambio, en las otras dos secciones restantes (fig. 4.27 y fig. 4.28), existe una mayor diferencia entre las distribuciones, destacando el caso del Eppler 211. Sin embargo a pesar de estos efectos, los resultados siguen la misma tendencia que lo ilustrado para la sección NACA. Las secciones obtenidas asociadas a los flujos reales son muy buenas, con pequeñas desviaciones puntuales localizadas en el borde de salida. Ante estos resultados, se confirma que las series de Fourier son el mejor método para resolver el problema de diseño inverso. Importante destacar que esta capacidad de aceptar como datos de entrada las presiones reales supone una enorme ventaja, porque no se tratan de valores teóricos estimados analíticamente de las ecuaciones simplificadas, si no de datos reales propios de un túnel de viento. Aunque bien es cierto que estas distribuciones se corresponden con flujos subsónicos bajos, donde la capa límite es laminar sobre la mayor parte del perfil. No obstante, para este proyecto esto supone un primer paso hacia un método viable para el diseño inverso de perfiles. 66 4. Análisis de los resultados 0 0.2 0.4 0.6 0.8 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Cp x Cp experimental Cp potencial (a) Distribución de presiones 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Cp exp Cp pot (b) Perfiles Figura 4.26 S. Fourier: Clark Y , α=0°. 0 0.2 0.4 0.6 0.8 1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 1.2 Cp x Cp experimental Cp potencial (a) Distribución de presiones 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Cp exp Cp pot (b) Perfiles Figura 4.27 S. Fourier: Eppler 211 , α=0°. 0 0.2 0.4 0.6 0.8 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Cp x Cp experimental Cp potencial (a) Distribución de presiones 0 0.2 0.4 0.6 0.8 1 −0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4 z x Cp exp Cp pot (b) Perfiles Figura 4.28 S. Fourier: Göttingen 394 , α=0°. 4.4 Flujo real vs Flujo potencial 67 Esta última idea sugiere que, desde el punto de vista de aplicación práctica, sería interesante comparar los resultados aquí obtenidos con los ofrecidos por otros métodos reales para el diseño inverso. De los métodos presentados en el capítulo 1, se decide escoger el Ball-Spine Algorithm (BSA) porque está desarrollado unicamente para solucionar el problema de diseño, y no como los otros dos que son rutinas para resolver de forma general cualquier problema no lineal. La comparación entre ambos métodos se hace tomando el perfil BOEING 727 estudiado en [ 4 ] con ángulo de ataque nulo. Debe puntualizarse que aunque las condiciones no son las mismas para ambos casos, sirve para ilustrar la precisión de ambos programas. Viendo la fig. 4.29, es evidente que no existe ni punto de comparación entre ellos. A base de realizar muchas iteraciones resolviendo el problema directo, BSA es capaz de corregir la diferencia con la distribución de presiones objetivo, mientras que Fourier invierte el problema y lo resuelve en un único paso, consiguiendo peores aproximaciones. La ventaja de BSA sobre Fourier, está en que separa el problema de diseño en dos partes independientes: un bloque que resuelve el campo fluido alrededor de un objeto, y otro bloque que corrige la forma de la sección según las diferencias entre las presiones. De esta forma consigue que el cálculo de las presiones sea más preciso, pudiendo considerar efectos viscosos en diferentes regímenes como subsónico o transónico. Fourier se ve limitado por las simplificaciones hechas que permitieron invertir el problema de diseño. (a) Ball-Spine Algorithm 0 0.2 0.4 0.6 0.8 1 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 z x Cp exp Cp pot (b) S. Fourier Figura 4.29 Comparación de métodos de diseño inverso , Boeing 727, α=0°. Conclusiones El desarrollo de superficies sustentadoras pasa inevitablemente por el estudio de los perfiles aerodinámicos. Como se ha demostrado, existe un vínculo entre el campo fluido y la forma de las secciones. Es vital su comprensión porque las fuerzas aerodinámicas tienen su origen en la deflexión de las líneas de corriente provocadas por el perfil. No todos los cuerpos son aptos para generar sustentación, sólo aquellos que sean esbeltos. También hay que poner especial atención a los esfuerzos viscosos confinados en la capa límite. A pesar de ser una región muy estrecha adyacente a la superficie, tiene gran influencia en el comportamiento aerodinámico de los perfiles, y en particular en la entrada en pérdida. Dependiendo de la aplicación se buscan unas características aerodinámicas concretas, de ahí que surjan las diferentes formas, clasificadas en familias. Fruto de la experimentación surgieron los primeros catálogos, que sirvieron como base para la selección de perfiles durante los primeros años. El desarrollo de métodos de diseño de perfiles aerodinámicos está justificado por tres motivos: el diseño de secciones que quedan fuera del alcance de los catálogos, diseño eficiente y específico para cada aplicación, y el ahorro económico en concepto de pruebas y tiempo. Muchas personas han trabajado en este campo, pero todos los procesos son agrupados en dos filosofías: métodos directos y métodos inversos. Los directos toman como punto de partida la forma del perfil aerodinámico y calculan la distribución de presiones sobre ellos, para posteriormente optimizar la sección. En cambio los inversos, comienzan definiendo la distribución de presiones objetivo y determinan la forma del perfil. La ventaja de los métodos inversos sobre los directos es que las propiedades aerodinámicas como sustentación, resistencia y capa límite pueden ser transformadas en una distribución de presiones. En esta materia destacó Eppler, al implementar un código que permitía diseñar perfiles aerodinámicos con las características de capa límite especificadas y analizar secciones de geometría conocida. Actualmente, se investigan técnicas iterativas como BSA, o incluso otras pertenecientes a la rama de inteligencia artificial como ANN y EA. En este proyecto se ha tratado de dar los primeros pasos en busca de un método de diseño de perfiles. Para ello se ha partido de las ecuaciones generales de Navier-Stokes de Mecánica de fluidos y se han aplicado las simplificaciones clásica para flujo potencial considerando pequeñas perturbaciones, llegando al problema de Laplace que resuelven el flujo exterior alrededor de un objeto en la Aerodinámica potencial. Dividiendo el problema en espesor y curvatura se encuentran las ecuaciones que relacionan directamente la geometría con las presiones. Invirtiendo estas expresiones, se consigue un método inverso que permite conocer la forma del perfil a partir de las distribuciones de presiones impuestas en una única iteración. En un primer momento puede pensarse que el flujo potencial no tiene validez en aplicaciones reales, pero a lo largo del proyecto se ha justificado que esto no es así. Cierto es que no tiene en cuenta los efectos viscosos aunque estos si influyen sobre las presiones, pero como en todos los avances, se comienza estudiando los casos sencillos. Dada la complejidad de las expresiones invertidas para el cálculo de la geometría, se opta por probar tres métodos diferentes para resolver dichas ecuaciones: un método basado en paneles de Vortex-Lattice, resolución mediante series de Fourier, y usando la fórmula de inversión de Söhngen. Así que el proyecto pasa por evaluar y determinar cual es el mejor de los tres. Para comprobar la integridad de estos programas se decide probar si son capaces de reproducir perfiles reales a partir de sus distribuciones de presiones. Las secciones consideradas han sido perfiles NACA de 4 dígitos de todo tipo a diferentes ángulos de ataque. 69 76 A. XFLR5 Figura A.4 Pantalla inicial del análisis de perfiles. El análisis a realizar será un Batch analysis, que se encuentra en la pestaña Analysis. Sobre este tipo de procesos, Fernández[ 10 ] explica que se llama Batch cuando las instrucciones insertadas en él, son realizadas de forma secuencial y automática por el ordenador. A continuación lo que se debe hacer es configurar el análisis que se va ha realizar sobre el perfil. La ventana de configuración (fig. A.5) se encuentra dividida en varias partes, lo que hace más cómodo la preparación. Figura A.5 Configuración del tipo de análisis. En primer lugar se elige el perfil a analizar. Puede ser aquel que fue seleccionado con el desplegable, o elegir desde un listado varios de los cargados previamente. El siguiente paso consisten en definir el tipo de análisis que se quiere. XFLR5 cuenta con 4 tipos, diferenciándose entre sí por las variables empleadas en el barrido principal. El más convencional y usado en este proyecto es el Tipo 1 porque permite un control directo sobre el número de Reynolds. El siguiente bloque permite forzar puntos de transición de la capa límite tanto en extradós como en intradós. El último de los bloques determina la variable sobre la que realizar un barrido secundario una vez fijada la variable del barrido principal. 77 Por ejemplo, para el caso aquí ilustrado se ha elegido el perfil NACA 1424 para analizar. Posteriormente se selecciona un análisis para un rango de Re que va desde 100.000 hasta 150.000 , con un incremento de 10.000 , fijando M=0,044 y considerando que la corriente no está perturbada N=9 . No se fuerzan los puntos de transición. Por último, para cada valor de Re se analizará el perfil a diferentes ángulos de ataque α que van desde −5°hasta 5°con un paso de 1°. Se pulsa el botón Analyze y comienza el análisis, dando como resultado las diferentes gráficas. Para estudiar los resultados XFLR5 dispone de dos interfaces diferentes. Figura A.6 Polar view: características globales. Por defecto la primera que aparece es la Polar view, donde se representan todos los datos obtenidos para los diferentes perfiles (fig. A.6). Da una visión general de la evolución de las características de los perfiles según las condiciones definidas en el análisis. La comparación entre las diferentes geometrías es inmediata. Las gráficas que ofrece de forma predeterminada son: la polar del perfil ( Cl/Cd ), evolución de la sustentación ( Cl/α ), comportamiento del punto de transición en extradós ( α/Xtr t ), coeficiente de momento ( Cm/α ) y evolución de la eficiencia ( Cl Cd/α ). Sin embargo, si se quiere se pueden cambiar las variables representadas en la gráfica haciendo click derecho en Current Graph/ Define Graph Setting. Por otra parte la segunda interfaz, OpPoint view, permite ver la distribución de presiones Cp sobre el perfil para cada una de las situaciones en las que fue analizado (fig. A.7). Sin duda alguna, esta es la gráfica que resulta más interesante para los propósitos del proyecto. En la barra superior, mediante los tres desplegables se eligen el perfil, las condiciones del barrido principal y el valor de la variable en el barrido secundario. Siguiendo con el ejemplo, se escoge el perfil NACA 1424 a Re =100.000 , M=0.044 y N=9 , y en el tercer desplegable se selecciona el ángulo de ataque α=3°. Situado en la esquina inferior izquierda se ven datos relacionados con la geometría del perfil, mientras que en el lado derecho se ven las características aerodinámicas del perfil para la configuración seleccionada. Mediante el menú de la derecha, en la parte de display, se pueden marcar las opciones para que muestre sobre el perfil la capa límite (Show BL), la presión y el centro de presiones (Show Pressure), o animar la gráfica para ver su evolución respecto a la variable secundaria. No obstante lo que se necesita para el proyecto son los valores de la presión asociados a un flujo potencial (Show Inviscid Curve). Esta opción se encuentra en el menú Operating Points/Cp Graph. Al marcarla, aparece sobre la gráfica una curva puntea representando los valores correspondientes a flujo no viscoso. 78 A. XFLR5 Figura A.7 Distribucion de presiones Cppara NACA 1424 a α=3°. Por último, sólo queda exportar estos valores. Se vuelve al menú Operating Points/Cp Graph y se elige Export Graph. Esto genera un archivo .txt que posteriormente puede ser importado a Matlab. Si se abre (fig. A.8), se observa que los datos están agrupados en columnas, y sin más que borrar la primera línea, pueden ser leídos como una matriz en Matlab. Figura A.8 Presentación de los datos exportados. Anexo B Método de Green Es un método numérico usado ampliamente en la resolución de problemas ingenieriles basado en la fórmula de Green y la discretización del entorno. Según esta fórmula, dadas dos funciones armónicas Φ y Ψ acotadas en un dominio conexo Ωcuya frontera Σes continua a trozos, se cumple ec. (B.1). [6] ZΣ(Ψ∇Φ −Φ∇Ψ)·−→ n dσ=ZΩΨ∇2φ−Φ∇2ψdω=0(B.1) Aplicado al problema aerodinámico tratado en este proyecto, resuelve el problema (ec. (B.2) ) expresando el valor del potencial sobre cualquier punto de la superficie del objeto como la superposición de una distribución continua de dos tipos de soluciones básicas: fuentes situadas sobre la superficie, ψj=ln −→ x−−→ xj  , donde −→ xj∈Σs y dobletes situados sobre la superficie del objeto y orientados según la normal exterior al mismo, −→ n·∇ψj . Para llegar a esa expresión (ec. (B.3) ), se parte de la integral de Green y se evalúa en la superficie Σc=Σ∞SΣεSΣsSΣe+SΣe− , mostrada en la fig. B.1. La demostración completa del cálculo de la integral de superficie Σc para el potencial de velocidades φ , puede ser consultada en Introducción a la aerodinámica potencial [11]. Figura B.1 Superficies para la integral de Green [11].                    ∇2φ=0 |−→ x|→∞,φ→cosαx+sinαz −→ x∈Σs,−→ n·∇φ=0 + Condición de Kutta-Joukowski                    (B.2) La ecuación de Green para el cálculo del potencial se expresa como: φ(xj,zj) = 2(xjcosα+zjsinα)+ 1 πZΣs φ(x,z)(−→ x−−→ xj)·−→ n −→ x−−→ xj 2dσ+Γ πZΣe (−→ x−−→ xj)·(−−→ kcosα+−→ isinα) −→ x−−→ xj 2dσ (B.3) donde −→ xj=xj−→ i+zj−→ k es el punto sobre la superficie del perfil Σs sobre el que se sitúa el doblete y por tanto, es un punto fijo durante la integración, mientras que −→ x=x−→ i+z−→ k es un vector de posición sobre Σs y Σe que sí varía en la integración. Como puede comprobarse, el potencial está expresado como la suma de los potenciales de una distribución continua de dobletes orientados según las normal de las superficies Σs y Σe de intensidades φ(x,z) y Γ , respectivamente. Para el caso de flujo potencial, la intensidad de las fuentes es nula, debido a que se cumple la condición de impenetrabilidad.[11] 79 80 B. Método de Green Figura B.2 Discretización de la geometría en N paneles. Definición de punto medio y vectores tangente y normal de un panel k[11]. Las N ecuaciones de la expresión (B.3) , una por cada doblete, forman el sistema de ecuaciones que se utiliza para resolver numéricamente el flujo potencial alrededor de un sólido de geometría dada. En primer lugar se discretiza la superficie del objeto Σs en N paneles. Cada panel k , tiene una longitud ∆k definida desde el nodo k hasta el nodo k+1 (fig. B.2). Se definen N+1 nodos, donde el nodo N+1 se hace coincidir con el primer nodo, consiguiendo así los N paneles. Conociendo el vector posición de estos nodos, −→ rk=xk−→ i+zk−→ k , se puede calcular el punto medio y la longitud de cada panel kcomo: −→ mk=mkx−→ i+mkz−→ k=xk+1+xk 2−→ i+zk+1+zk 2−→ k ∆k=q(xk+1−xk)2+(zk+1−zk)2 , además de los vectores tangente y normal: −→ tk=xk+1−xk ∆k −→ i+zk+1−zk ∆k −→ k −→ nk=zk+1−zk ∆k −→ i−xk+1−xk ∆k −→ k , e incluso la distancia entre puntos medios: ∆0 k=q(m(k+1)x−mkx)2+(m(k+1)z−mkz)2 Colocando los N dobletes en los puntos medios de cada panel k , y suponiendo que el valor del potencial φkes constante a lo largo de cada panel, la ecuación integral pasa a expresarse como: φj=2(mjx cosα+mjz sin α)+ 1 π N ∑ k=1 φkZ∆k 0 (−→ x−−→ mj)·−→ nk −→ x−−→ mj 2ds+Γ πZ∞ 0 (−→ x−−→ mj)·(−−→ kcosα+−→ isinα) −→ x−−→ mj 2ds (B.4) donde el valor de las integrales es: si k=j,Z∆k 0 (−→ x−−→ mj)·−→ nk −→ x−−→ mj 2ds =0,(−→ x−−→ mk)⊥−→ nk si k6=j,Z∆k 0 (−→ x−−→ mj)·−→ nk −→ x−−→ mj 2ds =arctan"∆k+(−→ rk−−→ mj)·−→ tk (−→ rk−−→ mj)·−→ nk#−arctan"(−→ rk−−→ mj)·−→ tk (−→ rk−−→ mj)·−→ nk# Z∞ 0 (−→ x−−→ mj)·(−−→ kcosα+−→ isinα) −→ x−−→ mj 2ds ="π 2sign[(1−mjx)sinα+mjz cosα]−arctan (1−mjx)cosα−mjz sinα (1−mjx)sinα+mjz cosα!# 81 Cabe destacar que se cuenta con N ecuaciones, pero el sistema tiene N+1 incógnitas, los N potenciales φj y el valor de la circulación Γ . Por lo tanto es necesario buscar una ecuación adicional para que el sistema quede cerrado. Dicha condición adicional es la condición de Kutta, que impone que la velocidad en el borde de salida tiene que ser igual cuando se calcule por el intradós que por el extradós. Considerando que el borde de salida se encuentra en el nodo M+1, se debe cumplir que: φM−φM−1 ∆0 M−1 =φM+1−φM+2 ∆0 M+1 (B.5) En resumen, el método de Green permite hallar el potencial sobre la superficie del perfil mediante la resolución del siguiente sistema lineal de N+1 ecuaciones ( j∈[1,N] ) para las N+1 incógnitas [φ1,φ2,...,φN,Γ] , que son los valores del potencial φken los puntos medios de cada panel y la circulación Γ. [11]          ∑N k=1(δjk −Rjk)φk−Γfj=2(mjx cosα+mjz sinα) φM−φM−1 ∆0 M−1−φM+1−φM+2 ∆0 M+1=0 (B.6) donde: Rjk =     1 πarctanhbs jk+∆k bn jk i−arctan hbs jk bn jk i k6=j 0k=j      con bs jk = (−→ rk−−→ mj)·−→ tkybs jk = (−→ rk−−→ mj)·−→ nk y fj=1 π"π 2sign[(1−mjx)sinα+mjz cosα]−arctan (1−mjx)cosα−mjz sinα (1−mjx)sinα+mjz cosα!# Una vez resuelto el sistema, los coeficientes de presiones sobre el perfil se obtienen como: Cp(−→ rk) = 1−φk−φk−1 ∆0 k−12 (B.7) El código asociado a la implementación de la resolución de este sistema se encuentra disponible en el Apéndice E de Introducción a la aerodinámica potencial [11]. Anexo C Código Met. Vortex-Lattice % ========================================================================= % RESOLUCION FLUJO POTENCIAL ---> VORTEX LATTICE % ========================================================================= % Función que resuelve el problema simétrico y antisimétrico con flujo % potencial, de forma inversa. Partiendo de los coeficientes de presiones % de extradós e intradós, el objetivo es conocer cual es la forma del perfil. % % Tanto para resolver el problema de ESPESOR como el de CURVATURA, se % emplea Vortex-Lattice para calcular las derivadas de espesor y curvatura, % lo que permite determinar la geometria del perfil. % % En el caso del problema de ESPESOR se calcula la intensidad de las % fuentes, ya que están relacionadas con el cambio de la geometría. Como % dato conocido, se cuenta con la velocidad de perturbación horizontal, % gracias a su relación con la distribucion de presiones. % % Por otro lado, para el problema de CURVATURA, la distribución de % presiones nos indica cual es la intensidad de los torbellinos, teniendo % que calcular las velocidades de perturbación verticales, las cuales están % relacionadas con la derivada de la curvatura. % %========================================================================== % % Parametros de entrada %-------------------------------------------------------------------------- % Re Reynolds % s Vector adimentsional de coordenadas horizontal donde están % medidos Cp [0 - 1] % Cp_extrados Distribución de presiones en el extrados % Cp_intrados Distribución de presiones en el intrados % % Parámetros de salida %-------------------------------------------------------------------------- % zpE Geometría de extradós % zpI Geometría de intradós % X Vector adimensional de coordenadas horizontales donde esta % calculada la geometría del perfil % %========================================================================== function [zpE,zpI,X]=fp_VL(Re,s,Cp_extrados,Cp_intrados) 83 84 C. Código Met. Vortex-Lattice %% Parametros inciales del problema %Parametros del perfil c=1; % Cuerda del perfil [m] N=200; % Número de paneles h=c/N; % Long de cada panel [m] disp([' Vortex-Lattice --> N º Paneles = ' num2str(N) ]) %Puntos de centro de torbellinos/fuentes Xgj de cada panel xg(1:N)=((1:N)-1)*h+h/4; %Puntos de colocacion Xci xc(1:N)=xg+h/2; %Condiciones de la corriente incidente [Sistema internacional] rhoinf=1.225; % Densidad [kg/m^3] [Aire] nu=1.5e-5; % Viscosidad/Densidad [m^2/s] [Aire] Uinf=Re*nu/c; % Velocidad [m/s] %% Problema Simétrico o de Espesor CpI_xc=interp1(s,Cp_intrados,xc); % Interpolación de Cp para intradós y CpE_xc=interp1(s,Cp_extrados,xc); % extradós en los puntos de colocación Xci u_sim=-Uinf/4*(CpI_xc+CpE_xc); % Velocidad de perturbación simétrica % Matriz Aij for i=1:N for j=1:N A(i,j)=h/pi*Uinf/(xc(i)-xg(j)); end end %Término independiente bi b(1:N,1)=u_sim'; %Resolucion del sistema Dze(1:N,1)=A\b; % Derivada del espesor en los puntos Xgj FDze = @(x)interp1(xg,Dze,x); % Función para interpolar la Dze for j=1:N ze_g(j)=quad(FDze,0.0025,xg(j)); % Espesor en los puntos Xgj end ze=interp1(xg,ze_g,xc); % Espesor en los puntos Xci ze(N)=ze_g(N); %% Problema Antisimétrico o de Curvatura CpI_xg=interp1(s,Cp_intrados,xg); % Interpolación de Cp para intradós y CpE_xg=interp1(s,Cp_extrados,xg); % extradós en los puntos de colocación Xgj gamma_ext=-Uinf/2*(CpE_xg-CpI_xg); % Densidades de circulación en los puntos Xgj % Matriz Bij for i=1:N for j=1:N B(i,j)=h/2/pi/(xg(j)-xc(i)); end 85 end Dzc=B*gamma_ext'/Uinf; % Derivada de la curvatura en los puntos Xci FDzc = @(x)interp1(xc,Dzc,x); % Función para interpolar la Dzc for i=1:N zc(i)=quad(FDzc,0.0075,xc(i)); % Curvatura en los puntos Xci end %% Problema completo. Determinación de la geometría X=xc; zpE=zc+ze; % Extradós zpI=zc-ze; % Intradós Bibliografía [1] The naca airfoil series,http://people.clarkson.edu/~pmarzocc/AE429/The%20NACA%20airfoil% 20series.pdf . [2] Aircraft Aerodynamics and Design Group Stanford University, Airfoil design methods,http://adg. stanford.edu/aa241/airfoils/airfoildesign.html. [3] , History of airfoil development,http://adg.stanford.edu/aa241/airfoils/airfoilhistory.html. [4] M. Nili Ahmadabadi, F. Ghadak, and Mohammadi M., Subsonic and transonic airfoil inverse design via ball-spine algorithm, Elsevier, Computers and Fluids (2013). [5] J.D Anderson, Fundamentals of aerodynamics, 4ªed., McGraw-Hill, 2007. [6] Antonio Barrero Ripoll and Miguel Pérez-Saborid Sánchez Pastor, Fundamentos y aplicaciones de la mecánica de fluidos, McGraw-Hill. [7] L. Di Angelo and P. Di Stefano, An evolutionary geometric primitive for automatic design synthesis of functional shapes: The case of airfoils, Elsevier, Advances in Engineering Software (2014). [8] Milton Van Dyke, An album of fluid motion, 11ªed., The parabolic press, 2007. [9] R. Eppler and Dan M. Somers, A computer program for the design and analysis of low-speed airfoils, NASA TM 80210, 1980. [10] Lucas Fernández-Peña Mollá, Xflr5 v4.17 software de análisis de perfiles, planos y aviones: estudio y diseño de un cn-235, Universidad Carlos III de Madrid, 2012. [11] José Manue Gordillo Arias de Saavedra and Guillaume Riboux Acher, Introducción a la aerodinámica potencial, 1ªed., Paraninfo, 2012. [12] Ira H. Abbott and Albert E. Von Doenhoff, Summary of airfoil, NACA Report Nº824, 1945. [13] Anibal Isidoro Carmona, Aerodinámica y actuaciones del avión, 12ªed., Thomsom-Paraninfo, 2004. [14] Fernando Izaurieta and Carlos Saavedra, Redes neuronales artificiales, Charlas de fisica (1999). [15] Joseph Katz and Allen Plotkin, Low-speed aerodynamics, 2ªed., Cambridge University Press. [16] Jeffrey Keith Jepson, Enhancement to the inverse design of low speed natural laminar flow airfoils, Master’s thesis, North Carolina State University. Aerospace Engineering. MSc Thesis, 2003. [17] Alfonso Mateos Andaluz, Algoritmos evolutivos y algoritmos genéticos, Universidad Carlos III de Madrid. [18] Miguel Ángel Muñoz, Manual de vuelo,http://www.manualvuelo.com. [19] UIUC Applied Aerodynamics Group. Department of Aerospace Engineering, Uiuc airfoil coordinates database,http://m-selig.ae.illinois.edu/ads/coord_database.html. 93 94 Bibliografía [20] Gustavo H. Scarpin, Aerodinamica de perfiles, Instituto Universitario Aeronáutico. [21] Dan M Somers, Eppler airfoil design and anaysis code,http://www.airfoils.com/eppler.htm. [22] , Subsonic airfoil design,http://www.airfoils.com/airfoil.htm. [23] Theodore Theodorsen, Theory of wing section of arbitrary shape, NACA Report Nº411, 1931.