Study of the influence of protuberances in the trailing edge of airfoils and determination of their aerodynamic efficiency through CFD using Ansys Fluent

==Resumen==

En la presente investigación se estudió el efecto de incorporar protuberancias en el borde de ataque al perfil aerodinámico NACA-0012. Para las simulaciones se diseñaron dos perfiles: perfil paramétrico (NACA 0012) y perfil modificado (0012-MOD), con ayuda del programa CAD Inventor Professional 2016. El perfil 0012-MOD incorpora modificaciones en el borde de ataque, cuya inspiración proviene de las malformaciones presentes en las aletas de la ballena jorobada (Megaptera novaeangliae). Las simulaciones computacionales se realizaron con el software comercial Ansys Fluent®, utilizando el modelo de viscosidad turbulenta k-ε realizable, con algoritmo de cálculo SIMPLE. El fluido de estudio es aire en el rango incompresible con Ma: 0,29. Con el objeto de obtener las curvas de sustentación y arrastre para ambos perfiles, se realizaron simulaciones para ángulos de ataque α entre los 0° y los 40°, de esta forma será posible comparar los perfiles y determinar su comportamiento aerodinámico.

'''Palabras clave''': Perfil aerodinámico, sustentación, CFD, protuberancias

==Abstract==

In the present research, the effect of incorporating protuberances in the leading edge to the aerofoil NACA-0012 was studied. For the simulations, two profiles were designed: parametric profile (NACA 0012) and modified profile (0012-MOD), with help of CAD Program Inventor Professional 2016. The 0012-MOD profile incorporates modifications in the leading edge, whose inspiration comes from malformations present on the flippers of the humpback whale (Megaptera novaeangliae). The computational simulations were performed in Ansys Fluent, using the turbulent viscosity model k-ε realizable, with SIMPLE calculation algorithm. The study fluid is air in the incompressible range with Ma: 0.11. To obtain the lift and drag curves for both profiles, they were simulated for angles of attack α: between 0° and 40°, this way it will be possible to compare the profiles and determine their aerodynamic behaviour.

'''Keywords''': Airfoil, lift, CFD, protuberances

==1. Introducción==

La ballena jorobada (Megaptera novaeangliae) ha sido nombrada como la más acrobática de las ballenas ya que puede realizar giros con ángulos de ataque superiores a los 40°. La alta maniobrabilidad de esta especie puede estar conectado directamente con su ecología de alimentación. Estas ballenas consumen grandes cantidades de plancton, arenque y capelín, y para capturarlos, la ballena necesita de complejas técnicas que requieren de habilidades de viraje, surfeo y buceo.

Un método ampliamente efectivo utilizado por las ballenas jorobadas se denomina alimentación por embestida y destaca por su estilo de nado a una velocidad aproximada de 2,6 m/s con un ángulo de ataque que va de los 30° a los 90° [1].

Otro método se conoce como burbujeo. La ballena genera una columna de burbujas al exhalar aire desde su orificio nasal mientras nada en ascenso. La red de burbujas resultantes encierra a su presa y la concentra de tal forma que la ballena puede consumir una gran cantidad en una maniobra de nado hacia arriba a través del centro de la red [2]. Un pequeño radio de giro implica que mayor cantidad de presas podrían ser concentradas en la red de burbujas, lo cual podría potencialmente incrementar la efectividad del método de alimentación. El radio de giro se puede determinar por la ecuación de fuerza inercial actuando en la ballena como función de la fuerza de sustentación generada por sus aletas [3].

Las protuberancias pueden ser vistas físicamente como generadores de vorticidad, permitiéndole a las ballenas transportar una mayor cantidad de energía entre las capas de fluido cercanas a sus aletas. Un mayor intercambio de energía permite desarrollar una capa límite turbulenta que es más resistente a flujo reverso si se compara con la capa límite laminar. Dicha resistencia permitirá retrasar la entrada en pérdida que experimentarían sus aletas sin la presencia de las protuberancias en el borde de ataque [4], por consiguiente, desarrolla mayores fuerzas aerodinámicas.

