VECTORIZACIÓN DE LINEAS

CORRIENTE EN PERFILES
AERODINÁMICOS MEDIANTE EL
USO DEL METODO VLM EN
PYTHON
ADRIÁN BONO BAREA
MASTER DE MATEMÁTICAS
UNA BREVE INTRODUCCIÓN A LA
HISTORIA DE LA AERODINAMICA
LA AERODINAMICA COMO CAMPO DE ESTUDIO
El estudio aerodinámico se engloba dentro de la mecánica de fluidos como una rama más de la física que
pretende dilucidar el comportamiento de los fluidos cuando estos circulan en las proximidades de un
cuerpo.

Los primeros pasos realizados en este campo dieron a conocer la enorme complejidad de la que se reviste.
Hoy en día, gracias a los avances tecnológicos en las ciencias de la computación que nos permiten tanto
el análisis como la evaluación de modelos discretizando dominios, se abre un amplio abanico de
posibilidades hacia la realización de nuevos diseños.

A lo largo de este apartado se darán a conocer las diferentes soluciones que se han dado y sus aplicaciones
a la aviación.

1738 – LA ECUACIÓN DE BERNOULLI

Aunque el comportamiento de los fluidos se remonta a la época de Arquímedes, quien fue el primero en
describir el flujo como causa de un gradiente de presiones, los orígenes de la aerodinámica no se dan
hasta el siglo XVIII con el matemático Daniel Bernoulli.

Bernoulli describió por primera vez en su obra “Hydrodynamica” la relación existente entre presión y
velocidad hoy conocida como principio de Bernoulli. Este enunciado nos dice que la presión disminuye en
proporción que la velocidad de un determinado fluido aumenta.

𝑣12 𝑝1 𝑣22 𝑝2
+ = + (1)
2 ρ 2 ρ

Esta relación se incluye como una solución particular a un flujo incompresible, adiabático y no viscoso.
Bernoulli no pudo extender su solución fuera de este dominio.

2
1757 – LAS ECUACIÓNES DE EULER
Unos años más tarde de los enunciados de Bernoulli, Leonhard haría público su trabajo "Principes
généraux du mouvement des fluides" que extendía la solución al flujo compresible. Las ecuaciones están
descritas como un set de 3 ecuaciones quasilineares hiperbólicas que describen las variaciones que sufren
la velocidad del flujo, la densidad, la presión y la energía como propiedades del fluido.

