Optimización de formas aerodinámicas

Plataforma de optimización de formas aerodinámicas.

Como parte del tema de optimización del curso de "Modelado, Simulación y Optimización de Sistemas Aeroespaciales" [1], se desarrolló una plataforma en Scilab, la cual está enfocada a optimización de geometrías aerodinámicas y vincula diversos paquetes computacionales. Se diseñó de manera modular para permitir emplear distintos algoritmos de optimización, y usa la simulación de la dinámica de fluidos para la estimación de aptitud de forma aerodinámica. El objetivo de optimización puede ser tanto el coeficiente de sustentación como el de arrastre.

La plataforma integra diversos paquetes computacionales de libre uso, los cuales son un software de diseño mecánico asistido por computadora - FreeCAD, un software de dinámica de fluidos - OpenFOAM, y un software de cálculo científico – Scilab. Los paquetes anteriores tienen una interfaz de línea de comandos, que les permite ser utilizados en la optimización automática sin intervención humana.

La plataforma funciona como se muestra en la Figura 1. El programa comienza con la inicialización del caso de OpenFOAM, el cual es una estructura de archivos, mostrada en la Figura 2, que contienen toda la información necesaria para realizar la simulación CFD.

Figura 1: Diagrama de flujo del programa de optimización.

Figura 2: Estructura de archivos de OpenFOAM.

La carpeta 0.orig, contiene los archivos necesarios para definir las condiciones iniciales. Los archivos k, nut, omega, p y U, definen las condiciones iniciales para la energía cinemática turbulenta, la viscosidad, la frecuencia de turbulencia, la presión y la velocidad del flujo respectivamente. Dentro de la carpeta include, se tiene el archivo fixedInlet, que define que el flujo a la entrada será constante (fixedValue). El archivo wallsPatches, le indica al software que las paredes del volumen, tendrán condiciones de deslizamiento (walls type slip). El archivo initialConditions define los valores de las condiciones iniciales.

La siguiente carpeta, constant, tiene los archivos transportProperties y turbulenceProperties, generalmente empleados para casos que no incluyen la ecuación de la energía. En el primero se define el modelo de la viscosidad, mismo que puede ser Newtoniano, Bird-Carreau, Ley Potencial, Hershel-Bulkey, entre otros. En el segundo archivo se especifica el modelo para la turbulencia, que puede ser laminar (sin turbulencia), RAS (Raynolds Average Simulation), o LES (Large-Eddy Simulation).

Dentro de la carpeta triSurface es en donde está la geometría creada, en formato STL (STereoLithography).

Finalmente, la carpeta system, contiene los diccionarios para la generación de mallas (blockMeshDict y snappyHexMeshDict), para el cálculo de los coeficientes (forceCoeffs), y otros archivos y diccionarios que se analizarán más adelante.

En el archivo blockMeshDict, también se define la entrada del flujo (inlet), la salida del flujo (outlet), y el conjunto de paredes: frente, fondo, suelo, techo (walls). Estas partes conforman las fronteras del caso.

Después de la inicialización del caso, el programa entra en el ciclo principal de optimización. En este ciclo, se realiza lo siguiente:

1. El programa genera nuevos valores para los parámetros de geometría a optimizar y produce un archivo de comandos en lenguaje Python, que es utilizado por el programa FreeCAD para crear la geometría del objeto de optimización. FreeCAD ejecuta la secuencia de comandos de Python para generar y guardar en el disco un archivo de la geometría del objeto en formato STL.
2. El programa ejecuta las utilidades blockMesh y snappyHexMesh de OpenFOAM para crear y mallar el volumen circundante a la geometría del objeto, que fungirá como túnel de viento virtual.
3. El programa usa el solucionador CFD llamado simpleFoam para resolver el caso.
4. El programa lee el archivo con los resultados del cálculo de los coeficientes aerodinámicos y toma el último valor del coeficiente de arrastre, cuyo valor mínimo es el objetivo de optimización.
5. El programa verifica la condición de parada, que consta de una tolerancia mínima para el coeficiente de arrastre, o bien, de un máximo de iteraciones. Si la condición de parada se cumple, el programa guarda la geometría óptima en el disco, de lo contrario, el programa repite el ciclo de optimización.

