Predicción de arrastre para el Modelo Común ... - Uniandes
Transcript of Predicción de arrastre para el Modelo Común ... - Uniandes
Predicción de arrastre para el Modelo Común de Investigación desarrollado por la NASA (CRM), por medio de CFD.
Juan Sebastian Velandia Rodriguez
Proyecto de grado presentado para optar al título de Ingeniero Mecánico
Asesor del Proyecto Omar López,
MSc. Ingeniería Mecánica Universidad de los Andes PhD. Ingeniería Mecánica Universidad de Texas
Universidad de los Andes Facultad de Ingeniería
Departamento de Ingeniería Mecánica Bogotá D.C., Colombia
Mayo de 2013
Tabla de Contenido
Introducción ....................................................................................................................... 4 Fuerzas aerodinámicas ............................................................................................................. 4 Sustentación ............................................................................................................................................... 5 Arrastre ........................................................................................................................................................ 5 Ondas de choque ...................................................................................................................................... 7 Eficiencia aerodinámica ........................................................................................................................ 7 Control longitudinal ................................................................................................................................ 7
Dinámica de fluidos computacional. .................................................................................... 8 Aplicaciones ............................................................................................................................................... 8 Métodos numéricos y modelamiento. ............................................................................................. 8 Predicción de fuerzas .......................................................................................................................... 10
Objetivos ............................................................................................................................ 10 Metodología ...................................................................................................................... 10 Revisión geometría ................................................................................................................. 11 Selección malla ......................................................................................................................... 11 Modelo de turbulencia ........................................................................................................... 12 Valores de predicción ......................................................................................................................... 13 Residuales ................................................................................................................................................ 15 Selección ................................................................................................................................................... 16
Procedimiento .......................................................................................................................... 17 Resultados y análisis ...................................................................................................... 17 Convergencia de malla ........................................................................................................... 17 Validación resultados .......................................................................................................................... 17 Análisis de convergencia ................................................................................................................... 18 Estabilidad en coeficientes ............................................................................................................... 19
Predicción de arrastre ........................................................................................................... 19 Fuerzas aerodinámicas ...................................................................................................................... 20 Desempeño aerodinámico ................................................................................................................ 21 Condiciones de flujo ............................................................................................................................ 22
Ondas de choque ...................................................................................................................... 27 Control aerodinámico ............................................................................................................ 29 Costo computacional .............................................................................................................. 31
Conclusiones ..................................................................................................................... 32
Referencias ........................................................................................................................ 34 Anexos ................................................................................................................................. 35 Anexo A: Archivos de MATLAB. ........................................................................................... 35 Anexo B: Ejemplo journal utilizado. .................................................................................. 39 Anexo C: Manual de configuración ..................................................................................... 40
Lista de figuras Figura 1. Diagrama de cuerpo libre en un avión. (Moran, 1984) ................................... 4 Figura 2. Fuentes de arrastre en un avión comercial. (Bertin, 2002) ............................. 6 Figura 3. Onda de choque en vuelo transónico. (Moran, 1984) ........................................ 7 Figura 4. Curva típica Cm contra ángulo de ataque. (Yechout, 2003) ............................ 8 Figura 5. Exactitud en predicción viscosa para distintos tipos de malla. (Baker,
2005) ............................................................................................................................................. 12 Figura 6. Convergencia del coeficiente de sustentación con diferentes modelos de
turbulencia. ............................................................................................................................... 13 Figura 7. Convergencia del coeficiente de arrastre con diferentes modelos de
turbulencia. ............................................................................................................................... 14 Figura 8. Convergencia del coeficiente de momento con diferentes modelos de
turbulencia. ............................................................................................................................... 14 Figura 9. Residuales de simulación con modelo k-‐e .......................................................... 15 Figura 10. Residuales de simulación con modelo Spalart-‐Allmaras. ......................... 16 Figura 11. Residuales de simulación con modelo SST-‐kw. ............................................. 16 Figura 12. Dependencia de coeficientes con respecto a tamaño de malla. .............. 18 Figura 13. Convergencia de Cd con respecto a iteraciones. ........................................... 19 Figura 14. Coeficiente de sustentación contra ángulo de ataque. ............................... 20 Figura 15. Coeficiente de arrastre contra ángulo de ataque. ......................................... 20 Figura 16. Curvas polares para diferentes configuraciones. ......................................... 21 Figura 17. Eficiencia aerodinámica para diferentes configuraciones ........................ 22 Figura 18. Coeficiente de fricción para diferentes posiciones del estabilizados .. 25 Figura 19. Coeficiente de presiones para diferentes posiciones del estabilizador.
........................................................................................................................................................ 26 Figura 20. Coeficientes de presión en diferentes planos. ............................................... 28 Figura 21. Número de mach y líneas de corriente. ............................................................ 29 Figura 22. Coeficiente de momento contra ángulo de ataque. ...................................... 30
Lista de tablas Tabla 1. Detalles juego de mallas utilizado. .......................................................................... 12 Tabla 2. Comparación de resultados obtenidos y reportados. ..................................... 18 Tabla 3. Valores de Cd por elementos. .................................................................................... 23 Tabla 4. Angulo de trim para diferentes configuraciones. .............................................. 30 Tabla 5. Costo computacional de primer caso de estudio. ............................................. 31 Tabla 6. Costo computacional segundo caso de estudio.2 ............................................... 31
Introducción La mecánica de fluidos es una de las áreas de estudio de la mecánica que más fenómenos físicos busca modelar y entender. Distintos enfoques han sido utilizados para la compresión de los fluidos, hasta hace algunos años la experimentación y el análisis teórico eran los más importantes. Por medio de estas dos herramientas se buscaba modelar matemáticamente los comportamiento que se veían en el diario vivir. Sin embargo, hay eventos que son difíciles de observar o de medir, por lo que la experimentación parece no ser suficiente para generar una base teórica después de cierto nivel de detalle. Hace algunas décadas se comenzó a trabajar en un nuevo enfoque para la mecánica de fluidos, conocido como dinámica de fluidos computacional. El enfoque consiste en utilizar los modelos matemáticos conocidos y reproducirlos a través del tiempo, buscando que los resultados obtenidos sean comparables con la experimentación realizada. En particular, este proyecto de grado busca recrear simulaciones computacionales realizadas en el marco del 4th Drag Prediction Workshop (organizado por la AIAA) con el fin de generar un manual de usuario para este tipo de simulaciones. El documento hará una breve descripción de alguna teoría necesaria, recordará los objetivos planteados en un principio, mostrará la metodología utilizada, los resultados obtenidos y, finalmente, algunas conclusiones y sugerencias.
Fuerzas aerodinámicas En condiciones de vuelo estacionarias, cualquier vehículo aéreo se ve afectado por las mismas fuerzas: Empuje (T), sustentación (L), arrastre (D) y peso (W). Algunas de estas fuerzas se encuentran a favor del movimiento del vehículo en la dirección deseada (sustentación y empuje) mientras que otras se oponen al mismo (arrastre y peso). El diagrama de cuerpo libre de la Figura 1 muestra la dirección de cada una de las fuerzas mencionadas, en condiciones de vuelo normal.
Figura 1. Diagrama de cuerpo libre en un avión. (Moran, 1984)
De las fuerzas mencionadas, el interés se centra en las aerodinámicas, es decir, sustentación y arrastre. El peso es una fuerza dependiente de la masa del avión, se busca que sea mínimo para garantizar la condición de vuelo con niveles de sustentación más bajo. El empuje es proporcionado por los motores (turbinas) del vehículo, se busca que sea mínimo para que el consumo de combustible no sea mayor al necesario.
Sustentación La sustentación es la fuerza que se genera por fenómenos de diferencias de presiones presentes en el ala. Debido a que los perfiles alares tienden a ser delgados y alargados, la sustentación se debe mayoritariamente a los efectos que afectan la geometría de forma normal a la superficie, es decir, el diferencial de presiones. Los esfuerzos cortantes serán considerados más a fondo en términos de arrastre, no tanto en sustentación pues el área donde influirían en dirección “vertical” es muy pequeña. (Bertin, 2002) Al integrar el campo de presiones generado sobre el ala del avión se encuentra una fuerza resultante, en dirección ascendente, esta fuerza es la que se conoce como sustentación (L). Con el fin de comparar la sustentación en distintos tipos de vehículos, se define un coeficiente de sustentación adimensional, se la siguiente manera (Moran, 1984):
𝐶! =𝐿
12𝜌!𝑉!
!𝑆
Para definir este coeficiente sólo se necesita conocer la densidad del flujo en la condición de flujo libre, la velocidad del flujo libre y un área de referencia (S) del cuerpo sobre el que se calcula. Este coeficiente será utilizado más adelante. En un avión comercial, los valores de sustentación dependen de sus condiciones de vuelo, sin embargo, valores entre 0.1 y 0.7.
Arrastre La fuerza de arrastre corresponde al arrastre producido por distintos componentes en el avión, debido a distintos fenómenos del fluido a su alrededor. En particular, los fenómenos que más aportan al arrastre de un vehículo de este tipo son la fricción en la superficie y el arrastre inducido por sustentación, como se puede ver en la Figura 2. La fricción superficial es un fenómeno que depende directamente de la viscosidad del fluido por el que pasa el vehículo. Se genera un esfuerzo cortante en toda la superficie del avión. De manera similar a la sustentación, el arrastre por esta fuente se obtiene integrando el valor del esfuerzo cortante en cada punto de la superficie del avión. Los efectos de fricción son responsables de cerca del 50% del arrastre del avión. Con el fin de comparar este efecto en puntos discretos, se utiliza un parámetro adimensional conocido como el coeficiente de fricción:
Figura 2. Fuentes de arrastre en un avión comercial. (Bertin, 2002)
𝐶! =𝜏!
12𝜌𝑉!
!
Se le conoce como arrastre inducido al generado por los elementos del avión que están dispuestos para generar sustentación. De cierto modo, se puede interpretar que al integrar las presiones sobre la superficie del avión, la fuerza resultante tiene una componente que se opone al movimiento del avión, conocida como arrastre inducido (Moran, 1984). Desde un punto de vista menos detallado, este tipo de arrastre se puede interpretar como el generado por la interrupción del flujo por parte del cuerpo que se encuentra en él. Esta interrupción genera una diferencia de presiones (más alta flujo arriba, más baja flujo abajo) que implica una fuerza que se opone al movimiento, conocida como arrastre inducido (Hoerner, 1965). Los efectos de este tipo de arrastre llegan casi al 30% del arrastre total del avión. Las comparaciones a niveles de presión se hacen por medio de un parámetro adimensional, el coeficiente de presión, definido como:
𝐶! =𝑝 − 𝑝!12𝜌𝑉!
!
De manera similar a la sustentación, se define un coeficiente de arrastre para los vehículos aéreos con el fin de comparar su comportamiento en términos de arrastre. (Hoerner, 1965)
𝐶! =𝐷
12𝜌!𝑉!
!𝑆
En un avión comercial, los valores de arrastre dependen de sus condiciones de vuelo, sin embargo, valores inferiores a 0.05.
Ondas de choque En condiciones de vuelo transónico (velocidades cercanas a la velocidad del sonido) un fenómeno a tener en cuenta que la generación de ondas de choque. La geometría del vehículo atravesando el flujo a las velocidades mencionadas genera efectos de compresibilidad en el flujo donde las propiedades del aire se vuelven discontinuas. Esta compresibilidad genera zonas donde el flujo pasa de ser subsónico a transónico, conocida como onda de choque, mostrada en la Figura 3.
Figura 3. Onda de choque en vuelo transónico. (Moran, 1984)
La formación de estas ondas altera la presión en la superficie del elemento aumentándola en el sentido opuesto al flujo, generando un aporte adicional a los términos de arrastre, conocido como wave drag.
Eficiencia aerodinámica La eficiencia aerodinámica es una figura de mérito utilizada comúnmente en análisis aerodinámicos. Se trata simplemente de la razón entre el coeficiente de sustentación y el de arrastre. Entre más alta sea la sustentación de un vehículo aéreo, mayor será su eficiencia aerodinámica. Del mismo modo, entre menos arrastre genere, su eficiencia será mayor. (Torenbeek, 2009)
𝐸 =𝐶!𝐶!
El comportamiento de esta figura de mérito con respecto al ángulo de ataque suele ser una parábola donde el punto máximo representa la condición de vuelo ideal, en términos aerodinámicos. La importancia de este punto es muy evidente en vuelos largos, donde el consumo de combustible es vital en términos energéticos y ambientales, la condición mencionada disminuye considerablemente el consumo del vehículo.
Control longitudinal Si bien la evaluación de control longitudinal en una aeronave es un tema basto que no está directamente relacionado con las simulaciones realizadas, se hará una pequeña aproximación a la relación entre la posición del estabilizador y la estabilidad estática de la aeronave. Dado que el estabilizador es un elemento relacionado directamente al control de la aeronave, su posición altera considerablemente el modo de vuelo de la misma.
Este efecto se logra por medio de la alteración del coeficiente de momento (Cm), por lo que es necesario generar una gráfica de este valor con respecto al ángulo de ataque de la aeronave, para encontrar la condición de vuelo en la que este elemento cumple su función de anular los momentos sobre el vehículo.
Figura 4. Curva típica Cm contra ángulo de ataque. (Yechout, 2003)
El corte de la curva con el eje horizontal muestra el valor de interés, denominado ángulo de trim. Al modificar la inclinación del estabilizador (iH) se puede ajustar este ángulo al valor deseado, una inclinación positiva genera que el ángulo de trim disminuya mientras, una positiva lo aumenta. Esto se debe a que inclinaciones positivas disminuyen el momento generado por el elemento, mientras las negativas, lo aumentan. (Yechout, 2003)
Dinámica de fluidos computacional. La dinámica de fluidos computacional (CFD) es tal vez una de las ramas de la mecánica que evoluciona con mayor velocidad hoy en día. El estado del arte de esta herramienta es difícil de describir pues su evolución es constante. Cerrando un poco el campo a las aplicaciones aerodinámicas, algunos autores hacen un breve recuento de lo que se ha hecho y el estado actual de esta aplicación (Vos, 2002).
Aplicaciones El papel de la CFD en el diseños aerodinámicos es cada vez más importante. Inicialmente, se le consideró simplemente una herramienta que supliera la necesidad de experimentos en túneles de viento a escala, donde no se tenían los túneles. Con el tiempo, se ha convertido en una herramienta necesaria para todos los campos del diseño de aviones, no sólo la aerodinámica. Las aplicaciones del CFD al diseño de aeronaves se puede dividir en 6 campos: aspecto del vehículo, desempeño, integración de sistemas, control y estabilidad, cargas para diseño estructural y diseño aeroelástico. Cada una de estas áreas se apoya en esta herramienta para conocer y optimizar los parámetros que hacen de lo vehículos a diseñar seguros, eficientes y confiables.
Métodos numéricos y modelamiento. Lo anterior se logra gracias a, en parte, una gran evolución en los métodos numéricos, la capacidad de cómputo y las mejores en el modelamiento físico de los fenómenos aerodinámicos. En particular, se puede notar que en la última década los solucionadores basados en las ecuaciones de Navier-‐Stokes son los
más utilizados (tratando estas ecuaciones como diferencias finitas), debido a que su predicción del flujo es bastante cercana a la realidad (Vos, 2002). Su formulación diferencial se muestra a continuación, para todas las direcciones (White).
𝜌𝑔! −𝜕𝑝𝜕𝑥 + 𝜇
𝜕!𝑢𝜕𝑥! +
𝜕!𝑢𝜕𝑦! +
𝜕!𝑢𝜕𝑧! = 𝜌
𝑑𝑢𝑑𝑡
𝜌𝑔! −𝜕𝑝𝜕𝑦 + 𝜇
𝜕!𝑣𝜕𝑥! +
𝜕!𝑣𝜕𝑦! +
𝜕!𝑣𝜕𝑧! = 𝜌
𝑑𝑣𝑑𝑡
𝜌𝑔! −𝜕𝑝𝜕𝑧 + 𝜇
𝜕!𝑤𝜕𝑥! +
𝜕!𝑤𝜕𝑦! +
𝜕!𝑤𝜕𝑧! = 𝜌
𝑑𝑤𝑑𝑡
Además, se deben resolver las ecuaciones de conversación de masa, momento y energía, mostradas en ese orden. Para su explicación más detallada se recomienda consultar la guía teórica para el usuario del software FLUENT (ANSYS, 2009), donde se encuentra el fondo teórico de cada uno de los términos.
𝜕𝜌𝜕𝑡 + ∇. 𝜌𝑣 = 𝑆!
𝜕𝜕𝑡 𝜌𝑣 + ∇. 𝜌𝑣𝑣 = −∇.p+ ∇. 𝜏 + 𝑝𝑔 + 𝐹
𝜕𝜕𝑡 𝜌𝐸 + ∇. 𝑣 𝜌𝐸 + 𝑝 = ∇. 𝑘!""∇.T− ℎ!𝐽!
!
+ 𝜏!"" . 𝑣 + 𝑆!
Por otro lado, el modelamiento de la turbulencia y del desprendimiento de capa parecen fenómenos que aún están por fuera de las capacidades de predicción de este solucionador. Se han generado distintos modelo de turbulencia y aún se estudia cuál de todos es más acorde con al realidad. Algunos autores (Spalart, 2010) describen lo que ha sido el modelamiento por medio de RANS (Reynolds-‐Averaged-‐Navier-‐Stokes) en los últimos años, y porque sigue siendo tan importante hoy en día, después de la aparición de nuevos métodos de modelamiento de turbulencia. Se puede encontrar una comparación de los 3 modelos de turbulencia más utilizados en aplicaciones aerodinámicas, es decir, Spallat-‐Allmaras, SST k-‐w y k-‐e (Arbeláez, 2011). Los 3 modelos mencionados hacen parte de los modelos por RANS, otro tipo de modelos son descritos en (Ferziger, 2002), donde se encuentra información acerca de modelos basados en LES (Large Eddy Simulation) y RSM (Reynolds Stress Model). Sin embargo, no se profundizará mucho en este tema, en términos teóricos, más adelante se mostrarán los criterios evaluados para utilizar uno u otro modelo.
Predicción de fuerzas Resulta evidente que una de las aplicaciones más importantes de estas herramientas computacionales es la predicción del comportamiento del fluido en estudio. En términos de aerodinámica, el diseño se basa (entre otras cosas) en la estimación de una fuerza de arrastre y otra de sustentación, parámetros necesarios para poder dimensionar otros elementos del avión como, por ejemplo, las turbinas. De forma general, las herramientas actuales son bastante acertadas prediciendo los campos de presiones que, como se mencionó anteriormente, son los que generan casi en totalidad la fuerza de sustentación. Por lo anterior se puede decir que la predicción de sustentación es bastante buena. Por el contrario, la fuerza de arrastre depende no sólo de las presiones sino también de las fuerzas de fricción. Estas fuerzas, debido al carácter turbulento del régimen en que se encuentra el vehículo dependen fuertemente del modelo de turbulencia seleccionado. La predicción de arrastre se encuentra aún lejos de ser acertada, las variaciones en predicción pueden ser hasta de 20 drag counts (dependiendo, por ejemplo, del modelo de turbulencia utilizado) mientras pruebas experimentales en túneles de viento muestran variaciones de sólo 5 (Vos, 2002).
Objetivos Objetivo general Diseñar y aplicar un protocolo para la predicción de arrastre en CFD con proyección a una futura participación de la Universidad de los Andes en el DPW, organizado por la AIAA. Objetivos específicos
1. Revisión bibliográfica de temas relacionados con la problemática de predicción de arrastre en CFD.
2. Selección de malla propuesta y utilizada en el cuarto DPW. 3. Estudio de convergencia de malla, para las condiciones propuestas por el
caso. 4. Curvas polares de arrastre para diferentes ángulos de ataque, así como
diferentes ángulos de incidencia de cola (Tail incidence angle). 5. Comparación de los resultados obtenidos con los reportados por el cuarto
DPW y por la literatura consultada.
Metodología Los procesos de predicción de arrastre por CFD cuentan con unos pasos relativamente establecidos. De acuerdo con (Dam, 1999), se debe hacer una
revisión de la geometría, selección del tamaño y la forma de la malla, una revisión del nivel de convergencia de la misma, seleccionar modelos de transición y de turbulencia. Para este caso en particular, no se necesita un modelo de transición pues se sabe que las condiciones de flujo son completamente turbulentas. Las demás etapas se describen brevemente a continuación. Esta sección del documento debe apoyarse fuertemente en el Anexo C pues se especificarán en él detalles que serán pasados por alto en este documento.
Revisión geometría La geometría utilizada para estas simulaciones fue diseñada con el fin de minimizar los errores numéricos. Es una práctica común escoger primero los modelos y el solucionador a utilizar para después generar una geometría acorde con las limitaciones de los modelos seleccionados (Seebass, 1990). Sin embargo, para este workshop, se partió de un rediseño de la geometría debido a los problemas encontrados en la geometría utilizada en las ediciones anteriores del taller (Denominada DLR-‐F6). El proceso de diseño de la nueva geometría, llamada Common Research Model se puede encontrar en la documentación del workshop. (Vassberg, 2008). Los parámetros relevantes, así como algunos comentarios acerca de la geometría, se pueden encontrar en su respectiva sección en el Anexo C.
Selección malla Como es común en este tipo de talleres donde se busca, a través del trabajo colaborativo, la comparación de resultados, las mallas utilizadas en el workshop se encuentran abiertas al público. En este orden de ideas, no se entró a generar un enmallado para la geometría a analizar, simplemente se aplicaron ciertos criterios de selección y se utilizó la que mejor se ajustaba a ellos. El enmallado de la geometría es uno de los pasos más críticos en las prácticas de CFD, debe haber un balance entre el tamaño de la malla (debe ser suficientemente fina para captar los fenómenos en la capa límite) y el costo computacional de utilizarla (entre más fina sea, mayor memoria necesita para su procesamiento). En términos de predicción de fuerzas aerodinámicas se ha estudiado la dependencia de los valores de sustentación y arrastre del tipo y el tamaño de malla. Aparentemente, la predicción de los campos de presión alrededor de la geometría son acertados sin importar mucho el refinamiento de la malla. Por el contrario, los efectos viscosos sólo se capturan de forma aceptable desde cierto tamaño de elementos sobre las superficies. La predicción de arrastre depende del tamaños de la malla, mientras la de sustentación parece no hacerlo (S. H. Peng, 2007). Teniendo esto en cuenta, se ha mostrado que la opción que mejor captura los fenómenos es esta zona es una malla estructurada en multi-‐bloque. Esta
aproximación permite refinar mucho las zonas aledañas a las superficies, para conocer el comportamiento del fluido en la capa límites, sin comprometer el tamaño de los elementos lejos de las superficies.
Figura 5. Exactitud en predicción viscosa para distintos tipos de malla. (Baker, 2005)
La Figura 5 muestra la exactitud en la predicción de la viscosidad en función de la facilidad de uso de los distintos tipos de malla presentes. Como se mencionó anteriormente, el mejor tipo de malla para estos efectos es la estructurada en multi-‐bloque, sin embargo, su facilidad de uso no es la mejor. La dificultad en el uso se refiere, en gran parte, a la dificultad en la generación de la misma, sin embargo, dado que en este caso se cuenta con las mallas ya generadas, se puede hacer caso omiso de esta restricción. Finalmente se seleccionó un juego de mallas denominado Hexa_multiblock_ansys1 más detalles sobre el juego de mallas se encuentran en la Tabla 1.
Tamaño Configuración Número de elementos Gruesa Wing-‐Body-‐Tail, iH=0º, “Gruesa” 3.5 Millones
Media
Wing-‐Body-‐Tail, iH=0º 10 Millones Wing-‐Body-‐Tail, iH=-‐2º
Wing-‐Body-‐Tail, iH=2º Wing-‐Body 8 Millones
Fina Wing-‐Body-‐Tail, iH=0º, “Fina” 35 Millones Tabla 1. Detalles juego de mallas utilizado.
Modelo de turbulencia En el mismo orden de ideas, se necesita un modelo de turbulencia que replique de forma acertada los fenómenos presentes en la capa límite. Distintas sugerencias se encontraron en la literatura, sin embargo, ninguna aproximación parecía concluyente. Lo único en que coincidían todos los autores es que la predicción de arrastre depende de la selección de este modelo (Rumsey, 2005)
1 Disponible en: ftp://cmb24.larc.nasa.gov/outgoing/DPW4/hexa_multiblock_ANSYS/
Con el fin de seleccionar adecuadamente este modelo, se corrieron 3 simulaciones, una con cada uno de los modelos tratados en (Arbeláez, 2011), con el fin de revisar convergencia en la simulación, así como tiempo de estabilización de resultados.
Valores de predicción Inicialmente, se corrieron 3 simulaciones en la malla gruesa, dejando todos los parámetros iguales, excepto el modelo de turbulencia utilizado. El valor de los coeficientes de interés, en función del número de iteraciones, se muestra en las siguientes figuras.
Figura 6. Convergencia del coeficiente de sustentación con diferentes modelos de turbulencia.
Los valores del coeficiente de sustentación (Cl) para las tres simulaciones realizadas, en función del número de iteraciones, se muestra en la Figura 6. Para el modelo S-‐A se ve como el valor aumenta drásticamente después de las 400 iteraciones, mientras que para el modelo k-‐e, después de las 500. El único valor que parece estable es el encontrado por medio del modelo SST-‐kw.
Figura 7. Convergencia del coeficiente de arrastre con diferentes modelos de turbulencia.
Al revisar el comportamiento de los distintos modelos en predicción del coeficiente de arrastre (Figura 7), las diferencias entre los 3 parecen más estrechas. De todas formas, el modelo k-‐e se aleja bastante de los otros dos modelos. S-‐A y SST k-‐w muestran un comportamiento relativamente estable después de las 150 iteraciones.
Figura 8. Convergencia del coeficiente de momento con diferentes modelos de turbulencia.
Finalmente, se revisó el comportamiento del coeficiente de momento para ambos modelos (Figura 8). Este es tal vez, de los tres coeficientes observados, el más contundente en términos de selección de modelo de turbulencia. Para los modelos S-‐A y k-‐w, los valores de este coeficiente llega a un valor muy lejano del esperado después de unas 200 iteraciones. Se puede decir que el análisis de convergencia de coeficiente sugiere la utilización del modelo SST-‐kw.
Residuales Las siguientes gráficas muestran los residuales de continuidad, de velocidad en X, Y y Z, de energía y de los parámetros k y épsilon donde es necesario. De forma general, los residuales son un indicador de las diferencias en estos valores entre un elemento y otro justo al lado del primero. Para que se pueda considerar el resultado de una simulación como adecuado, los residuales deben tener cierta estabilidad después de un número de iteraciones y deben oscilar alrededor del valor más bajo posible.
Figura 9. Residuales de simulación con modelo k-‐e
Con los residuales para cada uno de los modelos de turbulencia ensayados se refuerza la decisión que se había tomado al ver los valores a predecir con cara modelo. La Figura 9 muestra como los residuales con el modelo k-‐e se estabilizan alrededor de un valor de 1E5. Si bien muestran estabilidad, el valor alrededor del cual oscilan implica que la solución carece de sentido. La Figura 10 muestra similitudes en le comportamiento de los residuales entre el modelo Spalart-‐Allmaras y k-‐e. Los residuales parecen relativamente estables después de 100 iteraciones, sin embargo el valor numérico (1e8) sugiere nuevamente una solución cuya precisión no es aceptable. Finalmente, la Figura 11 muestra como los residuales con el modelo SST k-‐w se estabilizan en valores cercanos a 5e-‐2. Si bien el criterio de aceptación en el valor de residuales es definido por el usuario, es común tener como referencia un valor alrededor de 1e-‐3, aunque, lógicamente, lo ideal sería tener un valor de 0 (Ferziger, 2002). En este caso, se encuentran valores superiores a la referencia común, sin embargo, como se trata de un análisis comparativo, el análisis de los
residuales sugiere la utilización del modelo SST k-‐w para las simulaciones necesarias.
Figura 10. Residuales de simulación con modelo Spalart-‐Allmaras.
Figura 11. Residuales de simulación con modelo SST-‐kw.
Selección La selección del modelo de turbulencia a utilizar se llevó a cabo gracias a los resultados obtenidos en las simulaciones comentadas y los criterios analizados previamente. El análisis de resultados no es muy profundo pues la divergencia de dos de los modelos hace que la comparación sea sencilla. Del mismo modo, un valor tan alto en los residuales de las simulaciones sugiere que los resultados no
son confiables y que puede haber algún problema numérico en la implementación de ese modelo. Del análisis anterior se descarta entonces el uso de los modelos S-‐A y k-‐e. Luego, parece que el modelo que mejor desempeño tiene es el SST-‐kw. Independientemente de las simulaciones hechas, parece que este modelo es bueno prediciendo el coeficiente de sustentación, sin embargo, parece que en predicción de arrastre no es tan acertado (Arbeláez, 2011). Lo anterior parece una consecuencia del hecho de que le modelo es bueno prediciendo la zona de separación de flujo, tanto su ubicación como el comportamiento del fluido en esta región (Baker, 2005). Se utilizará entonces este modelo para las simulaciones a desarrollar.
Procedimiento El procedimiento detallado que se llevó a cabo para la obtención de resultados se encuentra descrito en el Anexo C, los parámetros utilizados, así como la configuración de los casos se encuentran descritas en este manual. De forma general se realizaron simulaciones para replicar los datos de los casos 1 y 2 y del cuarto DPW. El primer caso es un análisis de convergencia de malla mientras que el segundo es la generación de curvas polares para distintas configuraciones de la geometría del avión, variando el ángulo de incidencia del estabilizador horizontal.
Resultados y análisis
Convergencia de malla El análisis de convergencia busca verificar que la solución numérica encontrada por la simulación sea idéntica a la solución que se encontraría por medio de la ecuación diferencial, cuando el tamaño de los volúmenes de la malla tiende a cero (Ferziger, 2002). Sin embargo, dado que las ecuaciones a resolver son fuertemente no lineales, se utiliza un criterio distinto para validar la convergencia de malla. Lo que se busca es que el resultado obtenido no sea dependiente del tamaño de la malla utilizada, que la solución sea independiente y tienda a un mismo valor con cualquier tamaño de malla (Ferziger, 2002).
Validación resultados La Tabla 2 muestra los resultados obtenidos por medio de las simulaciones realizadas (columna Sim.), así como de los obtenidos del análisis estadístico del DPW 4 (Morrison J. ) y los obtenidos experimentalmente para cada uno de los coeficientes aerodinámicos analizados. En términos de sustentación (Cl), se debía garantizar un valor de 0.5 para que los resultados fueran considerados válidos, sin embargo, las simulaciones realizadas
se encuentran por fuera del rango aceptado (±0.001). El valor de Cl pudo haber sido ajustado por medio de un cambio en el ángulo de ataque de las simulaciones realizadas (2.209º, 2.26º y 2.308º para las mallas gruesas, media y fina respectivamente) pero el costo computacional era demasiado alto y se optó por dejar esta corrección como un trabajo futuro.
Para el coeficiente de arrastre (Cd), sólo el resultado de la malla gruesa se encuentra por fuera del intervalo sugerido durante el DPW, por lo que se consideran unos resultados bastante aceptables. La comparación con resultados experimentales sigue teniendo un margen bastante amplio, prueba de que la predicción de arrastre aún debe evolucionar. Finalmente el coeficiente de momento (Cm) predicho se encuentra dentro del rango de valores encontrado en el DPW pero muy lejos del experimental (Excepto en el caso de la malla fina). Lo anterior se puede explicar por el hecho de que lo valores experimentales se ven afectados fuertemente por el soporte que se utiliza para sostener la geometría en el túnel de viento. Por esta razón, el valor de Cm experimental no es muy confiable y se hace difícil determinar si la simulación se ajusta a la realidad o no.
Análisis de convergencia
Figura 12. Dependencia de coeficientes con respecto a tamaño de malla.
CL CD CM Malla Sim. DPW4 Exp. Sim. DPW4 Exp. Sim. DPW4 Exp.
Gruesa 0.4952 0.5 (±0.001)
0.5 (±0.001)
0.0284 0.02701 (±0.0081)
0.02751 (±0.0001)
-‐0.0469 -‐0.04025 (±0.0161)
-‐0.03785 (±0.0002) Media 0.4967 0.0270 -‐0.0439
Fina 0.4985 0.0268 -‐0.0379 Tabla 2. Comparación de resultados obtenidos y reportados.
Como se mencionó anteriormente, el análisis de convergencia no se puede hacer contra la solución de las ecuaciones diferenciales sino, simplemente, con respecto a la variación de los coeficientes entre una malla y otra. La Figura 12 muestra el valor de cada coeficiente aerodinámico contra el tamaño de la malla (a la izquierda más fina, a la derecha más gruesa), la gráfica muestra un comportamiento relativamente invariable, por lo que se puede hablar de convergencia. Sin embargo, un análisis más detallado es sugerido en este caso pues, como lo menciona (Salas, 2006), esta condición no es suficiente para hablar de convergencia como tal.
Estabilidad en coeficientes
Figura 13. Convergencia de Cd con respecto a iteraciones.
La Figura 13 muestra el comportamiento del valor del coeficiente de arrastre con respecto al número de iteraciones realizadas, para cada malla. En los 3 casos se puede ver que las primeras 250 iteraciones tienen grandes cambios en el valor a predecir, sin embargo, a partir de este momento, se adquiere cierta invariabilidad (menor a 5 drag counts). La malla gruesa deja de oscilar en 300 iteraciones, mientras la media se toma 450 y la fina aún oscila un poco en 500 . Las desviaciones en los 3 casos son inferiores a 3.4e-‐4, lo cual se puede considerar insignificante. Como se mencionó anteriormente, aunque esta condición no sea suficiente para hablar de convergencia de malla, se considera un comportamiento aceptable para seguir con el proceso con estas mallas.
Predicción de arrastre Entrando en materia con los valores que se buscan predecir por medio del DPW, se muestran resultados de fuerzas aerodinámicas (arrastre y sustentación), así como de desempeño aerodinámica para diferentes ángulos de ataque (entre 0º y 4º) y, un breve análisis del flujo alrededor del ala para diferentes configuraciones de estabilizador (para un ángulo de ataque de 2º).
Fuerzas aerodinámicas De forma general, los valores de sustentación para la configuración con el estabilizador a -‐2º son los más bajos para todos los ángulos de ataque probados. Luego, la configuración con el estabilizador sin ninguna inclinación (iH=0º) es la segunda peor condición de sustentación de la geometría.
Figura 14. Coeficiente de sustentación contra ángulo de ataque.
Figura 15. Coeficiente de arrastre contra ángulo de ataque.
Las otras dos configuraciones probadas (iH=2º y sin estabilizador) muestra un desempeño muy similar en términos de sustentación, siendo las que mejor coeficiente presentan para todo el barrido de ángulos ensayados (Figura 14).
Figura 16. Curvas polares para diferentes configuraciones.
La Figura 15 muestra la tendencia esperada para el coeficiente de arrastre con respecto al ángulo de ataque. Dada la relación cuadrática entre el arrastre y la sustentación (Bertin, 2002), el crecimiento del arrastre debe ser similar a una parábola, como se muestra en la gráfica. A diferencia de la sustentación, el arrastre no tiene una configuración que sea mejor o peor que otras en todos los ángulos probados. A bajos ángulos de ataque (entre 0º y 1º) parece que la mejor es la configuración iH=-‐2º, pero a altos ángulos de ataque (entre 2º y 4º) el comportamiento es mejor para la configuración iH=2º.
Desempeño aerodinámico El desempeño aerodinámico del vehículo utilizado en estas simulaciones se evaluará por medio de las curvas polares obtenidas (producto del segundo caso de estudio del DPW), así como de la eficiencia aerodinámicas de cada configuración. Las curvas polares de las diferentes configuraciones simuladas se muestran en la Figura 16. De forma general, se espera que la mejor configuración se encuentre hacia la izquierda y en la parte superior de la gráfica pues en esta zona el arrastre es mínimo y la sustentación máxima, que corresponde a la configuración sin estabilizador. La configuración más ineficiente se encuentra por lo tanto en la esquina opuesta (inferior, derecha) que corresponde a iH=-‐2º. Estas curvas polares se pueden obtener por varios modelos que relacionan Cl con Cd de forma cuadrática, sin embargo, se deben determinar ciertos factores en situaciones que no fueron simuladas (por ejemplo Cl=0). Por lo anterior, no se hace una comparación con alguno de los modelos teóricos mencionados. (Hull, 2007).
Figura 17. Eficiencia aerodinámica para diferentes configuraciones
La eficiencia aerodinámica para cada una de las configuraciones analizadas, en función del ángulo de ataque, se muestra en la Figura 17. Con el fin de tener un orden de magnitud en la escala vertical de la figura, se debe mencionar que el valor de eficiencia aerodinámica para un vehículo aéreo de transporte con una velocidad de crucero cercana a 0.8 Mach es de 18 (Hale, 1984). Este valor teórico se alcanza en la configuración con iH=0º, lo cual corrobora que los datos obtenidos por la simulación no son muy lejanos a los reales. Resulta interesante notar que el valor de la eficiencia es más alto, en todo el rango de ángulos de ataque, para la configuración sin estabilizador trasero que para cualquier inclinación del mismo, llegando a valores de 20. Por el contrario, configuraciones como inclinaciones de -‐2º en el estabilizador disminuyen esta eficiencia a un valor máximo de 17, afectando considerablemente factores como el consumo de combustible de una aeronave. Esta figura de mérito es tal vez la que mejor caracteriza la aerodinámica de la geometría analizada y permite suponer que el papel del estabilizados trasero no es más que el de un elemento de control de vuelo. En términos aerodinámicos, los vehículos aéreos serían más eficientes sin este tipo de dispositivos, sin embargo, el control de vuelo necesita de un elemento que permita, por ejemplo, variar el ángulo de ataque del vehículo.
Condiciones de flujo Debido a la eficiencia mostrada en la Figura 17, las condiciones de flujo se analizaron únicamente para el vuelo a 2º de ángulo de ataque. La cantidad de datos simulados permite analizar hasta 28 condiciones de flujo distintas (una por punto en las gráficas de la sección anterior), sin embargo, se tomaron únicamente las 4 diferentes posiciones del estabilizador, al ángulo de ataque mencionado.
Del mismo modo, se analiza únicamente el flujo en la zona del estabilizador pues es el único parámetro que varía de una condición de vuelo a la otra. El flujo en la parte delantera del fuselaje o en el ala del vehículo se asume igual para todos los casos, las diferencias se deben encontrar en la zona alterada, es decir, el estabilizador. Para validar esta suposición se compara el arrastre generado por componentes, el resultado se muestra en la Tabla 3. Configuración Elemento iH=0º iH=+2º iH=-‐2º No Tail Fuselaje 0.0097 0.0097 0.0098 0.0096 Ala 0.0146 0.0146 0.0145 0.0145 Estabilizador 0.0008 0.0020 0.0007 NA
Tabla 3. Valores de Cd por elementos.
Los datos muestran claramente que el fuselaje y el ala del avión no presentan cambios en términos de arrastre (apenas 2 drag counts para el primero y 1 para el segundo) al variar la posición del estabilizador. Esto corrobora el supuesto hecho y verifica que el flujo, corriente arriba del estabilizados, es el mismo en todas las simulaciones que comparten el mismo ángulo de ataque. Por el contrario, el aporte del estabilizador varía hasta en 13 drag counts. Dado que le interés del DPW es el arrastre, el análisis de flujo se hará en términos de presiones y fricción en la zona mencionada inicialmente, en caso de ser necesario se evaluarán otras posibles fuentes de arrastre durante el análisis de resultados
Coeficiente de fricción de pared La Figura 18 muestra los contornos de coeficiente de fricción en el estabilizador para las 4 configuraciones simuladas, tanto en la parte superior como inferior del estabilizador. La parte superior del estabilizador parece no mostrar grandes variaciones en el coeficiente, todas muestras una zona de alta fricción en el borde de ataque, como es de esperarse, que se va disminuyendo a medida que se avanza con respecto a la cuerda del estabilizador. Hacia el borde de salida los valores son tan bajos que se puede suponer el inicio de una zona de desprendimiento, sin embargo los datos no parecen concluyentes en esta sección. A diferencia del estabilizados, el fuselaje muestra claramente una zona de separación donde el coeficiente tiene un valor de 0. La zona es muy similar en la zona cercana al borde de ataque, para las diferentes configuraciones, sin embargo, a cerca al borde de salida, la separación es más grande en la configuración iH=+2º mientras que la configuración iH=-‐2º muestra la menos separación de las 3. La configuración sin estabilizador no muestra esta zona mencionada. La parte inferior muestra un comportamiento similar, con una zona de alta fricción en el borde de ataque y la disminución de este valor a lo largo de la cuerda, sin embargo, se encuentra un comportamiento ligeramente diferente
para la configuración iH=-‐2º, donde una zona de baja fricción aparece entre el borde de ataque y la zona media del estabilizador. Más adelante se buscará una explicación para este comportamiento. El comportamiento en la parte baja del fuselaje es similar en todas las configuraciones, hay una zona de separación (marcada en azul) de aproximadamente el mismo tamaño en todas las configuraciones.
Coeficiente de presiones Los coeficientes de presiones para las diferentes posiciones del estabilizador se muestran en la Figura 19. La condición de referencia debe ser en la que el estabilizador no tiene ninguna inclinación particular (iH=0º), en este caso se ve el comportamiento esperado para las presiones sobre un ala en condiciones transónicas. El borde de ataque presenta una línea roja de alta presión y, justo después, se tiene una zona azul donde el coeficiente es negativo, es decir, se genera una zona de “succión” (Davari & Soltani, 2011) que se extiende por la parte central del estabilizador. La parte inferior del estabilizador muestra un contorno azul más oscuro, es decir, la presión que se genera sobre esta zona es mayor, sugiriendo un fenómenos de reattachment de flujo, más adelante se verificará si se trata de esto o no. La siguiente configuración (iH=+2º) muestra un comportamiento similar, sin embargo, la coloración del estabilizador muestra que la zona de succión se extiende más hacia el fuselaje. Así mismo, el hecho de que el coeficiente de presión sea más grande (en magnitud) implica que el fenómenos es más intense que en la configuración anterior. Para la cara inferior del estabilizador por el contrario, la zona que se veía en la configuración anterior se recoge hacia el borde de ataque, afectando una menor porción de la superficie. El coeficiente de presión del estabilizador con un ángulo de incidencia de –2º muestra un comportamiento opuesto al mencionado anteriormente. La zona superior muestra una distribución de presiones relativamente uniforme en todo el estabilizador (a excepción del borde de ataque, donde es más alta). Sin embargo, la cara interior del estabilizador muestra que la zona de succión se corrió hacia la parte central del estabilizador, dejando una zona de presión moderada entre esta zona y el borde de ataque. En todos los casos se evidencia una pequeña zona de alta presión en la zona más alejada del fuselaje, en la cara superior del estabilizador. Esta zona se le puede atribuir al efecto denóminado downwash, el paso del estabilizador por el flujo genera un vórtices alrededor de la punta del estabilizador que recae en dirección normal a la superficie, aumentando la presión allí. El efecto es más notable en la primera configuración (iH=0º), en la segunda parece ser absorbido por algo de mayor magnitud y en la última, es apenas perceptible.
Vista Superior Inferior
iH=0º
iH=+2º
iH=-‐2º
No Tail
Figura 18. Coeficiente de fricción para diferentes posiciones del estabilizados
Superior Inferior
iH=0º
iH=+2º
iH=-‐2º
No Tail
Figura 19. Coeficiente de presiones para diferentes posiciones del estabilizador.
La configuración sin estabilizador muestra una distribución de presiones normal, sólo se diferencia por una zona de presión menor a las otras configuraciones en donde se une con el borde de ataque del estabilizador. Con el fin de comprender mejor el fenómenos asociado a las zonas azules de los contornos, se hicieron cortes a 25%, 50%, 75% y 95% de la cuerda, a partir del fuselaje y se graficaron los valores del coeficiente de presión en toda la superficie del estabilizador. La escala de vertical se encuentra invertida con el fin de ver la superficie superior en la parte alta de la gráfica y la inferior, en la baja. Las gráficas se muestran en la Figura 20. Las distribuciones de presión en la superficie del estabilizador muestra que su comportamiento hacia la zona del borde de salida es relativamente constante en todas las configuraciones y en todos los planos analizados. Sin embargo, las Figura 18 y Figura 19 dejan ver una aparente zona de separación de flujo en la zona cercana al fuselaje. Esta zona no se evidencia en la Figura 20 porque incluso el plano más cercano al fuselaje se encuentra muy alejado para capturar el fenómeno. Por lo demás, se puede corroborar que todas las secciones tienen un máximo de presión en aproximadamente x/b=0.1. Al comparar esto con los contornos mostrados antes, se corrobora que esta en esta región es que se presentan las zonas de alta presión (azules). De hecho, las gráficas correspondientes a la configuración iH=+2º muestran las presiones más altas y una zona donde parece estabilizarse en el primer cuarto del estabilizador. Este comportamiento es típico de un ala en flujo transónico, donde se encuentra una onda de choque producto del paso de la velocidad del flujo por la barrera del sonido. Además, en las curvas que representan la zona inferior se encuentran regiones donde se pasa el valor de -‐0.5, suficiente para verse en color azul en los contornos. Lo anterior sucede en los cuadros A y B para la configuración iH=0º y en el cuadro C para la misma configuración así como para la configuración iH=-‐2º. En este caso no es evidente si de trata de una zona de presión debido a la onda de choque identificada o a separación de flujo, por lo que se busca ver la onda de choque para entender el fenómeno encontrado.
Ondas de choque Para estudiar la onda de choque encontrada y evaluar la influencia de la separación de flujo en las distribuciones de presiones encontradas, se hizo un corte en la mitad del estabilizador (x/b=0.5) y se generó un contorno con el número de MACH y algunas líneas de corriente para visualizar el flujo. La Figura 21 muestra el resultado del proceso mencionado, para las 3 configuraciones con estabilizador.
A
B
C
D
Figura 20. Coeficientes de presión en diferentes planos.
Inicialmente se puede decir , por las líneas de flujo, que en esta sección del estabilizador no se genera ninguna separación de flujo considerable, los cambios de presión considerables parecen ser producto entonces de condiciones de vuelo transónicos, donde hay zonas de aceleración de flujo. Al evaluar la velocidad del flujo cerca al perfil, se puede ver claramente que hay zonas donde la velocidad es superior a la del sonido que corresponden a todas aquellas que no son de color verde. Sin embargo, la velocidad del flujo cambia considerablemente de una configuración a la otra. La configuración inicial (iH =0º) muestra un flujo que aún se puede considerar con velocidad ligeramente inferior a la del sonido, en la parte superior, sin embargo, la parte inferior del perfil muestra una pequeña onda de choque (M=1.2), que corresponde con el aumento de presión visto en los contorno. La configuración correspondiente a iH=-‐2º muestra un comportamiento muy similar a la anterior, la zona superior muestra aceleración del flujo de forma moderada, mientras la zona inferior presenta una aceleración mucho más importante que la mostrada anteriormente en tamaño y en intensidad (M>1.3).
Las dos configuraciones mencionadas son las que muestran zonas de presiones azules en la parte inferior del estabilizador.
iH=0º
iH=-‐2º
iH=+2º
Figura 21. Número de mach y líneas de corriente.
Finalmente, la otra configuración analizada (iH=+2º) muestra una pequeña zona de flujo a M=1.1 en su parte superior. Esta zona se puede considerar responsable de que esta sea la única configuración que muestra zonas azules en la parte superior a aproximadamente x/b=0.25. La zona inferior también muestra una pequeña zona donde el flujo se acelera a la velocidad del sonido, sin embargo no es tan importante como la de las configuraciones anteriores.
Control aerodinámico La Figura 22 muestra el valor del coeficiente de momento para cada configuración a diferentes ángulos de ataque. Esta gráfica es útil a la hora de determinar las condiciones de vuelo a las que el estabilizado neutraliza los momentos generado por la aeronave, para esto se busca el corte de la línea con el eje horizontal.
Al evaluar las curvas obtenidas, se puede ver que la única condición que cruza el eje horizontal de la gráfica es en la que el estabilizador no tiene inclinación alguna (iH=0º). Como es de esperarse, un ángulo de incidencia negativo aumenta el coeficiente de momento y uno positivo lo disminuye, sin embargo, ninguna de estas condiciones genera estabilidad estática en el rango de valores simulados.
Configuración iH=0º iH=-‐2º iH=+2º No Tail
Ángulo de trim [º] 1.3 4.1* -‐1.4* -‐2* Tabla 4. Angulo de trim para diferentes configuraciones.
Figura 22. Coeficiente de momento contra ángulo de ataque.
Al evaluar las curvas obtenidas, se puede ver que la única condición que cruza el eje horizontal de la gráfica es en la que el estabilizador no tiene inclinación alguna (iH=0º). Como es de esperarse, un ángulo de incidencia negativo aumenta el coeficiente de momento y uno positivo lo disminuye, sin embargo, ninguna de estas condiciones genera estabilidad estática en el rango de valores simulados. Dado que la pendiente de las curvas es relativamente estable entre 0º y 3º de ángulo de ataque, se utilizó esta pendiente para estimar un valor de corte, en la zona de ángulos de ataque negativos. Los resultados encontrados se muestran en la Tabla 4, los valores estimados tienen un asterisco (*). Los valores encontrados sugieren entonces que, en condicione de vuelo estables, sin ninguna perturbación externa, las configuraciones con el estabilizador a 0º y a -‐2º tienden a generar un vuelo ascendente a 1.3º y 4.1º ángulos de ataque respectivamente. Por el contrario, si el estabilizador de ubica a 2º de inclinación
o no se coloca estabilizador, el vuelo tiende a ser descendente, con ángulos de ataque de -‐1.4º y -‐2º respectivamente.
Costo computacional La Tabla 5 muestra la memoria y el número de CPU’s utilizados para las simulaciones del primer caso de estudio. Del mismo modo, se muestra una aproximación del tiempo utilizado. La Tabla 6 muestra lo mismo, para el segundo caso de estudio.
Malla Gruesa Media Fina Memoria [GB] 132 132 264 CPU's 8 8 32 Tiempo aproximado (500 iteraciones) [h] 5.5 12.5 20
Tabla 5. Costo computacional de primer caso de estudio.2
De este caso es necesario destacar el requerimiento en memoria para procesar la malla “Fina”, los 132 GB utilizados previamente no fueron suficientes por lo que fue necesario utilizar dos máquinas con esa configuración. Los tiempos de simulación entre la “Gruesa” y la “Media” parecen comportarse de manera lineal con el tamaños de las mallas (ver Tabla 1). No es posible comparar el tiempo necesario para la malla “Fina” pues se utilizaron muchas más unidades de procesamiento. Dado que para este caso, las mallas de todas las configuraciones son de tamaño similar (ver Tabla 1), bajo las mismas condiciones de número de procesadores, los tiempos de simulación son muy similares. Se puede ver como, al multiplicar los procesadores por 4, se tiene una reducción del tiempo de simulación a aproximadamente una cuarta parte. Esto comprueba que existe una relación lineal entre el número de procesadores utilizados y el tiempo de cómputo. En términos de memoria, todas las configuraciones se simularon con 132 GB sin ningún problema.
Malla iH=0º iH=+2º iH=-‐2º Sin estabilizador Memoria [GB] 132 132 132 132 CPU's 32 8 8 6
Tiempo aproximado (500 iteraciones) [h]
2.5 12 13.5 14.5
Tabla 6. Costo computacional segundo caso de estudio.2
2 Los tiempos se estimaron de las simulaciones realizadas en el servidor master-‐hpc-‐mox de la Universidad de los Andes. Las unidades de procesamiento y la memoria también corresponden a las disponibles en este servidor.
Conclusiones El desarrollo de este proyecto cumplió con los objetivos planteados inicialmente, replicando los datos obtenidos en el 4th DPW por medio de CFD, dando un análisis un poco más profundo de los visto en las presentaciones de este evento. Inicialmente, se verificó que la selección de la malla realizada cumplía con los requisitos de convergencia necesarios para llevar a cabo las simulaciones posteriores. Al mismo tiempo, se desarrolló y documentó una metodología para la ejecución de procesamiento en ANSYS FLUENT. En términos de resultados, se encontraron los comportamientos esperados de arrastre y sustentación al variar el ángulo de ataque. Además, se observó que la presencia de un estabilizador no genera ningún beneficio desde el punto de vista de desempeño aerodinámico. De hecho, la configuración que no tiene estabilizador genera una eficiencia aerodinámica por encima de las reportadas como “normales” para un vehículo de este tipo. El flujo se analizó en detalle únicamente en la condición de vuelo donde la eficiencia aerodinámica era máxima (AoA=2º), sin embargo, el análisis hecho se puede extender a otras condiciones de vuelo en trabajos futuros. Para esta condición de vuelo se encontró que el arrastre generado en el ala y el fuselaje era el mismo para todas las configuraciones, mientras el estabilizador mostraba cambios considerables (13 drag counts). Los contornos de fricción y de presión sobre el estabilizador mostraron que el cambio en fricción era mínimo, por lo que no se puede atribuir un gran cambio al arrastre por este componente. Por el contrario, se encontraron varios cambios en las presiones a lo largo de la superficie del estabilizador. Las gráficas de presión a los largo de la superficie, para cuatro planos de corte sobre el elemento mostraron comportamiento que sugerían la presencia de una onda de choque. El número de Mach en el flujo cercano al estabilizador y algunas líneas de corriente evidenciaron que el arrastre adicional es producto de una condición de vuelo transónica, más que de un fenómeno de separación de flujo. El arrastre adicional encontrado se puede atribuir al conocido como wave drag. Finalmente, dado que el estabilizador parece un elemento dedicado únicamente al control de la aeronave, se encontró que la única condición que garantiza un ángulo de estabilidad en el rango de ángulos de ataque simulado es la condición normal, sin ninguna inclinación. Cualquier variación que se haga en el estabilizador hace imposible mantener el vehículo con una orientación fija en este rango de ángulos de ataque, por el momento que genera. Dado el gran número de simulaciones realizadas, diferentes análisis se pueden sugerir como trabajo futuro. La presencia de la onda de choque y su influencia en el arrastre para ángulos de ataque distintos al utilizado parece la primera opción
para seguir el trabajo de este proyecto. De todas formas, el análisis de las zonas de separación cerca al fuselaje, tanto en el estabilizador con en el ala del vehículo sería otro análisis de interés. Sería interesante también simular condiciones de donde el coeficiente de sustentación sea 0, para comparar el comportamiento aerodinámico del vehículo con alguno de los modelos teóricos disponibles en la literatura. En términos de control, validar las estimaciones de estabilidad realizadas, simulando condiciones de vuelo con ángulos de ataque negativos podría complementar el trabajo realizado en esta dirección.
Referencias ANSYS, F. (2009). ANSYS FLUENT 12.0 Theory Guide. Arbeláez, D. (2011). Comparación de 3 modelos computacionales de turbulencia en aplicaciones aerodinámicas. Baker, T. J. (2005). Mesh generation: Art or science? Progress in aerospace science , 29-‐63. Bertin, J. J. (2002). Aerodynamics for engineers. Estados Unidos: Prentice Hall. Dam, C. V. (1999). Recent experience with different methods of drag prediction. Progress in Aerospace Sciences . Davari, A., & Soltani, M. R. (2011). Effects of wing geometry on wing-‐body-‐tail interference in subsonic flow. Scientia Iranica . Ferziger, J. H. (2002). Computational methods for fluid dynamics. Nueva York: Springer. Frink, N. T. (2006). 3rd AIAA CFD Drag Prediction Workshop. San Francisco: AIAA. Hale, F. J. (1984). Introduction to aircraft performance, selection and design. Estados Unidos: John Wiley & Sons, Inc. Hoerner. (1965). Fluid-‐dynamic drag. Estados Unidos. Hull, D. (2007). Fundamentals of Airplane Flight Mechanics. Berlin: Springer Berlin Heidelberg. Moran, J. (1984). An introduction to theoretical and computational aerodynamics. Nueva York: John Wiley & Sons, Inc. Morrison, H. J. 4th AIAA CFD Drag Prediction Workshop. San Antonio: AIAA. Morrison, J. (n.d.). Statistical Analysis of CFD Solutions from the Fourth AIAA Drag Prediction Workshop. Rumsey, C. L. (2005). Study of CFD variation on transport configuration from the second drag-‐prediction workshop. Computer & Fluids , 785-‐816. S. H. Peng, P. E. (2007). Some remarks on CFD drag prediction of an aircraft model. Proceedings of the fifth international conference on fluid dynamics . Salas, M. D. (2006). Some observations on grid convergence. Computers & Fluids , 688–692. Seebass, A. R. (1990). Applied computational aerodynamics. Washington: AIAA. Spalart, P. (2010). Reflections on RANS modelling. In P. Spalart, Progress in Hybrid RANS-‐LES Modelling (pp. 7-‐24). Berlin: Springer Berlin Heidelberg. Torenbeek, E. (2009). Flight Physics. Netherlands: Springer Netherlands. Vassberg, J. (2008). Development of a Common Research Model for Applied CFD Validation Studies. Velandia, J. S. (2013). Manual para configuración de casos para el cuarto Drag Prediction Workshop. Vos, J. B. (2002). Navier-‐Stokes solvers in European aircraft design. Progress in Aerospace Sciences , 601-‐697. White, F. Fluid Mechanics. McGraw-‐Hill. Yechout, T. (2003). Introduction to aircraft flight mechanics. Reston: AIAA.
Anexos
Anexo A: Archivos de MATLAB. Lectura de archivos. %% Lectura de archivos. % Se importa el archivo de texto que contiene el historial % de los coeficientes. Arch=importdata('NOMBRE DE ARCHIVO',' '); % Se calculan el numero de lineas que tiene. Este valor debe limitarse % manualmente cuando el archivo contiene valores para varios ?ngulos de % ataque. max=length(Arch.data(:,1)); % Recorrido para sumar el aporte de todas las superficies. for i=1:max Coef(i)=sum(Arch.data(i,2:end)); end % Generaci?n de vector de iteraciones. xo=[1:max]; % Generaci?n de gr?fica de coeficiente contra iteraciones plot(xo,Ct) Comparación modelo de turbulencia. %% Analisis modelo turbulencia % Se importan los archivos con los valores de los coeficientes ClSA=importdata('clAoA0SA.txt',' ',2); CmSA=importdata('cmAoA0SA.txt',' ',2); CdSST=importdata('cdAoA0SST.txt',' ',2); ClSST=importdata('clAoA0SST.txt',' ',2); CmSST=importdata('cmAoA0SST.txt',' ',2); CdKE=importdata('cdAoA0SSTke.txt',' ',2); ClKE=importdata('clAoA0SSTke.txt',' ',2); CmKE=importdata('cmAoA0SSTke.txt',' ',2); % Se adecuan los valores importantes para comparaci?n CdSA=CdSA.data(:,2); ClSA=ClSA.data(:,2); CmSA=33.667/7.005*CmSA.data(:,2); CdSST=CdSST.data(:,2); ClSST=ClSST.data(:,2); CmSST=33.667/7.005*CmSST.data(:,2); CdKE=CdKE.data(:,2); ClKE=ClKE.data(:,2); CmKE=33.667/7.005*CmKE.data(:,2); % Se genera el vector de iteraciones its1=[1:500]; % Se generan 3 gr?cias, una para cada coeficiente con los resultados % obetnidos por cada modelos de turbulencia. figure() plot(its1,ClSA,its1,ClSST,its2,ClKE(1:700),'LineWidth',2) grid on
xlabel('Iteraciones','FontSize',13,'FontWeight','bold') ylabel('Cl','FontSize',13,'FontWeight','bold') title('Convergencia Cl','FontSize',13,'FontWeight','bold') legend('S-A','SST-kw','k-e','Location','NorthWest') axis([0 600 -0.3 0.6]) figure(2) plot(its1,CdSA,its1,CdSST,its2,CdKE(1:700),'LineWidth',2) grid on xlabel('Iteraciones','FontSize',13,'FontWeight','bold') ylabel('Cd','FontSize',13,'FontWeight','bold') title('Convergencia Cd','FontSize',13,'FontWeight','bold') legend('S-A','SST-kw','k-e','Location','NorthWest') axis([0 600 -0.1 0.4]) figure(3) plot(its1,CmSA,its1,CmSST,its2,CmKE(1:700),'LineWidth',2) grid on xlabel('Iteraciones','FontSize',13,'FontWeight','bold') ylabel('Cm','FontSize',13,'FontWeight','bold') title('Convergencia Cm','FontSize',13,'FontWeight','bold') legend('S-A','SST-kw','k-e','Location','SouthWest') axis([0 600 -2 0.5]) Generación de gráficas. %% Generaci?n de gr?ficas % Se utiliz? el archivo de lectura de archivos para obtener los % datos en los siguientes arreglo. % AoA representa el ?ngulo de ataque, ClIHX representa la configuraci?n, % reemplazando X por 0, m2 (menos 2) o M2 (mas2) % Adem?s se cacula la eficiencia aerodin?mica EfHX, con la misma % nomenclatura AoA=[0,1,1.5,2,2.5,3,4]; ClIH0=[0.165,0.3087,0.3805,0.4565,0.5331,0.5961,0.6715]; CdIH0=[0.01814,0.02036,0.02246,0.02509,0.0292,0.0353,0.0508]; CmIH0=[0.0642,0.01351,-0.0105,-0.03215,-0.0549,-0.066,-0.0487]; EfH0=ClIH0./CdIH0; ClIHm2=[0.1284,0.2699,0.3429,0.4186,0.4951,0.5592,0.6360]; CdIHm2=[0.0198,0.0211,0.0228,0.0250,0.0287,0.0345,0.0494]; CmIHm2=[0.2019,0.1543,0.1298,0.1080,0.0857,0.0720,0.0868]; EfHm2=ClIHm2./CdIHm2; ClnT=[0.2070,0.3390,0.4082,0.4822,0.5549,0.6130,0.6825]; CdnT=[0.0160,0.0195,0.0214,0.0240,0.0282,0.0343,0.0497]; CmnT=[-0.0847,-0.1024,-0.1121,-0.1189,-0.1291,-0.1297,-0.1030]; EfnY=ClnT./CdnT; ClIHM2=[0.2018,0.3448,0.4183,0.4923,0.5667,0.6268,0.7009]; CdIHM2=[0.0182,0.0209,0.0229,0.0263,0.0306,0.0369,0.0525]; CmIHM2=[-0.0753,-0.1270,-0.1478,-0.1662,-0.1849,-0.1925,-0.1688]; EfHM2=ClIHM2./CdIHM2; % El archivo genera 5 gr?ficas, las curvas polares, Cl vs AoA, Cd vs Aoa, % eficiencia aerodin?mica vs AoA y Cm vs AoA.
figure(1) plot(CdIH0,ClIH0,'o-',CdIHm2,ClIHm2,'o-',CdnT,ClnT,'o-',CdIHM2,ClIHM2,'o-') grid on xlabel('Cd','FontSize',16,'FontWeight','bold') ylabel('Cl','FontSize',16,'FontWeight','bold') title('Curva Pola','FontSize',16,'FontWeight','bold') legend('iH=0','iH=-2','noTail','iH=+2','Location','NorthWest') figure(2) plot(AoA,ClIH0,'o-',AoA,ClIHm2,'o-',AoA,ClnT,'o-',AoA,ClIHM2,'o-') grid on xlabel('AoA','FontSize',16,'FontWeight','bold') ylabel('Cl','FontSize',16,'FontWeight','bold') title('Cl Vs AoA','FontSize',16,'FontWeight','bold') legend('iH=0','iH=-2','noTail','iH=+2','Location','NorthWest') figure(3) plot(AoA,CmIH0,'o-',AoA,CmIHm2,'o-',AoA,CmnT,'o-',AoA,CmIHM2,'o-') grid on xlabel('AoA','FontSize',16,'FontWeight','bold') ylabel('Cm','FontSize',16,'FontWeight','bold') title('Cm Vs AoA','FontSize',16,'FontWeight','bold') legend('iH=0','iH=-2','noTail','iH=+2') figure(4) plot(AoA,EfH0,'o-',AoA,EfHm2,'o-',AoA,EfnY,'o-',AoA,EfHM2,'o-') grid on xlabel('AoA [o]','FontSize',16,'FontWeight','bold') ylabel('Cl/Cd','FontSize',16,'FontWeight','bold') title('Eficiencia aerodinamica','FontSize',16,'FontWeight','bold') legend('iH=0','iH=-2','noTail','iH=+2') figure(5) plot(AoA,CdIH0,'o-',AoA,CdIHm2,'o-',AoA,CdnT,'o-',AoA,CdIHM2,'o-') grid on xlabel('AoA','FontSize',16,'FontWeight','bold') ylabel('Cd','FontSize',16,'FontWeight','bold') title('Cd Vs AoA','FontSize',16,'FontWeight','bold') legend('iH=0','iH=-2','noTail','iH=+2','Location','NorthWest') Convergencia de malla %% Convergencia de malla %Lectura de archivos CdG=importdata('cd-historyGAoA.txt',' ',2); CdM=importdata('cd-historyMAoA.txt',' ',2); CdG=CdG.data(:,2); CdM=CdM.data(:,2); CdM(end)=[]; ClG=importdata('cl-historyGAoA.txt',' ',2); ClM=importdata('cl-historyMAoA.txt',' ',2); ClG=ClG.data(:,2); ClM=ClM.data(:,2); ClM(end)=[]; CmG=importdata('cm-historyGAoA.txt',' ',2); CmM=importdata('cm-historyMAoA.txt',' ',2);
CmG=CmG.data(:,2); CmM=CmM.data(:,2); CmM(end)=[]; % Generaci?n gr?fica valores contra iteraciones its1=[1:500]; plot(its1,CdG,'-',its1,CdM,'-',its1,CdF,'-') grid on xlabel('Iteraciones','FontSize',13,'FontWeight','bold') ylabel('Cd','FontSize',13,'FontWeight','bold') title('Convergencia de malla','FontSize',13,'FontWeight','bold') legend('Gruesa','Media','Fina') axis([50 500 0.025 0.07]) % Tama?os de malla y valores para cada valor en cada una. xTam=[3516705,10793559,35159816]; xTam=xTam.^(-2/3); yCd=[CdG(end),CdM(end),mean(CdF(end-50))]; yCm=[CmG(end)*33.667/7.005,CmM(end),mean(CmF(end-50))*33.667/7.005]; yCl=[ClG(end),ClM(end),mean(ClF(end-50))]; % Generaci?n gr?fica valores contra n?mero de elementos. plot(xTam,yCd,'-o',xTam,yCm,'-o',xTam,yCl,'-o') grid on xlabel('Tamano malla (N\^(-2/3))','FontSize',13,'FontWeight','bold') ylabel('Coeficiente','FontSize',13,'FontWeight','bold') title('Analisis convergencia de malla','FontSize',13,'FontWeight','bold') legend('Cd','Cm','Cl','Location','West')
Anexo B: Ejemplo journal utilizado.
El siguiente cuadro muestra un ejemplo de journal para la configuración de monitores y condiciones de flujo a un ángulo de ataque. Si se necesitan simular varios ángulos (o alterar algún otro parámetro), uno después del otro (como fue el caso en este proyecto), basta con copiar el texto y alterar lo que sea necesario (componentes vectoriales y nombres de archivos). define bc pff far no 0 no 0.85 no 310.93 yes no 0.999390827 no 0 no 0.034899497 no no yes 5 10 solve monitors force sdm cdAoA2 yes aircraft_body_fairing aircraft_body_fuse aircraft_body_taill aircraft_body_tailte aircraft_body_tailu aircraft_body_wingl aircraft_body_wingte aircraft_body_wingu () no yes cdAOA2 no yes 0.999390827 0 0.034899497 solve monitors force slm clAoA2 yes aircraft_body_fairing aircraft_body_fuse aircraft_body_taill aircraft_body_tailte aircraft_body_tailu aircraft_body_wingl aircraft_body_wingte aircraft_body_wingu () no yes clAOA2 no yes -‐0.034899497 0 0.999390827 solve monitors force smm cmAoA2 yes aircraft_body_fairing aircraft_body_fuse aircraft_body_taill aircraft_body_tailte aircraft_body_tailu aircraft_body_wingl aircraft_body_wingte aircraft_body_wingu () no yes cmAOA2 no yes 33.6779 0 4.51993 0 1 0 solve iterate 1000 wcd AoA2
Anexo C. Manual de configuración
Tabla de contenidos
Introducción. ...................................................................................................................... 3 Nomenclatura. .................................................................................................................... 3
Descripción del caso. ....................................................................................................... 3 Caso de estudio 1 : Convergencia de malla. ....................................................................... 4 Caso de estudio 2 : Estudio de Downwash .......................................................................... 4
Geometría. ........................................................................................................................... 5
Configuración ..................................................................................................................... 6 Enmallado ..................................................................................................................................... 6 Archivos base ............................................................................................................................................ 6 Adecuación y corrección. ...................................................................................................................... 7
Solver .............................................................................................................................................. 8 Modelos .......................................................................................................................................... 8 Materiales ...................................................................................................................................... 9 Condiciones de operación ..................................................................................................... 10 Condiciones de frontera ........................................................................................................ 10 Valores de referencia ............................................................................................................. 12 Métodos de solución ............................................................................................................... 13 Otras configuraciones ............................................................................................................ 13
Metodología ...................................................................................................................... 14 Transferencia de archivos .................................................................................................... 14 Simulación .................................................................................................................................. 14 Ejecutar FLUENT ................................................................................................................................... 14 Cargar caso .............................................................................................................................................. 15 Modificar condiciones de frontera ................................................................................................ 15 Modificar monitores ............................................................................................................................ 16 Inicializar flujo ....................................................................................................................................... 16 Simular ...................................................................................................................................................... 17 Puntos de control y archivo de datos ........................................................................................... 17
Interpretación de resultados ............................................................................................... 17 Referencias ........................................................................................................................ 18
Lista de Figuras Figura 1. Geometría utilizada en las simulaciones. ............................................................... 5 Figura 2. Cuadro “Mesh” en la “Configuración general” de FLUENT. ............................ 7 Figura 3. Resultado de la revisión de la malla. ........................................................................ 7 Figura 4. Diálogo "Scale Mesh". ..................................................................................................... 8 Figura 5. Configuración del solver. .............................................................................................. 8 Figura 6. Configuración de modelos. ........................................................................................... 9 Figura 7. Configuración modelo de turbulencia. .................................................................... 9 Figura 8. Propiedades del aire. ................................................................................................... 10 Figura 9. Zonas para condiciones de frontera. ..................................................................... 11 Figura 10. Configuración pressure-‐far-‐field. ........................................................................ 12 Figura 11. Valores de referencia. ............................................................................................... 12 Figura 12. Métodos de solución. ................................................................................................ 13 Figura 13. Coeficiente de sustentación contra iteraciones. ............................................ 18
Lista de Tablas Tabla 1. Condiciones de vuelo generales. (Morrison) ........................................................... 3 Tabla 2. Condiciones de simulación caso 1. (Morrison) ....................................................... 4 Tabla 3. Condiciones de simulación caso 2. (Morrison) ....................................................... 5 Tabla 4. Parámetros geométricos relevantes. ......................................................................... 6 Tabla 5. Sugerencias mallas a utilizar. ........................................................................................ 6 Tabla 6. Condiciones de frontera en zonas. ........................................................................... 11
Introducción. En Junio del 2009 se llevó a cabo el cuarto taller de predicción de arrastre (Drag Predictions Workshop, DPW) organizado por el Instituto Americano de Astronáutica y Aeronáutica (AIAA) en San Antonio, Texas. El objetivo de estos talleres es conocer el estado del arte de la predicción de arrastre por medio de dinámica de fluidos computacional (CFD) buscando identificar las áreas donde se necesita profundizar en términos de investigación y desarrollo. Con el fin de preparar la participación de la Universidad de los Andes en las próximas ediciones, se escribe este manual de configuración. Se replicaron las simulaciones hechas por el grupo de Ansys FLUENT durante su participación en el cuarto DPW, por medio del proceso descrito a lo largo de este documento.
Nomenclatura.
-‐ AIAA Instituto Americano de Astronáutica y Aeronáutica -‐ Cl Coeficiente de sustentación -‐ Cd Coeficiente de arrastre -‐ Cm Coeficiente de momento -‐ Cref Cuerda de referencia -‐ CRM Common Research Model -‐ DLR Centro aeroespacial alemán. -‐ DPW Drag Prediction Workshop -‐ iH Ángulo inclinación de cola (estabilizador) -‐ Re Número de Reynolds -‐ Sref Área de referencia -‐ SST Shear stress transport -‐ Xref Centro de momento en X -‐ Yref Centro de momento en Y -‐ Zref Centro de momento en Z
Descripción del caso. El cuarto DPW busca predecir el arrastre en una geometría similar a la de un avión comercial en condiciones de vuelo de crucero. Una breve descripción de la geometría se hará en la siguiente sección. En cuanto a las condiciones de vuelo, por parámetros relevantes se muestran en la Tabla 1.
Parámetro Valor Mach 0.85 M Re 5E6 Temperatura de referencia 100º F
Tabla 1. Condiciones de vuelo generales. (Morrison)
Las condiciones mencionadas anteriormente son las que se deben cumplir a lo largo de las simulaciones a realizar, de todas formas, se desprenden dos casos de estudio con ligeras variaciones, a continuación se explican ambos.
Caso de estudio 1 : Convergencia de malla. Como es común en la práctica de simulaciones por métodos numéricos, es necesario realizar un análisis de convergencia de malla. El análisis busca corroborar que los resultados obtenidos por las simulaciones son similares, independiente del tamaño de la malla utilizada. Con lo anterior se espera encontrar un rango de tamaño de malla en que los resultados sigan siendo invariables o independientes y el costo computacional sea lo más bajo posible. Además de las condiciones mencionadas en la Tabla 1, este análisis de convergencia se llevará a cabo bajo la siguiente configuración:
Parámetro Valor Cl 0.500 (±0.001) iH 0º
Tabla 2. Condiciones de simulación caso 1. (Morrison)
Las condiciones de vuelo deben ser simuladas para, por lo menos, tres mallas distintas, denominadas “Gruesa”, “Media” y “Fina”, más adelante se entrará en detalle acerca de las mallas utilizadas para este caso. Finalmente, lo que se busca con este estudio es generar una gráfica de convergencia del parámetro de interés (en este caso, el coeficiente de arrastre) contra el número de iteraciones para validar la invariabilidad del valor encontrado. Del mismo modo, se comparará este valor contra el número de elementos de la malla, buscando el número de elementos adecuado para realizar el siguiente caso de estudio.
Caso de estudio 2 : Estudio de Downwash Este caso de estudio busca comparar las condiciones de arrastre para diferentes condiciones de vuelo de la geometría analizada. Su nombre se debe al fenómeno del flujo que se genera entre la zona de alta y baja presión alrededor del ala, generando un flujo de aire en dirección contraria a la de la fuerza de sustentación. El objetivo de este caso de estudio es generar cuatro curvas polares a las condiciones mostradas en la Tabla 3. Una vez se tengan las curvas polares, se podrá analizar las diferencias en términos de sustentación y arrastre para cada una de las condiciones simuladas.
Elemento Ángulos de variación (º) Cola horizontal -2, 0, 2, Sin elemento. Fuselaje-Ala 0, 1, 1.5, 2, 2.5, 3, 4. Tabla 3. Condiciones de simulación caso 2. (Morrison)
Las condiciones de vuelo son similares a las mostradas en la Tabla 1. Este estudio se llevará a cabo con la malla denominada “Media”, de todas formas, esto será sólo una referencia del tamaño de la misma pues, debido a las diferencias en los casos de cada curva polar, se deben usar enmallados distintos para cada una.
Geometría. Hasta el DPW anterior (tercera edición (Frink, 2006)) la geometría que se había utilizado era un modelo conocido como DLR-‐F6. Su nombre hace referencia al instituto que lo desarrolló (Centro Aeroespacial Alemán, DLR por sus siglas en alemán) y al perfil alar utilizado (F6). Sin embargo, esa geometría generaba algunos problemas numéricos en las soluciones, por lo que se necesitaba un cambio. Con el fin de tener un modelo común tanto para prácticas de CFD como para pruebas en diferentes túneles de viento, buscando hacer comparables resultados tanto numéricos como experimentales, se generó un nuevo modelo para la cuarta edición del DPW. El modelo se denominó Common Research Model (CRM), diseñado por un grupo de colaboración del centro de investigación de Langley de la NASA y la compañía Boeing. La geometría generada es fiel a un avión comercial, desde el fuselaje hasta las dimensiones necesarias en los estabilizadores para efectos de control y maniobrabilidad. La Figura 1 muestra la geometría mencionada.
Figura 1. Geometría utilizada en las simulaciones.
Los archivos para generar la geometría (y sus diferentes configuraciones) se encuentran disponibles en la página del DPW (Morrison) Para efectos de este manual, sólo nos interesarán algunas propiedades geométricas del modelo, listadas en la Tabla 4.
Parámetro Valor Sref 594720 in2 (383.69 m2) Cref 275.8 in (7.005 m) Xref 1325.9 in (33.68 m) Yref 468.75 in (12.42 m) Zref 177.95 in (4.52 m) Tabla 4. Parámetros geométricos relevantes.
Configuración
Enmallado
Archivos base Las directivas del DPW publicaron una serie de recomendaciones para la generación de las mallas a utilizar, no se entrará muy en detalle en las recomendaciones. 1 De forma general, las mallas deben ser semi-‐esferas cuyo diámetro sea equivalente a aproximadamente 100 cuerdas del modelo. Los tamaños sugeridos de las mallas, así como las diferentes configuraciones, se muestran en la Tabla 5.
Tamaño Configuración Número de elementos Gruesa Wing-‐Body-‐Tail, iH=0º, “Gruesa” 3.5 Millones
Media
Wing-‐Body-‐Tail, iH=0º 10 Millones Wing-‐Body-‐Tail, iH=-‐2º
Wing-‐Body-‐Tail, iH=2º Wing-‐Body 8 Millones
Fina Wing-‐Body-‐Tail, iH=0º, “Fina” 35 Millones Tabla 5. Sugerencias mallas a utilizar.
Debido al carácter colaborativo de estos eventos, las mallas utilizadas por todos los participantes se encuentran disponibles en la página del sitio (Morrison). Después de revisar el contenido disponible, se encontró que había un juego de mallas utilizado directamente en ANSYS Fluent, denominado hexa_multiblock_ANSYS. Este juego de mallas cuenta con los tamaños y configuraciones sugeridos por los organizadores (Tabla 5), es necesario mencionar que las mallas fueron generadas en pulgadas, por lo que deberán ser escaladas antes de usarlas. Cada una de las mallas necesarias se puede descargar en un archivo comprimido que cuenta con formatos .tin (geometría), .blk (bloques de generación de malla) y .cgns (“Sistema de notación general en CFD”). Una vez se obtienen los archivos es necesario pasar las mallas por un proceso de adecuación antes de configurar y simular los diferentes casos.
1 Los detalles se pueden consultar en: http://aaac.larc.nasa.gov/tsab/cfdlarc/aiaa-‐dpw/Workshop4/gridding_guidelines_4.html 2 TUI hace referencia a la Interfaz de Texto de Usuario, el cuadro de comandos a la derecha en la
Adecuación y corrección. Para el juego de mallas mencionado, se pueden importar directamente los archivos .CGNS a Ansys FLUENT y trabajar sobre ellos escalándolos a las dimensiones necesarias y corrigiendo los errores encontrados. Una vez se abre FLUENT y se le asigna el directorio de trabajo adecuado para encontrar los archivos descargados, se debe importar el archivo .CGNS por medio del menú: FILE > IMPORT > CGNS > MESH Se selecciona el archivo mencionado en la ventana de diálogo y se esperan unos minutos a que el programa lea la malla importada.
Figura 2. Cuadro “Mesh” en la “Configuración general” de FLUENT.
Una vez se importó la malla, es necesario revisarla, por medio del botón “Check” que se ve en la Figura 2. Al correr esta verificación se encontrará un error por algunas (la cantidad depende de la malla que se esté utilizando) “Left-‐handed faces”, el procedimiento para arreglar este problema, aparece en la ventana de comandos y se muestra en la Figura 3.
Figura 3. Resultado de la revisión de la malla.2
Luego, teniendo en cuenta que la generación de las mallas fue hecha en pulgadas, es necesario escalarlas por medio del botón “Scale” en la Figura 2. En el diálogo que se abre es necesario utilizar, en la parte izquierda el panel “Scaling”, indicando que la malla fue creada en pulgadas. El panel mencionado, así como las dimensiones finales necesarias se muestran en la Figura 4.
2 TUI hace referencia a la Interfaz de Texto de Usuario, el cuadro de comandos a la derecha en la visualización por defecto de FLUENT.
Figura 4. Diálogo "Scale Mesh".
Una vez se llevan a cabo estos dos procesos (escalado y corrección de “left-‐handed faces”) se puede comenzar con la configuración del modelo.
Solver Inicialmente se debe seleccionar el tipo de solucionador que se va a utilizar, para este caso se utilizará uno basado en densidad. Además, como no se va a simular nada en estado transitorio, el tiempo se tomará estacionario. El cuadro de configuración se muestra en la Figura 5.
Figura 5. Configuración del solver.
Modelos Para estas simulaciones es necesario activar dos modelos, el de energía y el de viscosidad, los demás deben estar apagados, como se muestra en la Figura 6.
Figura 6. Configuración de modelos.
El modelo de energía basta con encenderlo mientras que para el de viscosidad hay que escoger qué tipo de modelo se va a usar, así como definir sus parámetros. La comparación de varios modelos para esta aplicación específica sugirió que el que mejor garantizaba la convergencia de la simulación es el modelo SST k-‐omega, sus parámetros se muestran en la Figura 7.
Materiales Dentro de los materiales que se deben definir sólo es necesario tratar con el fluido Air, definido por FLUENT. Sin embargo, es necesario ajustar sobretodo la propiedad de viscosidad pues, para lograr el Re necesario en la simulación, se debe asumir que la viscosidad del aire es constante en el valor mostrado en la Figura 8, donde también se pueden ver las otras propiedades relevantes.
Figura 7. Configuración modelo de turbulencia.
Figura 8. Propiedades del aire.
Condiciones de operación Dentro de las condiciones de operación (en el panel “Cell Zone Conditions”) el único parámetro relevante a configurar es la Presión de operación, que debe ser igual a 1 atmósfera (101325 Pa).
Condiciones de frontera Las condiciones de frontera representan una de las partes más importantes de cualquier modelo a simular. En este caso, se cuenta con 4 condiciones para distintos tipos de zona.
Figura 9. Zonas para condiciones de frontera.
Las condiciones de frontera para cada una de las zonas (Figura 9) presentes en las mallas utilizadas se muestran en la Tabla 6
Zona Condición de frontera aircraft_body_XX Wall far pressure-‐far-‐field int_fluid Interior symm Symmetry
Tabla 6. Condiciones de frontera en zonas.
En particular, la configuración de la zona far va a ser importante pues es la zona donde se variará las condiciones de flujo, permitiendo cambios en el ángulo de ataque del avión. Su configuración se muestra en la Figura 10. Los valores de X-‐Component of Flow Direction y Z-‐Component of Flow Direction son los únicos que se van a variar a lo largo de las simulaciones. Las condiciones de turbulencia, el número de Mach y las condiciones térmicas (en la pestaña Thermal, temperatura de 310.93 K) serás las mismas para todos los casos.
Figura 10. Configuración pressure-‐far-‐field.
Valores de referencia La sección de los valores de referencia se puede leer directamente de la zona far, de todas formas se muestran los valores utilizados en la Figura 11.
Figura 11. Valores de referencia.
Además de los valores leídos desde la zona indicada, se tienen valores geométricos que vienen del modelo utilizado. El área y la longitud característica
se toman de la Tabla 4. El área se divide entre 2 pues se está simulando únicamente la mitad de la geometría en cada caso.
Métodos de solución Distintos tipos de solución pueden ser configurados en ANSYS Fluent, los utilizados para estos casos se muestran en la Figura 12.
Figura 12. Métodos de solución.
Además de las opciones mostradas, es necesario activar la opción “Convergence Acceleration for Stretched Meshes” en este panel de configuración.
Otras configuraciones Las demás secciones de configuración del caso no tienen grandes cambios con respecto a la opciones por defecto. El panel Monitors se configurará durante la metodología de simulación pues se utilizará un monitor diferentes para cada condición de flujo. El panel Initialization también se manejará durante la metodología, sin embargo, es necesario dejar activada la opción de Compute from far y Relative to Cell Zone. El panel Solution Controls se deja por defecto, así como el de Dynamic Mesh.
Metodología Esta sección busca explicar los pasos para llevar a cabo las simulaciones, una vez se tenga configurado el caso como se describió en la sección anterior. Se debe tener en cuenta que, los recursos computacionales de la Universidad de los Andes obligan a hacer la configuración en un servidor distinto al que correrá las simulaciones como tal. Por lo anterior, se debe hacer un proceso de transferencia de archivos, uno de simulación y otro de interpretación de resultados.
Transferencia de archivos La transferencia de archivos de un servidor a otro se puede hacer por medio de un cliente SCP cualquiera (Por ejemplo, WinSCP). Las credenciales que brinda el MOX para conectarse funcionan con este tipo de transferencia. El proceso es sencillo:
-‐ Descargar el archivo .cas guardado en la máquina donde se configuró el caso.
-‐ Cargar el archivo .cas en el servidor donde se va a simular el caso. Es necesario tener presente la ubicación donde se subió el archivo del caso configurado, pues desde esta misma carpeta se debe ejecutar FLUENT.
Simulación Los comandos que se mostrarán después de la sección de cómo ejecutar FLUENT se pueden colocar juntos en un archivo de texto y utilizarlo como journal. Luego se puede ejecutar el programa con este archivo y se llevarán a cabo la secuencia de comandos indicada, uno después del otro. Un ejemplo de este tipo de archivos se muestra en el Anexo B. Para tener acceso al servidor se debe iniciar una sesión por medio de SSH con algún cliente de este tipo, por ejemplo, putty.
Ejecutar FLUENT Los servidores utilizados por la Universidad cuentan con un sistema operativo Linux, por lo que la ejecución del programa se debe hacer por medio de comandos. Para ejecutar FLUENT hay que dirigirse a la carpeta donde se guardó el archivo .cas por medio del comando cd. Se debe escribir, en la ventana de putty donde se inició la sesión ssh el siguiente comando: cd ruta_del_directorio_de_trabajo En este punto, se puede ejecutar FLUENT, por medio del siguiente comando:
/share/apps/ansys/v140/fluent/bin/fluent -‐g 3ddp -‐tX -‐pinfiniband -‐cnf=compute-‐0-‐Y -‐ssh En el comando anterior, X debe ser reemplazado por el número de procesadores que se va a utilizar, debe ser concertado con los otros usuarios del servidor en ese momento, con el fin de dar buen manejar al número de licencias. El número debe variar entre 4 y 32 en condiciones normales. La Y se debe reemplazar, a la fecha, por 1 o 0, se refiere al nodo que se busca usar en el servidor. En principio no implica ningún problema hacerlo en uno o en otro, sin embargo se debe gestionar los recursos computacionales del servidor de tal forma que el trabajo se reparta en ambos nodos para evitar fallos de memoria. En este momento, se abre FLUENT y se encuentra en el equivalente a la TUI cuando hay interfaz gráfica.
Cargar caso Inicialmente se carga el caso que se configuró y se colocó en el directorio de trabajo. El comando para hacerlo es el siguiente. file read-‐case nombre_caso.cas Durante el proceso se reporta la lectura de la malla, así como la construcción de la misma, luego aparece el símbolo > en la pantalla, indicando que ya se cargó el caso y se pueden ejecutar más comandos.
Modificar condiciones de frontera Para variar el ángulo de ataque en cada simulación, se debe modificar el pressure-‐far-‐field. El comando a ejecutar se muestra a continuación. define bc pff far no 0 no 0.85 no 310.93 yes no 1 no 0 no 0 no no yes 5 10 De izquierda a derecha, los números que se encuentran en el comando representan:
-‐ Presión relativa. Se coloca en 0 pues la de operación se configuró en 1 atmósfera.
-‐ Número de Mach: Para todos los casos se utilizará 0.85 -‐ Temperatura: Para todos los casos se utilizarán 310.93 K -‐ Componente X del flujo: Corresponde al coseno del ángulo de ataque. -‐ Componente Y del flujo: Corresponde 0 en todos los casos. -‐ Componente Z del flujo: Corresponde a menos el seno del ángulo de
ataque.
-‐ Parámetros de turbulencia: 5 y 10 corresponden a la intensidad y la componente viscosa del modelo de turbulencia configurado. (Ver Figura 10)
Para cambiar el ángulo de ataque basta con alterar los valores de los componentes del flujo en el comando mostrado y volver a ejecutar la simulación.
Modificar monitores Los monitores son los valores de interés en la simulación, se debe configurar uno para el coeficiente de sustentación, otro para el coeficiente de arrastre y otro para el coeficiente de momento. solve monitors force sdm cdX yes lista_de_zonas () yes yes cdY no yes 1 0 0 solve monitors force slm clX yes lista_de_zonas () yes yes clY no yes 1 0 0 solve monitors force smm cmX yes lista_de_zonas () yes yes cmY no yes 33.6779 0 4.51993 0 1 0 En el comando, cdX, clX y cmX son el nombre que se le pone al monitor, se recomiendo que X sea reemplazado por el ángulo de ataque utilizado, para mayor facilidad en la manipulación de resultados. Del mismo modo, cdY, clY y cmY son los nombres de los archivos de texto donde se guardarán los datos, por lo que se recomienda la misma práctica. La frase lista_de_zonas hace referencia a las zonas donde se desea conocer el valor de drag, lift y momento, para la este caso se debe reemplazar lista_de_zonas por todas las zonas mostradas en la Figura 9 que se configuraron con condición de frontera tipo wall. Se deben escribir una zona después de la otra, separadas por un punto. Cualquier error en el nombre de estas zonas hace que el comando sea mal ejecutado y, por lo tanto, los monitores estén mal definidos. Los últimos 3 números de las líneas del comando que definen drag y lift (1 0 0) representan nuevamente las componentes en X, Y y Z del flujo que se está trabajando, por lo que deben ser reemplazados por el coseno del ángulo, 0 y el seno del ángulo respectivamente. La configuración del monitor de momento es la misma para todos los casos, no se debe modificar nada en función del flujo.
Inicializar flujo En este punto se han configurado todos los parámetros de forma adecuada, por lo que se debe inicializar el flujo. Debido al tipo de malla utilizado y los demás parámetros de configuración, no basta con inicializar el flujo para garantizar
convergencia en los resultados. Se deben realizar dos inicializaciones, una normal y una conocida como Fast-‐Multi-‐Grid (FMG), el comando es el siguiente. solve i if solve i fmg-‐i yes Los parámetros de la incialización FMG no se alterarán, el último “yes” del comando simplemente permite que se inicialice de esta forma con los parámetros por defecto. Este puede ser el paso más delicado antes de comenzar a simular el caso, puede tomar varios minutos.
Simular Una vez se inicializa el flujo, el proceso de simulación es realmente sencillo, se ejecuta el comando mostrado abajo y se espera que se ejecuten las iteraciones indicadas. solve iterate 1000 En este caso se colocan 1000 iteraciones pues los casos del DPW muestran cierta estabilidad en la solución después de este número. De todas formas, para algunos casos parece que 700 son necesarias y el ahorro computacional es considerable, mientras que para otras parecen necesarias hasta 1500. Una práctica común es correr un número de iteraciones menor e ir revisando el resultado para decidir si se simulan más iteraciones, hasta encontrar la estabilidad deseada.
Puntos de control y archivo de datos Otra práctica común es generar archivos de datos en un número de iteraciones intermedias, como una medida preventiva en caso de perder conexión con el servidor o tener algún otro error durante la simulación. El comando para guardar estos archivos de datos es el siguiente. wd nombre_archivo.dat Después se puede cargar este archivo (junto con el .cas que se generó) para continuar las iteraciones adicional necesarias, por medio del comando mostrado en la sección anterior. El mismo comando mostrado en esta sección es utilizado para guardar el archivo de datos con los resultados de la simulación.
Interpretación de resultados Además de los archivos de datos generado, son de especial importancia los datos guardados por los monitores. Dependiendo de lo que se desee visualizar, se debe
acceder a unos o a otros archivos. En este caso, el interés está sobre todo en los valores de arrastre y sustentación encontrados. Debido a la gran cantidad de datos que se pudieron haber registrado, es recomendable tratarlos por medio de un programa tipo MATLAB, en vez de una hoja de cálculo. El archivo de texto guardado reporta el valor de arrastre o sustentación para cada una de las zonas elegidas, por lo que se deben sumar todas las columnas dentro de una misma zona para encontrar el valor total del coeficiente. Al graficas, por ejemplo, el valor de sustentación con respecto al número de iteraciones, de las simulaciones realizadas para una de las curvas polares del caso de estudio 2, se encuentra algo de este estilo (Figura 13).
Figura 13. Coeficiente de sustentación contra iteraciones.3
Por la forma en que se llevaron a cabo las simulaciones, es recomendable promediar el valor de las últimas 100 iteraciones de cada monitor para estimar su valor. Sin embargo, una inspección gráfica de la oscilación de cada uno es recomendable con el fin de ajustar más el número de iteraciones a promediar, en una zona de poca oscilación.
Referencias ANSYS, F. (2009). ANSYS FLUENT 12.0 Theory Guide. Frink, N. T. (2006). 3rd AIAA CFD Drag Prediction Workshop. San Francisco: AIAA. Morrison, H. J. 4th AIAA CFD Drag Prediction Workshop. San Antonio: AIAA.
3 Los valores en la gráfica no se pueden tomar como referencia pues representan el valor en una inclinación dada, no variando conforme al ángulo de ataque.