𝜕𝜌
+ 𝑢 ∙ ∇𝜌 + 𝜌∇ ∙ 𝑢 = 0
𝜕𝑡
𝜕𝑢 ∇𝑝
+ 𝑢 ∙ ∇𝑢 + =𝑔 (2)
𝜕𝑡 𝜌
𝜕𝑒 𝑝
{ 𝜕𝑡 + 𝑢 ∙ ∇𝑒 + ∇ ∙ 𝑢 = 0
𝜌

Se puede demostrar que para el conjunto de ecuaciones se puede llegar al principio de Bernoulli y al igual
que este, Euler consiguió dar una solución particular donde no se tienen en cuenta los fenómenos viscosos.

1822 Y 1843 – LAS ECUACIÓNES DE NAVIER-STOKES

Casi un siglo después, el físico Claude-Louis Navier y el matemático George Gabriel Stokes consiguen
describir el comportamiento de los fluidos viscosos mediante la aplicación de la segunda ley de newton a
un fluido. Además, asumen que el estrés que se producido en el propio es el resultado de la suma de los
fenómenos viscosos que son proporcionales a su velocidad.

𝐷𝑢 𝜕𝑢 1
𝜌 = 𝜌( + 𝑢 ∙ ∇𝑢) = −∇𝑝 + 𝜇∇2 𝑢 + 𝜇∇(∇ ∙ 𝑢) + 𝜌𝑔 (3)
𝐷𝑡 𝜕𝑡 3

Las ecuaciones de Navier-Stokes logran describir el comportamiento de muchos fenómenos físicos, desde
modelos meteorológicos y oceánicos, hasta las fuerzas que se producen cuando un fluido circula a través
de un objeto. Sus aplicaciones alcanzan también la medicina con el estudio de modelos de diseminación
de medicamentos a través del torrente sanguíneo. Además, acopladas a las ecuaciones de Maxwell
pueden ser usadas en el estudio de la magneto-hidrodinámica.

No solo atraen la atención de los físicos e ingenieros, también de los matemáticos. A pesar de que estas
ecuaciones describen multitud de fenómenos, no se ha conseguido probar que existe siempre una
solución para cualquier problema en tres dimensiones, y si existe solución, es fluida. El “Clay Mathematics
Institute” lo considera como uno de los últimos problemas no resueltos de la física clásica.

3
1960 HASTA LA ACTUALIDAD – CFD
Las ecuaciones de Navier-Stokes sirven como fundamento a muchos modelos de la mecánica de fluidos
computacional.

Esta rama hace uso del enorme poder de cómputo que ofrecen los ordenadores para resolver el
comportamiento de fluidos en circunstancias determinadas. Uno de los pioneros de este campo es Francis
H. Harlow, quien, junto con su equipo, desarrollo multitud de algoritmos para el cálculo de variables sobre
flujos transitorios en 2-D.

En la actualidad existen modelos que describen con gran precisión flujos tridimensionales incluso en
régimen turbulento. Algunos más complejos pueden ser el método “Direct Númerical Simulation” o DNS
y el método “Reynold-Averaged Navier-Stokes” o RANS que son utilizados para la modelización de flujos
de chorro en cohetes y toberas.

“VORTEX LATTICE METHOD” o VLM

CONTRUCCIÓN DEL METODO A TRAVES DE NAVIER-STOKES
En este documento se construirá uno de los métodos más usados en la mecánica de fluidos computacional
dada su sencillez y su fácil computo, este es el “método de paneles” o “Vortex Lattice Method”. Además,
se hará un estudio comparativo del modelo simulando, mediante el uso de Python, el comportamiento
de distintos perfiles aerodinámicos y enfrentándolo a la teoría linealizada de perfiles.

El método de paneles nos permite evaluar las ecuaciones de Navier-Stokes usando la propiedad de la
superposición de soluciones elementales de tal modo que la suma de las mismas dé como resultado la
solución del campo fluido que deseemos hallar. Para ello se deben tener las siguientes consideraciones
partiendo de las ecuaciones del momento y de continuidad de Navier-Stokes, como la condición de flujo
no viscoso. Este hecho elimina los términos del tensor de fuerzas viscosas de la ecuación (3) quedando:

𝐷𝑢 𝜕𝑢 1
𝜌 = 𝜌( + 𝑢 ∙ ∇𝑢) = −∇𝑝 + 𝜇∇2 𝑢 + 𝜇∇(∇ ∙ 𝑢) + 𝜌𝑔 (4)
𝐷𝑡 𝜕𝑡 3

𝐷𝜌
= 𝜌∇ ∙ 𝑢 = 0 (5)
𝐷𝑡

4
Otra suposición que se hace es la de flujo estacionario con fuerzas másicas despreciables, reduciendo las
expresiones a:

𝜕𝑢 1
𝜌( + 𝑢 ∙ ∇𝑢) = −∇𝑝 + 𝜇∇2 𝑢 + 𝜇∇(∇ ∙ 𝑢) + 𝜌𝑔 (6)
𝜕𝑡 3

𝐷𝜌
= 𝜌∇ ∙ 𝑢 = 0 (7)
𝐷𝑡

Se pueden aprovechar las siguientes propiedades de los operadores diferenciales para expresar la
ecuación de la forma (9):

|𝑎2 |
𝑎 ∙ ∇𝑎 = (∇ × 𝑎) × 𝑎 + ∇ ( ) (8)
2

|𝑢2 | 𝑝
(∇ × 𝑢) × 𝑢 + ∇ ( ) = −∇ 𝜌 (9)
2

𝜌∇ ∙ 𝑢 = 0 (10)

El primer término de la expresión (9) tiene relación con el rotacional de un flujo, en un flujo sin vorticidad
se da:

∇×𝑢 =0 (11)

Por lo tanto, la expresión queda reducida de la siguiente forma:

|𝑢2 | 𝑝
∇( ) = −∇ 𝜌 (12)
2

𝜌∇ ∙ 𝑢 = 0 (13)

5
Podemos definir una función potencial del campo de velocidades de tal modo que:

∇φ = 𝑢 (14)

Sustituyendo en las expresiones:

|∇φ2 | 𝑝
∇( ) = −∇ 𝜌 (15)
2

∇2 𝜑 = 0 (16)

De la ecuación de la continuidad se ha llegado a la expresión (16) que no es más que la ecuación de

Laplace. Esta ecuación cumple la propiedad de la superposición que hemos mencionado antes.

APLICACIÓN DE LA SUPERPOSICIÓN DE SOLUCIÓNES A UN PERFIL

El objetivo de la superposición es el de proporcionar una función potencial que satisfaga las condiciones
del flujo sujeto a estudio. El método de paneles discretiza la geometría de un cuerpo en paneles al que
dispone torbellinos sobre sus vértices de la siguiente forma:

Figura 1: Disposición de potencial de torbellinos sobre una curva discretizada en paneles.

6
De este modo podemos evaluar la velocidad del flujo en un punto P(x,z) influida por los torbellinos γ . Si
llevamos dicho punto de control sobre la superficie del panel, obtenemos la siguiente expresión en
coordenadas polares dada por la suma de todos los torbellinos potenciales:

γ
𝜑𝑡 = ∑𝑛_𝑣𝑒𝑟𝑡𝑖𝑐𝑒𝑠
𝑖=1
𝑖
∑𝑛_𝑣𝑒𝑟𝑡𝑖𝑐𝑒𝑠
𝑗=1 𝜃𝑖𝑗 (17)
2𝜋

El potencial obtenido tiene a priori infinitas soluciones por lo que deberemos escoger aquellas que casen
con el modelo físico. De otro modo obtendríamos una solución que, a pesar de cumplir la expresión
anterior, sería de poca utilidad para su posterior análisis. En relación con lo anterior se imponen las
llamadas “condiciones de contorno”, una de ellas es la condición de no penetrabilidad y es que el flujo no
puede atravesar las paredes de un objeto sólido. Para ello se establece que la velocidad en los puntos de
control sea nula o lo que es lo mismo que el potencial de los torbellinos sobre el punto de control se anule.
Esto es:

𝜑𝑡 + 𝜑 𝑓 = 0 (18)

Desarrollando el sistema en forma matricial:

𝑓
𝑎11 𝑎12 ⋯ 𝑎1𝑁 𝑎1,𝑁+1 γ1 𝜑1
𝑎31 𝑎22 ⋯ 𝑎2𝑁 𝑎2,𝑁+1 γ2 𝑓
𝜑2
𝑎31 𝑎32 ⋯ 𝑎3𝑁 𝑎3,𝑁+1 γ3 𝑓
= 𝜑3 (19)
⋮ ⋮ ⋱ ⋮ ⋮ ⋮
γ𝑁 ⋮
𝑎𝑁1 𝑎𝑁2 𝑎𝑁3 𝑎𝑁𝑁 𝑎𝑁,𝑁+1 𝑓
( 1 ) ( γ 𝜑𝑁
0 0 0 1 𝑁+1 )
(0)
Siendo aij:

𝜃𝑖𝑗
𝑎𝑖𝑗 = (20)
2𝜋

En el sistema matricial se ha incluido la condición de Kutta que impone que la circulación se iguale en el
borde de salida, con ello obtendremos una circulación neta no nula que será la responsable de la aparición
de fuerzas sobre el perfil. Este “pequeño” arreglo permite aproximar la solución matemática a la realidad
dotándola de un sentido físico ya que cualquier cuerpo inmerso en un fluido experimenta una fuerza por
el mismo.

7
Aplicando la resolución del sistema matricial obtenemos la circulación de cada torbellino. Por la definición
de potencial, obtenemos la velocidad de todo el dominio fluido aplicando su la derivada. En el sentido
discreto del problema se debe recurrir a la definición de la derivada mediante el concepto de T.V.M en un
paso infinitesimal:

𝜕𝑓(𝑥) 𝑓(𝑥+ℎ)−𝑓(𝑥−ℎ)
~ (21)
𝜕𝑥 2ℎ

Una vez obtenido el campo de velocidades podemos calcular el campo de coeficiente de presiones
mediante la ecuación de Bernoulli, esta es:

𝜕𝜑𝑖 2
𝐶𝑝𝑖 = 1 − ( ) (22)
𝜕𝑥

Del campo de coeficiente de presiones se obtiene el coeficiente de sustentación y el de resistencia.

𝑝 𝑛 .𝑐𝑜𝑛𝑡𝑟𝑜𝑙
𝐶𝑙 = ∑𝑖=1 Cp𝑖 ∙ sen(𝛼) (23)

𝑝 𝑛 .𝑐𝑜𝑛𝑡𝑟𝑜𝑙
𝐶𝑑 = ∑𝑖=1 Cp𝑖 ∙ cos(𝛼) (23)

DESARROLLO DEL METODO EN PYTHON

Teniendo en cuenta el trabajo previamente realizado en el apartado anterior, se procede a transcribir este
desarrollo en Python.

Para comenzar, se imponen las condiciones iniciales de nuestro problema y se importa la geometría del
perfil que estará definida por puntos con coordenadas x,y.

8
 Figura 2: Carga de librerías, imposición de condiciones del problema e importación de la geometría.

A continuación, se adapta la geometría dada por el archivo de puntos de manera que se obtenga una
geometría cerrada si esta no lo es. Para ello, en el array p_perfil, se concatena con un punto más, que será
coincidente con el punto inicial o final cerrando la geometría del perfil.

En una operación posterior se procederá a recorrer la curva del perfil con el fin de obtener su tamaño
para poder ser luego dividida entre el número de nodos establecidos previamente, con ello se pretende
discretizar toda la superficie del perfil en paneles del mismo tamaño.

Figura 3: Geometría de puntos equiespaciados.

9
 Figura 4: Suma concatenada de la distancia entre puntos del array p_perfil e interpolación de la curva en puntos
equiespaciados.

Para finalizar las características geométricas de nuestro perfil se construyen los paneles, así como los
puntos de control situados a la mitad del panel y los vectores normales y tangenciales al panel.

Figura 5: Declaración de las características de los paneles

El siguiente paso será definir las matrices de nuestro sistema matricial. Para empezar, se construye la
matriz A recorriendo cada punto de control y calculando el coeficiente ai,j que dependerá de la situación
del vórtice sobre el que se esté calculando al punto de control.

Figura 6: Construcción de las matrices A,b e imposición de la condición de Kutta.

10
Por último, para la resolución del sistema matricial se hará uso de la librería scipy, en concreto de las
funciones LU_solver y LU_factor.

Figura 7: Resolución del sistema matricial y cálculo de la circulación sobre cada torbellino

RESULTADOS Y ANALISIS COMPARATIVO CON XFLR-5

Una vez halladas las circulaciones pueden calcularse multitud de parámetros físicos relacionados con la
influencia que produce el flujo sobre el perfil. Para este estudio se ha pretendido seleccionar aquellos de
mayor relevancia como el coeficiente de sustentación, resistencia y las líneas de corriente.

Para el cálculo de las líneas de corriente es necesario entender que son realmente, estrictamente se
definen como aquellas donde la velocidad es tangente en todo punto del dominio fluido. Dicho de una
manera mas sencilla, son el camino que recorre una partícula cuando esta atraviesa las inmediaciones de
un objeto.

Un modo de representar dicha solución es imponer un array de puntos situados aguas arriba del perfil; la
velocidad de estos vendrá dada por la evaluación en cada punto.

Figura 8: Creación del array de partículas

Para evaluar esta velocidad se ha creado una función que dependerá exclusivamente de dos variables x,y.
Dicha función dará como resultado el potencial calculado sobre el punto. Con el fin de aumentar la rapidez
con la que se realiza estos cálculos se ha acudido a la librería numba que se ayuda de la optimización de
operaciones mediante el uso de código C, de ejecución más rápida y el multithreading.

11
Figura 9: Definición de la función potencial, uso de la función JIT de numba y evaluación de la velocidad mediante la derivación
del potencial por diferencias finitas.

Para finalizar el primer análisis, se muestran las líneas de corriente en la siguiente figura:

Figura 10: Representación de las líneas de corriente.

12
Como se puede ver en la imagen 10, las líneas de corriente se adaptan a la superficie del perfil rodeándolo,
además se produce un leve aumento de la velocidad en el extradós del perfil permitiendo que las
partículas recorran mayor distancia. Este hecho casa perfectamente con el fenómeno físico que se
producen en los aviones y es que; al aumentar la velocidad en el extradós disminuye la presión haciendo
que la componente neta del campo de presiones se traduzca en una fuerza de componente mayoritario
vertical denominada sustentación.

El calculo de esta fuerza junto con la resistencia serán objeto de estudio para el siguiente análisis en el
que se graficará la influencia que tiene el ángulo de ataque en dichas fuerzas. Para ello se evaluará el
resultado para cada ángulo en un rango.

Figura 11: Evaluación de coeficientes para un rango de ángulos dado.

Para finalizar se graficará el resultado:

13
 Figura 12: Gráfica de Cl y Cd frente a Alpha

Como punto final a nuestro trabajo debemos estudiar la validez de nuestros resultados comparándolo con
un estudio in situ o mediante el uso de alguna herramienta CFD. En este caso nos ayudaremos del software
XFLR5 desarrollado por el MIT.

Figura 13: Interfaz de creación de perfiles.

14
Realizaremos un análisis simple teniendo en cuenta un numero de Reynolds bastante bajo, entorno a los
50.000, esta simulación se hará en un rango de valores del ángulo de ataque entre -6º y 10º, los mismos
que hemos usado en nuestro código.

Figura 13: Simulación realizada.

Finalizado el análisis, extraemos los siguientes resultados:

α= 6.63º Python (VLM) XFLR5 (LLT)

Cl 1,11 1,08
Cd 0,0598 0,0372

Como se puede observar en la tabla, para un ángulo de ataque determinado, en este caso el ángulo de
ataque ideal con valor de 6,63º se obtienen pocas diferencias con respecto a otros métodos si hablamos
en términos de sustentación. Sin embargo, sí que se produce una ligera variación entre el método
desarrollado para Python (VLM) y el método usado por XFLR5(LLT), esto puede deberse tanto a diferencias
en la discretización o a las suposiciones que hace cada metodología para realizar el estudio.

Sin duda, el código desarrollado además de ser simple es funcional y puede ser de una gran utilidad para
el cálculo de fuerzas aerodinámicas en regímenes de bajas velocidades.

15