Actualmente, la plataforma utiliza dos diferentes algoritmos de optimización: algoritmo genético y sección dorada. El método de sección dorada es un algoritmo simple para optimizar funciones unimodales y unidimensionales. La ventaja principal de este algoritmo es la facilidad de implementación y su velocidad de trabajo, ya que no requiere muchas evaluaciones de la función objetivo. El algoritmo genético no tiene limitaciones en la dimensión del espacio de búsqueda y la forma de la función, sin embargo es complejo y lento.

Experimento

Para probar el funcionamiento de la plataforma, se realizó un ejemplo de optimización de forma aerodinámica, tomando como base el Problema Aerodinámico de Newton Restringido a Conos Truncados y utilizando el software de CFD. La solución de este problema de optimización con el uso de OpenFOAM como un estimador, se realizó en tres etapas.

Primera etapa – Simulación de un caso

En la primera etapa se estableció un caso de OpenFOAM para la simulación CFD de flujo de aire (se tomó el aire como medio circundante) alrededor de un cono truncado. En particular, la geometría de un cono truncado en el formato STL se creó manualmente utilizando FreeCAD, como se muestra en la Figura 3.

Figura 3: Geometría del cono desarrollado en CAD en formato STL.

Luego, en el archivo blockMeshDict, se establecieron los límites del área de simulación (la dimensión del túnel de viento virtual), así como también la entrada del flujo (inlet), la salida del flujo (outlet), y el conjunto de paredes (walls): frente, fondo, suelo y techo. Las fronteras se representan como se muestre en la Figura 4.

Figura 4: Fronteras del caso.

Posteriormente se definieron las condiciones iniciales como se muestran en la Tabla 1.

Tabla 1: Condiciones iniciales y de frontera.

Una vez definidas las fronteras y establecidas las condiciones iniciales y de frontera para el caso, se procede a la simulación. Para ello, se recure al archivo Allrun. Este archivo, contiene la lista de instrucciones que serán ejecutadas secuencialmente en OpenFOAM (Figura 5).

Figura 5: Estructura del archivo Allrun.

El comando surfaceFeatureExtract, extrae y escribe en un nuevo archivo las propiedades de la superficie creada anteriormente (Figura 3).

El comando blockMesh, es un generador de malla multi-bloque. Esto genera una malla muy burda, como se muestra en la Figura 6.

Figura 6: Discretización del volumen con malla multi-bloque.

Después se utiliza el comando snappyHexMesh, el cual refina la malla y la ajusta a la superficie, para obtener una geometría más suave y llana, y con ello una discretización más apropiada del dominio, como se muestra en la Figura 7. A este proceso se le conoce como refinamiento de malla. Cabe recordar que mientras más fina sea la malla (mas pequeños los volúmenes discretizados), mejor será la aproximación del cálculo, sin embargo también será mayor el costo computacional, pues se requerirá un mayor número de cálculos.

Figura 7: Discretización del volumen con malla refinada.

El comando patchSummary, escribe los campos y las condiciones de frontera para cada contorno. Posteriormente, el comando potentialFoam, ejecuta un solucionador de flujo potencial, con el cual se generan los valores iniciales que serán utilizados después por solucinadores más complejos. Finalmente, simpleFoam es el solucionador de estado estable para flujo incompresible y turbulento, que utiliza el algoritmo SIMPLE (proceso numérico para resolver las Ecuaciones de Navier-Stokes, cuyas siglas provienen de Semi-Implicit Method for Pressure Linked Equations). Este solucionador utiliza los archivos contenidos en la carpeta system como directivas de solución. El archivo decomposeParDict, fragmenta el caso de forma que se pueda correr en paralelo (cuando se cuenta con más de un procesador). En el archivo forceCoeffs, se seleccionan los coeficientes que serán evaluados, tales como coeficiente de sustentación y coeficiente de arrastre entre otros, siendo este último el que se tomará en cuenta como parámetro de solución para evaluarse en el programa de optimización.