La [[#img-1|Figura 1]] permite esclarecer el concepto de protuberancias en las aletas de la ballena jorobada otorgando una idea del modelo conceptual que podría presentar un perfil aerodinámico al aplicar biomímesis en su proceso de diseño.

<div id='img-1'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Hinojosa_472673893-image1.jpeg|372px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 1'''. Protuberancias en aletas de ballena jorobada
|}


===1.1. Estado del arte===

Desde los años 70 en adelante que la eficiencia aerodinámica se a mejorado constantemente, principalmente piloteado por la industria aeronáutica. Desde hace algún tiempo se estudia la eficiencia en aplicaciones ligadas con las energías renovables: generadores eólicos, turbinas hidráulicas, etc. Nuevas tecnologías son necesarias, para llegar a nuevos límites en estas mejoras (perfiles modificables, control activo de flujo, materiales y estructuras más livianas, etc).

El estudio de las protuberancias en perfiles aerodinámicos no es algo nuevo, existen estudios numéricos y experimentales que analizan su influencia en el retraso de la separación de la capa límite para múltiples tipos de perfiles alares [5-15], en la reducción del ruido en aspas en generadores eólicos [16-18], en la mejora en la eficiencia en turbinas para aplicaciones de energía mareomotriz [19-21], entre otros.

En este trabajo en particular, analizamos los coeficientes de sustentación y arrastre para un perfil NACA-0012 con y sin protuberancias, obteniendo como resultado una mínima diferencia < 15% para estos parámetros. Sin embargo, se comprueba que las protuberancias retardan la aparición del flujo separado, permitiendo al perfil operar hasta valores del ángulo de ataque de 30°.

==2. Situación física==

El estudio se enfoca en el análisis tridimensional del comportamiento del flujo de aire sobre los perfiles aerodinámicos NACA 0012 Y 0012-MOD. Las propiedades del aire, así como el número de Reynolds (Re) y número de Mach (Ma) que desarrollan, se dan a conocer en la [[#tab-1|Tabla 1]].

<div class="center" style="font-size: 75%;">'''Tabla 1'''. Propiedades del fluido de estudio
</div>

<div id='tab-1'></div>
{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:55%;" 
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Tipo de material'''
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Fluido'''
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Aire'''
|-
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Propiedades'''
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Variabilidad'''
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Densidad (kg/m<sup>3</sup>)
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,225
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Constante
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Viscosidad (kg/ms)
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,789<math>\times</math>10<sup>-</sup><sup>5</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Constante
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Velocidad aire (m/s)
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|40
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Valor de ingreso
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Reynolds, Re
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|2,739<math>\times</math>10<sup>6</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|-
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|Número de Mach, Ma
|  colspan='2'  style="border: 1pt solid black;border-bottom: 1pt solid black;text-align: center;vertical-align: top;"|0,11
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|-
|}


Para dicho valor de Ma, se pueden despreciar los efectos de compresibilidad presentes en el flujo de aire, por lo cual se hace posible considerar al fluido de estudio como newtoniano [5]. Las consideraciones anteriores permiten establecer un grupo de ecuaciones gobernantes que rigen la situación física a ser simulada computacionalmente, mediante el uso del método de volúmenes finitos.

===2.1. Ecuaciones gobernantes===

Para un fluido incompresible, las ecuaciones gobernantes son la ecuación de continuidad y la ecuación de conservación de la cantidad de movimiento lineal [22-23]. La ecuación (1) representa la conservación de masa presente en un volumen de control, conocida en la literatura técnica como ecuación de continuidad.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math display="inline">\frac{\partial \rho }{\partial t}+\nabla (\rho \overset{\rightarrow}{u})=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}


Para evaluar la conservación de la cantidad de movimiento lineal, se deben considerar las ecuaciones de Navier-Stokes promediadas por Reynolds (RANS) [24].

La primera componente de las ecuaciones (2), (3) y (4) representa el gradiente de presiones que actúan en cada una de las direcciones principales; el segundo término asocia la viscosidad dinámica del fluido y el gradiente de velocidades al cuadrado, este producto se le denomina fuerza viscosa neta por unidad de volumen [25]. El último término se conoce con el nombre de esfuerzos de Reynolds.

En <math display="inline">\overset{\rightarrow}{x}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;"|<math>\frac{D\overset{\rightarrow}{u}}{Dt}=-\frac{\partial \overline{p}}{\partial x}+</math><math>\mu \left[ \frac{{\partial }^{2}\overline{u}}{\partial {x}^{2}}+\frac{{\partial }^{2}\overline{u}}{\partial {y}^{2}}+\right. </math><math>\left. \frac{{\partial }^{2}\overline{u}}{\partial {z}^{2}}\right] -\rho \left[ \frac{\partial \overline{{u}^{'}{u}^{'}}}{\partial x}+\right. </math><math>\left. \frac{\partial \overline{{u}^{'}{v}^{'}}}{\partial y}+\frac{\partial \overline{{u}^{'}{w}^{'}}}{\partial z}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