El diccionario fvSchemes, establece los esquemas numéricos que se utilizarán para calcular términos como las derivadas, gradientes o interpolaciones que haya dentro de las ecuaciones que se resuelvan durante la simulación. Entre ellos, los siguientes esquemas (o subdirectorios).

1. Esquema temporal.

   Se puede configurar para resolver problemas de estado transitorio o estado estable. En este caso, se configura para estado estable, lo cual implica igualar las derivadas temporales a 0.

2. Esquema de gradiente.

   Configura los términos de gradientes. En este caso, como “Gauss lineal”. Lo cual especifica la discretización de volumen finita estándar de la integración gaussiana, que requiere de una interpolación de valores desde centros celulares a centros frontales. El esquema de interpolación utilizado será lineal o de diferenciación central.

3. Esquema de divergencia.

   Incluye términos de divergencia, incluyendo términos de advección y otros términos difusivos.

4. Esquema de gradientes normales a superficies.

   Contiene términos de gradiente normal de superficie, los cuales se evalúan en las caras de las celdas, y son los componentes normales a ellas. La base del cálculo del gradiente en una cara, es restar el valor en el centro de la celda de un lado de la cara, del valor en el centro del el otro lado, y dividir por la distancia. El cálculo es de precisión de segundo orden para el gradiente normal a la cara, si el vector que conecta los centros celulares es ortogonal a la cara, es decir, están en ángulo recto.

5. Esquema Laplaciano.

   Este subdirectorio contiene los términos laplacianos como el de difusión en las ecuaciones de momento. El esquema de Gauss es la única opción de discretización y requiere una selección de un esquema de interpolación para el coeficiente de difusión. En todos los casos, el esquema de interpolación lineal se usa para la interpolación de la difusividad.

6. Esquema de interpolación.

   Están contenidos los términos que son interpolaciones de valores típicamente desde centros celulares a centros frontales, usados principalmente en la interpolación de centros de velocidad a cara para el cálculo del flujo. Típicamente se emplea la interpolación lineal, salvo en algunos casos particulares como el análisis de tensión, donde se usa la interpolación cúbica.

Los esquemas seleccionados para este caso particular son los que se muestran en la Tabla 2.

Tabla 2: Esquemas numéricos utilizados.

El diccionario fvSolution, controla los solucionadores de ecuaciones lineales, tolerancias y los algoritmos utilizados en la solución. Puede o no contener a su vez un sub-diccionario para el análisis de estrés. En este archivo se especifican cada uno de los solucionadores lineales que se usan para cada ecuación discreta. Estos solucionadores distinguen la simetría de las matrices (ya sean simétricas o asimétricas), cuya simetría depende de los términos de la ecuación a resolver. Por ejemplo las derivadas temporales y términos laplacianos forman matrices simétricas, mientras que las derivadas advectivas forman matrices asimétricas.

Entre algunos de los solucionadores están el PCG, que es para gradientes precondicionados estabilizados bi-conjugados, tanto para matrices simétricas como asimétricas; el GAMG, para multi-mallas geométrico-algebráica generalizadas; o el Diagonal, para solucionar sistemas explícitos diagonales. Los solucionadores, tolerancias y suavisadores utilizados para este caso se muestran en la Tabla 3.

Tabla 3: Solucionadores utilizados para el caso.

Tras correr el comando simpleFoam, se realiza el cálculo numérico correspondiente y se simula el caso. Posteriormente, se obtienen los resultados de la simulación CFD, que son: distribución de presión sobre el cono (Figura 8), velocidad del aire alrededor del cono (Figura 9) y contornos de velocidad (Figura 10). El parámetro de interés para este caso, es la distribución de presión, pues está directamente relacionada con el coeficiente de arrastre.

Figura 8: Distribución de presión sobre el cono.

Figura 9: Distribución de presión sobre el cono y velocidad del aire alrededor del cono.

Figura 10: Contornos de flujo de aire alrededor del cono.