En <math display="inline">\overset{\rightarrow}{y}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
|  style="text-align: center;"|<math>\frac{D\overset{\rightarrow}{u}}{Dt}=-\frac{\partial \overline{p}}{\partial y}+</math><math>\mu \left[ \frac{{\partial }^{2}\overline{v}}{\partial {x}^{2}}+\frac{{\partial }^{2}\overline{v}}{\partial {y}^{2}}+\right. </math><math>\left. \frac{{\partial }^{2}\overline{v}}{\partial {z}^{2}}\right] -\rho \left[ \frac{\partial \overline{{v}^{'}{u}^{'}}}{\partial x}+\right. </math><math>\left. \frac{\partial \overline{{v}^{'}{v}^{'}}}{\partial y}+\frac{\partial \overline{{v}^{'}{w}^{'}}}{\partial z}\right]</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

En <math display="inline">\overset{\rightarrow}{z}</math>:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
|  style="text-align: center;"|<math>\frac{\partial \left( \rho k\right) }{\partial t}+\frac{\partial \left( \rho k{u}_{j}\right) }{\partial {x}_{j}}=</math><math>\left[ \left( \mu +\frac{{\mu }_{t}}{{\sigma }_{k}}\right) \frac{\partial k}{\partial {x}_{j}}\right] +</math><math>{G}_{k}+{G}_{b}-\rho \epsilon -{Y}_{M}+{S}_{k}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
|   style="text-align: center;"|<math>\frac{\partial \left( \rho k\right) }{\partial t}+\frac{\partial \left( \rho k{u}_{j}\right) }{\partial {x}_{j}}=</math><math>\left[ \left( \mu +\frac{{\mu }_{t}}{{\sigma }_{\epsilon }}\right) \frac{\partial \epsilon }{\partial {x}_{j}}\right] \rho {C}_{2}\frac{{\epsilon }^{2}}{k+\sqrt{\nu \epsilon }\, }+</math><math>{C}_{1\epsilon }\frac{\epsilon }{k}{C}_{3\epsilon }{G}_{b}+{S}_{\epsilon }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
|  style="text-align: center;"|<math>{\mu }_{t}=\rho {C}_{\mu }\frac{{k}^{2}}{\epsilon \, }</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}


El modelo de viscosidad turbulenta <math>k-\varepsilon</math> realizable incorpora dos nuevas ecuaciones conservativas al modelo. La primera de ellas evalúa la energía cinética turbulenta <math>k</math> (ecuación (5)), mientras que la ecuación (6) determina la tasa de disipación de la energía cinética turbulenta <math>\varepsilon</math> [26].

El valor de la viscosidad turbulenta <math>{\mu }_{t}</math> se determina de la combinación de <math>k</math> y <math>\varepsilon</math>, acorde a la ecuación (7).

Las constantes del modelo son:

{| style="width: 53%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="text-align: center;vertical-align: top;"|<math>{C}_{1\epsilon }=1,44</math>
|  style="text-align: center;vertical-align: top;"|<math>{C}_{2}=1,9</math>
|  style="text-align: center;vertical-align: top;"|<math>{\sigma }_{k}=1,0</math>
|  style="text-align: center;vertical-align: top;"|<math>{\sigma }_{\epsilon }=1,2</math>
|}

===2.2. Determinación del dominio de análisis===