Como parte del post-procesamiento, se genera un archivo con los coeficientes de momento, arrastre y sustentación. Para el caso del coeficiente de arrastre, los resultados se muestran en la Figura 11. El programa toma el valor del coeficiente de arrastre de la iteración 100, pues como se puede ver en la figura, en ese punto el valor del coeficiente se mantiene estable.

Figura 11: Coeficiente de arrastre vs Iteración.

Con base en estos resultados, se pueden hacer comparativas de otras geometrías simuladas bajo exactamente las mismas condiciones, lo cual permite justificar la selección de una geometría sobre otra dependiendo de las cualidades y características buscadas.

Sin embargo, comparar manualmente cada geometría resulta un proceso laborioso y repetitivo. Por tal motivo, se planteó el desarrollo de la plataforma de optimización, que genere las geometrías, analice cada una de ellas y seleccione la más adecuada, empleando el software CFD como estimador.

Segunda etapa – Simulación de 10 casos

En la segunda etapa, se desarrolló un programa en Scilab que escanea el intervalo de búsqueda y modifica automáticamente el radio del cono (desde 0, para formar un cono con punta, hasta el radio de la base, para formar un cilindro) con un paso determinado. En la Figura 12 se pueden ver tres de diez distintas variaciones de la geometría del cono trunco.

Figura 12: Distintas geometrías dentro del intervalo de búsqueda.

Posteriormente, el programa ejecuta la simulación de CFD bajo exactamente las mismas condiciones, y calcula el coeficiente de arrastre para cada una de las diez variaciones de la geometría. El resultado del escaneo del radio se muestra en la Figura 13.

Figura 13: Coeficiente de arrastre vs radio de la parte trunca.

Como se puede ver en la Figura 13, al variar el radio de la parte trunca del cono, también irá variando el coeficiente de arrastre. De las 10 geometrías generadas y simuladas, la que opone menor resistencia al flujo del aire al desplazarse a velocidad constante, es el cono con un radio de la parte trunca cercano a 1.05 m.

Los resultados obtenidos en la segunda etapa distan del resultado analítico del problema aerodinámico de Newton (mínimo de la función mostrada en la Figura 13). Esto se debe a que el problema simulado y el problema Aerodinámico de Newton son esencialmente distintos, ya que este último considera choques perfectamente elásticos entre las partículas del medio circundante y el sólido. En contraparte, el caso simulado contempla aire como medio circundante y lo trata como un medio continuo, sin partículas. Estas diferencias afectan las condiciones de frontera del problema y con ello el resultado. Estas diferencias afectan las condiciones de frontera del problema y con ello el resultado.

Tercera etapa – Optimización de la geometría

Finalmente, en la tercera etapa, se optimizó la forma aerodinámica del cono truncado utilizando la plataforma descrita anteriormente. Inicialmente, se utilizó la función estándar del algoritmo genético de Scilab como algoritmo de optimización, sin embargo, el tiempo de ejecución de optimización excedió 24 horas. Por lo tanto, se realizó la optimización con el algoritmo de sección dorada, con el cual se obtuvieron resultados en menos de 1 hora de simulación.

En el proceso de optimización, el programa encontró en 17 iteraciones que el radio de la parte trunca del cono que genera el mínimo arrastre es de 1.0439 m, con una tolerancia de 0.00068. Al comparar los resultados obtenidos de la segunda y tercera etapa del experimento, se puede concluir que la plataforma de optimización funciona correctamente.

En la Figura 14 se muestra el coeficiente de arrastre en relación a las iteraciones del programa de optimización, y se puede ver claramente que el programa converge a un valor de CD = 0.2604, el cual corresponde al cono con radio de la parte trunca de 1.0439 m.

Figura 14: Cambio del coeficiente de arrastre contra iteraciones de optimización.

Referencias

1. Kovalenko, Y., Escobedo, R. O., Prokhorov, Y. (2018). Software de código abierto para un curso de Modelado, Simulación y Optimización de Sistemas Aeroespaciales.
https://somim.org.mx/memorias/memorias2018/articulos/A5_40.pdf