Para determinar el dominio de análisis se consideraron recomendaciones de la NASA que establecen que el volumen de control debe ser lo suficientemente grande como para no influir en los resultados y a su vez optimizar los recursos computacionales de los cuales se dispone. El valor mínimo del dominio se estableció en 12,5 veces ''c''; donde ''c'' representa la cuerda del perfil. La [[#img-2|Figura 2]] muestra un esquema en términos generales de la forma que debe tener el dominio.

<div id='img-2'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Hinojosa_472673893-image2.png|294px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 2'''. Dominio recomendado para perfiles aerodinámicos
|}


===2.3. Diseño de las geometrías===

El diseño del perfil NACA 0012 se realizó utilizando el programa CAD Inventor® Professional 2016. Los puntos que permiten crear el perfil son generados por la herramienta en línea AirfoilTools®, que entrega un archivo en formato .csv el cual contiene 162 puntos coordenados que caracterizan al perfil NACA 0012 ([[#img-3|Figura 3]]), el cual fue extruido de 1,8 m.

<div id='img-3'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|[[Image:Draft_Hinojosa_472673893-image3-c.png|336px]] 
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 3'''. Perfil NACA 0012 en tres dimensiones
|}


Para el caso de la geometría modificada ([[#img-4|Figura 4]]) se requiere repetir el procedimiento realizado en el apartado anterior, variando el valor de la cuerda de forma de generar una serie de perfiles NACA 0012 con distintos valores de cuerda, que permitirán representar las protuberancias presentes en las ballenas jorobadas. Estas fueron modeladas con una ecuación sinusoidal de la forma <math>f(x)=80sin(x)</math> en mm.

<div id='img-4'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[Image:Draft_Hinojosa_472673893-image4.png|300px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 4'''. Perfil modificado de la forma <math>f(x)=A sin(x)</math>
|}


===2.4. Mallado del volumen a analizar===

Se realizaron tres refinamientos a la malla generada por defecto. El primero de ellos consiste en un refinado de cuerpo de tipo cuerpo de influencia. El tamaño de elementos definido es 30 cms.

El siguiente consiste en un refinamiento de cara. Se seleccionan las caras del perfil aerodinámico contenido en el dominio de análisis bajo la categoría de geometría, luego se define el tipo fijándolo en tamaño de elemento, estableciendo éste en 1 cm. Para evitar conflictos posteriores con la creación de la malla, es recomendable definir el comportamiento de estos elementos como firme.

Posteriormente, se realiza un nuevo refinado de cuerpo, esta vez seleccionando el tipo de refinamiento como esfera de influencia. Esta se define en función de un nuevo sistema de coordenadas y deben fijarse tanto el radio como el tamaño de elementos.

Una vez realizados los refinamientos, se consigue una malla de 12.556.440 celdas y calidad ortogonal 0,628675 para el perfil paramétrico. Para el caso del perfil 0012-MOD, la malla difiere tanto en calidad como en número de celdas, debido a que el módulo dedicado de Fluent™ para dichos efectos adapta la malla a la nueva geometría. La malla del perfil modificado se compone de 7.878.694 celdas y calidad ortogonal 0,456641.

Previo a la configuración del procedimiento de cálculo, se determinaron las condiciones de borde para cada zona del dominio ([[#tab-2|Tabla 2]]).

<div class="center" style="font-size: 75%;">'''Tabla 2'''. Condiciones de borde para las zonas del dominio
</div>
<div id='tab-2'></div>
{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;" 
|-
|  style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|'''Zona'''
|  style="border: 1pt solid black;border-left: 1pt solid black;border-right: 1pt solid black;text-align: center;"|'''Tipo'''
|  colspan='7'  style="border: 1pt solid black;border-left: 1pt solid black;text-align: center;"|'''Parámetros'''
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Airfoil
|  rowspan='2' style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|Pared
|  colspan='2'  style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|Movimiento de pared:
|  colspan='3'  style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|Estado de cizalladura:
|  colspan='2'  style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|Rugosidad de pared:
|-
|  colspan='2'  style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|''Estacionaria''
|  colspan='3'  style="bordert: 1pt solid black;border-right: 1pt solid black;text-align: center;"|''No deslizamiento''
|  colspan='2'  style="border: 1pt solid black;border-right: 1pt solid black;text-align: center;"|''constante''
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Far side 1
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Pared
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|Movimiento de pared
|  colspan='3'  style="border: 1pt solid black;text-align: center;"|Movimiento
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|Velocidad definida por:
|-
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|''En movimiento''
|  colspan='3'  style="border: 1pt solid black;text-align: center;"|''Relativo a zona de celdas adyacentes''
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|''Componentes''

''x, y, z''
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Far side 2
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Pared
|  colspan='2'  style="bordert: 1pt solid black;text-align: center;"|Movimiento de pared
|  colspan='3'  style="border: 1pt solid black;text-align: center;"|Movimiento
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|Velocidad definida por:
|-
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|''En movimiento''
|  colspan='3'  style="border: 1pt solid black;text-align: center;"|''Relativo a zona de celdas adyacentes''
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|''Componentes''

''x, y, z''
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Inlet
|  rowspan='2' style="border: 1pt solid black;text-align: center;"|Velocidad de entrada
|  colspan='5'  style="border: 1pt solid black;text-align: center;"|Método de especificación de turbulencia:
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|Velocidad definida por:
|-
|  style="border: 1pt solid black;text-align: center;"|''Intensidad''
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|
|  colspan='3'  style="border: 1pt solid black;text-align: center;"|''Relación de viscosidad''
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|''Componentes''

''x, y, z''
|-
|  style="border: 1pt solid black;text-align: center;"|Interior fluid domain
|  style="border: 1pt solid black;;text-align: center;"|Interior
|  colspan='7'  style="border: 1pt solid black;text-align: center;"|''Dominio interior del fluido de estudio''
|-
|  style="border: 1pt solid black;text-align: center;"|Outlet
|  style="border: 1pt solid black;text-align: center;"|Presión a la salida
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|Presión de inicio:
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|
|  colspan='3'  style="border: 1pt solid black;text-align: center;"|''0 (Pa) manométrico''
|-
|  style="border: 1pt solid black;text-align: center;"|Symmetry 1
|  style="border: 1pt solid black;text-align: center;"|Simetría
|  colspan='7'  style="border: 1pt solid black;text-align: center;"|''Cara de simetría''
|-
|  style="border: 1pt solid black;text-align: center;"|Symmetry 2
|  style="border: 1pt solid black;text-align: center;"|Simetría
|  colspan='7'  style="border: 1pt solid black;text-align: center;"|''Cara de simetría''
|}

===2.5. Configuraciones del procedimiento de cálculo===

Se seleccionó el modelo de viscosidad turbulenta k-ε realizable y fue considerada la función de pared estándar para el tratamiento de las celdas cercanas a la pared del perfil aerodinámico, recomendado en [25 y 27].

El algoritmo de cálculo adoptado corresponde al método semi-implícito para ecuaciones vinculadas a la presión (SIMPLE, por sus siglas en inglés), que es un procedimiento numérico que resuelve las ecuaciones de Navier-Stokes usando relaciones entre correcciones de velocidad y presión para forzar la conservación de masa y de esta forma obtener el campo de presiones y velocidades en todo el dominio.

Los recursos computacionales de los cuales se disponían para la investigación se detallan en la [[#tab-3|Tabla 3]].

<div class="center" style="font-size: 75%;">'''Tabla 3''': Recursos computacionales disponibles
</div>
<div id='tab-3'></div>
{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:50%;" 
|-
|  style="border: 1pt solid black;text-align: center;"|'''Procesador'''
|  style="border: 1pt solid black;text-align: center;"|'''Processor 6376'''
|-
|  style="border: 1pt solid black;text-align: center;"|Vel. procesador
|  style="border: 1pt solid black;text-align: center;"|2,30 GHz
|-
|  style="border: 1pt solid black;text-align: center;"|Núcleos físicos
|  style="border: 1pt solid black;text-align: center;"|14
|-
|  style="border: 1pt solid black;text-align: center;"|Memoria RAM
|  style="border: 1pt solid black;text-align: center;"|50 Gb
|}

==3. Resultados obtenidos==

===3.1. Conservación de masa y coeficientes aerodinámicos===

Una práctica recomendable para confirmar que el flujo de aire sigue la dirección presumida al inicio del análisis del problema, consiste en evaluar el flujo másico de aire en zonas de interés. Dadas las condiciones de borde establecidas para e el dominio, se consideraron cuatro zonas por las cuales podría escurrir aire. Se esperaba que todo el fluido que ingresa por la entrada (inlet), fuese posteriormente evacuado por la salida (outlet), mientras que para los lados lejanos 1 y 2 (far side 1 y 2) se monitoreó el comportamiento para así corroborar que todo el fluido sólo circulará por la entrada y la salida.

En todas las simulaciones no se observó entrada o salida de aire en los lados lejanos 1 y 2. Además, el flujo neto resultante es muy cercano a cero, más aún considerando el valor del residuo fijado en 10-3, lo que permite establecer que las dimensiones del perfil son adecuadas.

Los resultados obtenidos para el flujo másico neto se dan a conocer en la [[#tab-4|Tabla 4]].

<div class="center" style="font-size: 75%;">'''Tabla 4''': Conservación de masa y coeficientes aerodinámicos, perfil NACA 0012
</div>
<div id='tab-4'></div>
{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;font-size:85%;width:50%;" 
|-
|  rowspan='2' style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Ángulo de ataque α°'''
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Perfil NACA 0012'''
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Perfil 0012-MOD'''
|-
|  colspan='2'  style="border: 1pt solid black;text-align: center;vertical-align: top;"|'''Flujo neto (kg/s)'''
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|0°
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|5,01<math>\times</math>10<sup>-4</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,305<math>\times</math>10<sup>-5</sup>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|5°
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|-7,07<math>\times</math>10<sup>-6</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|3,11<math>\times</math>10<sup>-5</sup>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|10°
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|-1,24<math>\times</math>10<sup>-5</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,70<math>\times</math>10<sup>-5</sup>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|15°
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,28<math>\times</math>10<sup>-5</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|2,31<math>\times</math>10<sup>-5</sup>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|18°
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|-9,64<math>\times</math>10<sup>-3</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|2,25<math>\times</math>10<sup>-5</sup>
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|20°
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,02<math>\times</math>10<sup>-3</sup>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1,27<math>\times</math>10<sup>-5</sup>
|}

===3.2. Distribución del coeficiente de presión===

La [[#img-5|Figura 5]] muestra la distribución de presiones a lo largo de toda la superficie del perfil paramétrico para <math display="inline"> \alpha </math>: 18°, mostrando uniformidad en toda su envergadura. Destacan las zonas de presión negativa en la parte posterior cercana al borde de ataque y adicionalmente una región de alta presión bajo la línea que une el punto de estancamiento para un perfil perpendicular a la envergadura del perfil.
<div id='img-5'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"| [[File:Gonzalez_Hinojosa_2019a_7108_dist_pres_perfil_orig_18grados.png|414px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 5''': Distribución de presión a lo largo del perfil paramétrico, α: 18°
|}

De manera análoga, es posible observar en la [[#img-6|Figura 6]] la distribución de presiones a lo largo de toda la superficie del perfil modificado, considerando un α: 18°. En la superficie del perfil modificado se observa que la distribución de presión en el borde de ataque varía su magnitud debido a la incorporación de las protuberancias. Las zonas de baja presión se ven reducidas y concentradas en las concavidades de dichas modificaciones, mientras que para el resto del perfil se tiene una distribución más bien constante con valores del coeficiente de presión superiores al caso del perfil paramétrico.
<div id='img-6'></div>
{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|[[File:Gonzalez_Hinojosa_2019a_9995_dist_pres_perfil_modif_18grados.png|414px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 6''': Distribución de presión a lo largo del perfil modificado, α: 18°
|}

==4. Análisis de resultados==

La Figura 7 permite comparar las curvas de sustentación desarrolladas por los perfiles paramétrico (NACA 0012) y modificado (0012-MOD). Ambos perfiles engendran valores de sustentación similares para ángulos de ataque comprendidos entre 0° y 10°. Posterior al rango mencionado, se observa que el perfil NACA-0012 tiende a mantener constante el valor de sustentación desarrollado, etapa previa a la entrada en pérdida, que en la mayoría de los perfiles NACA ocurre para ángulos de ataque superiores a los 12°.

Por otra parte, el perfil 0012-MOD sigue incrementando el valor del coeficiente de sustentación con un comportamiento lineal, retrasando la entrada en pérdida aerodinámica para dicho perfil, tendencia que presenta un punto de inflexión para un ángulo de ataque de 30°, en donde éste experimenta la entrada en pérdida.

Las fuerzas de arrastre que ambos perfiles desarrollan presentan una evolución similar en función del ángulo de ataque (Figura 8). Como se esperaba, el coeficiente de arrastre para los dos perfiles sigue una tendencia cuadrática. Además, se observa que el perfil 0012-MOD tiende a desarrollan menores valores de arrastre que el NACA 0012, para ángulos de ataque comprendidos en el rango 5° ≤ α ≤ 15°.

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Hinojosa_472673893-image7.png|468px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 7''': Coeficiente de sustentación para perfiles paramétrico y modificado
|}


{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: auto;max-width: auto;"
|-
|style="padding:10px;"|  [[Image:Draft_Hinojosa_472673893-image8.png|468px]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" style="padding:10px;"| '''Figura 8''': Coeficiente de arrastre para perfiles paramétrico y modificado
|}

==5. Conclusiones==

En el trabajo investigativo realizado se estudiaron los efectos aerodinámicos debido a la incorporación de protuberancias en el borde de ataque del perfil aerodinámico NACA 0012.

Se representaron dos perfiles aerodinámicos con ayuda del programa de dibujo asistido Inventor Professional® 2016; el primero de ellos corresponde al perfil simétrico NACA 0012, denominado perfil paramétrico, el cual consistió en una extrusión del perfil bidimensional con envergadura de 1,8 m y cuerda de 1 m.

Para la construcción del segundo perfil, se utilizó el perfil paramétrico como base, sin embargo, se modificó el borde de ataque mediante el uso de la ecuación sinusoidal f(x)=Asin(Bx), donde A corresponde a la amplitud (80 mm) y B a la frecuencia de la función sinusoidal (en este caso, el valor es 1). La creación del perfil 3D se llevó a cabo mediante la opción de solevación.

Para ambos perfiles se realizaron simulaciones computacionales con ayuda del programa Ansys® Fluent, utilizando el modelo de turbulencia k-ε realizable. Se realizaron simulaciones a cada perfil, para ángulos de ataque entre 0 y 25° para el perfil paramétrico y 0° y 40° para el perfil modificado.

Para el perfil paramétrico el rango de estudio fue inferior ya que se presentaba una tendencia marcada del coeficiente de sustentación posterior a los 20° en la Figura 7, mientras que para el perfil modificado se analizó hasta encontrar un punto de inflexión característico de la curva de sustentación de los perfiles alares, que se observa para α: 30°.

Los resultados obtenidos de las simulaciones permitieron generar las curvas de sustentación y de arrastre para ambos perfiles, de esta forma fue posible analizar cuantitativamente el efecto que la incorporación de las protuberancias tendría sobre el perfil paramétrico.

Estos permiten establecer que el perfil modificado desarrolla hasta un 18,38% más de sustentación y 2,48% menos de arrastre que el perfil NACA 0012 para un ángulo de ataque de 18°.

Adicionalmente se observa una disminución de un 14,67% en el valor de las fuerzas de arrastre para el perfil modificado a un ángulo de ataque de 5°. Similares resultados se observan para α:10°, donde nuevamente el perfil modificado desarrolla menor arrastre, 13,58% inferior que el perfil paramétrico.

En términos generales, es posible establecer que el perfil modificado presenta mejoras significativas respecto del comportamiento aerodinámico observado para el perfil paramétrico en el rango de α > 12°, retrasando el fenómeno de pérdida aerodinámica.

===5.1. Limitaciones del estudio===

Si bien la calidad de la malla generada para ambos perfiles se encuentra por sobre las recomendaciones mínimas que aseguran la estabilidad en el proceso de cálculo y permiten confiar en los resultados obtenidos, se sugiere utilizar programas dedicados para el proceso de mallado.

En el proceso de búsqueda de soluciones para resolver las complicaciones experimentadas durante la etapa de diseño de la malla de los perfiles aerodinámicos, se encontraron programas dedicados, tales como ICEM CFD™, ANSA® y POINTWISE®. El uso de estos programas incrementa ostensiblemente la calidad ortogonal de la malla, lo que se traduce en resultados más confiables y disminución en los tiempos de cálculo requeridos para lograr convergencia del proceso iterativo. Sin embargo, es importante señalar que estos programas requieren de tal nivel de conocimientos y experiencia en su uso, que se requiere realizar cursos de formación para no incurrir en errores que se puedan propagar a los resultados de las simulaciones.

===5.2. Proyección de la investigación===

Debido a que sólo se realizaron simulaciones con una única variante del perfil NACA 0012, se deja manifiesta la necesidad de continuar con el proceso investigativo, variando los parámetros de la función sinusoidal, lo que permitirá crear más perfiles modificados y de esta forma contar con un universo más amplio de perfiles, lo que permitirá determinar la configuración final del perfil modificado.

Adicionalmente, se recomienda realizar mediciones experimentales, para lo cual se hace necesaria la construcción de prototipos de los perfiles modificados, los cuales deben ser estudiados en túneles de viento, lo que permitirá la creación de sus respectivas curvas aerodinámicas. Contrastar los datos experimentales y los obtenidos de las simulaciones computacionales permitirá validar científicamente las conclusiones de la presente investigación.

Seguido de esto, se recomienda considerar potenciales aplicaciones que tendría la incorporación de dichas protuberancias a perfiles aerodinámicos. En forma preliminar es posible pensar en aplicar tales como: nuevo diseño de palas de generadores eólicos; incorporación de modificaciones en las aspas del rotor principal de un helicóptero; diseño de las aspas de ventiladores y extractores; álabes de turbinas hidráulicas; alas de aviones; entre otros, todos en el rango subsónico.

==Referencias==

<div class="auto" style="text-align: left;width: auto; margin-left: auto; margin-right: auto;font-size: 85%;">

[1] Weinrich M., Mann J., Connor R., Tyack P., Whitehead H. Cetacean societies: field studies of dolphines and whales. Journal of Wildlife Management, 65(2):366-367, 2001.

[2] Acevedo J., Plana J., Aguayo A., Pastene L. Surface feeding behavior of humpback whales in the Magellan Strait. Rev. Biol. Mar. Oceanogr., 46(3):483-490, 2011.

[3] Weihs D. A hydrodynamical analysis of fish turning manoeuvres. Proceedings of the Royal Society of London. Series B: Biological Sciences, 182, 59-72, 1972.

[4] Van Nierop E.A., Alben S., Brenner M. How bumps on whale flippers delay stall: An aerodynamic model. Physical Review Letters, 100(054502), 2008.

[5] Wang Y., Hu W., Zhang S. Performance of the bio-inspired leading edge protuberances on a static wing and pitching wing. Journal of Hydrodynamics, 26(6):912-920, 2014.

[6] Lin S.-Y., Lin Y.-Y., Bai C.-J., Wang W.-C. Performance analysis of vertical-axis-wing-turbine blade with modified trailing edge through computational fluid dynamics. Renewable Energy, 99:654-662, 2016.

[7] Wang Z., Zhuang M. Leading edge serrations for performance improvement on a vertical-axis wind turbine al low tip-speed-ratios. Applied Energy, 208:1184-1197, 2017.

[8] Favier J., Pinelli A., Piomelli U. Control of the separated flow around an airfoil using a wavy leading edge inspired by humpback whale flippers. Comptes Rendus Mecanique, 340:107-114, 2012.

[9] Bolzon M.D.P., Kelso R.M., Arjomandi M. Formation of vortices on a turbercled wing, and their effects on drag. Aerospace Science and Technology, 56:46-55, 2016.

[10] Wei Z., New T.H., Cui Y.D. An experimental study on flow separation control of hydrofoils with leadiing-edge tubercles at low Reynolds number. Ocean Engineering, 108:336-349, 2015.

[11] Serson D., Meneghini J.R. Numerical study of wings with wavy leading and trailing edges. Procedia IUTAM, 14:563-569, 2015.

[12] Yoon H.S., Hung P.A., Jung J.H., Kim M.C. Effect of the wavy leading edge on hydrodynamic characteristics for flow around low aspect ratio wing. Computers & Fluids, 4(1):276-289, 2011.

[13] Aftab S.M.A., Razak N.A., Mohd Rafie A.S., Ahmad K.A. Mimicking the humpback whale: An aerodynamic perspective. Progress in Aerospace Sciences, 84:48-69, 2016.

[14] Sudhakar S., Karthikeyan N., Venkatakrishnan  L. Influence of leading edge tubercles on aerodynamic characteristics of a high aspect-ratio UAV. Aerospace Science and Technology, 69:281-289, 2017.

[15] Huang C.-C., Bai C.-J., Shiah Y.C., Chen Y.-J. Optimal design of protuberant blades for small variable-speed horizontal axis wind turbine-experiments and simulations. Energy, 115:1156-1167, 2016.

[16] Lau A.S. H., Haeri S., Kim J.W. The effect of wavy leading edge on aerofil-gust interaction noise. Journal of Sound and Vibration, 332:6234-6253, 2013.

[17] Tong F., Qiao W., Xu K., Wang L., Chen W., Wang X. On the study of wavy leading-edge vanes to achieve low fan interaction noise. Journal of Sound and Vibration, 419:200-226, 2018.

[18] Chen W., Qiao W., Tong F., Wang  L., Wang X. Experimental investigation of wavy leading edges on rod-aerofoil interaction noise. Journal of Sound and Vibration, 422:409-431, 2018.

[19] Shi W., Rosli R., Atlar M., Norman R., Wang D., Yang W. Hydrodynamic performance evaluation of a tidal turbine with leading-adge tubercles. Ocean Engineering, 117:246-253, 2016.

[20] Shi W., Atlar M., Norman R., Aktas B., Turkmen S. Numerical optimization and experimental validation for a tidal turbine blade with leading-edge tubercles. Renewable Energy, 96(A):42-55, 2016.

[21] Shi W., Atlar M., Norman R. Detaied flow measurement of the field around tidal turbines with and without biomimetis leading-edge tubercles. Renewable Energy, 111:688-707, 2017.

[22] Anderson J. Fundamentals of Aerodynamics, 6th Edition, McGraw-Hill, USA, 2008.

[23] White F. Fluid Mechanics, McGraw-Hill, USA, 2011. 

[24] Kundu C., Cohen I. Fluid mechanics. 2nd Edition, Elsevier, USA, 2002. 

[25] Launder B., Spalding D. The numerical computation of turbulent flows. Computer Methods in Applied Mechanics and Engineering, 3(2):269-289, 1974.  

[26] Cengel Y. Mecánica de fluidos. McGraw-Hill, 2006.

[27] Versteeg H.K., Malalasekera W. An Introduction to computational fluid dynamics: The Finite Volume Method.  Pearson Educationan Limited, 2007.
</div>