Aunque no te lo creas, con lo que hemos visto hasta ahora eres capaz de hacer grandes cosas. Vale sí, un perfil de Yukovski no es gran cosa aerodinámicamente, pero si lo hacemos en Python... Echa un vistazo a la figura ¿no está mal, no? algo así intentaremos conseguir al final de esta clase.
Como no se trata de aprender (o reaprender) aerodinámica, te daremos las funciones matemáticas y los pasos a seguir así como la estructura del programa. Tú sólo tienes que preocuparte de programar cada bloque. Puedes leer en detalle todo lo relativo a la aerodinámica en el libro Aerodinámica básica de Meseguer Ruiz, J., Sanz Andrés, A. (Editorial Garceta).
Lo primero es lo primero, importemos los paquetes:
# Recuerda, utilizaremos arrays y pintaremos gráficas.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
La transformación de Yukovski es: $$\tau=t+\frac{a^{2}}{t}$$
Los parámetros del problema son los del siguiente bloque, puedes cambiarlos más adelante:
# Datos para el perfil de Yukovski
# Parámetro de la transformación de Yukovski
a = 1
# Centro de la circunferencia
landa = 0.2 # coordenada x (en valor absoluto)
delta = 0.3 # coordenada y
t0 = a * (-landa + delta * 1j) # centro: plano complejo
# Valor del radio de la circunferencia
R = a * np.sqrt((1 + landa)**2 + delta**2)
# Ángulo de ataque corriente incidente
alfa_grados = 0
alfa = np.deg2rad(alfa_grados)
#Velocidad de la corriente incidente
U = 1
Se trata de definir una función que realice la transformación de Yukovski. Esta función recibirá el parámetro de la transformación, $a$ y el punto del plano complejo $t$. Devolverá el valor $\tau$, punto del plano complejo en el que se transforma $t$.
def transf_yukovski(a, t):
"""Dado el punto t (complejo) y el parámetro a
a de la transformación proporciona el punto
tau (complejo) en el que se transforma t."""
tau = t + a ** 2 / t
return tau
#comprobamos que la función está bien programada
#puntos del eje real siguen siendo del eje real
err_message = "La transformación de Yukovski no devuelve un resultado correcto"
np.testing.assert_equal(transf_yukovski(1, 1+0j), 2+0j, err_message)
Ahora queremos transformar la circunferencia de radio $R$ con centro en $t_0$ usando la función anterior:
N
puntos de la circunferencia de modo que en Xc
estén las coordenadas $x$ y en Yc
estén las coordenadas $y$ de los puntos que la forman. Controla el número de puntos mediante un parámetro que se llame N_perfil
.
$$X_c = real(t_0) + R·cos(\theta)$$
$$Y_c = imag(t_0) + R·sin(\theta)$$Xc
e Yc
, píntalos mediante un scatter
para comprobar que todo ha ido bien.Deberías obtener algo así:
# Número de puntos de la circunferencia que
# vamos a transformar para obtener el perfil
N_perfil = 100
#se barre un ángulo de 0 a 2 pi
theta = np.linspace(0, 2*np.pi, N_perfil)
#se crean las coordenadas del los puntos
#de la circunferencia
Xc = np.real(t0) + R * np.cos(theta)
Yc = np.imag(t0) + R * np.sin(theta)
#lo visulaizamos
plt.figure("circunferencia", figsize=(5,5))
plt.title('Circunferencia', {'fontsize':20})
plt.scatter(Xc, Yc)
plt.scatter(np.real(t0), np.imag(t0), color='orange', marker='x', s=100)
plt.grid()
# Lo visulaizamos más bonito
plt.figure("circunferencia", figsize=(5,5))
plt.title('Circunferencia', {'fontsize':20})
# Esto no tienes por qué entenderlo ahora
p = plt.Polygon(list(zip(Xc, Yc)), color="#cccccc", zorder=3)
plt.gca().add_patch(p)
plt.ylim(-1.5, 2)
plt.xlim(-2, 1.5)
plt.grid()
Ahora estamos en condiciones de transformar estos puntos de la circunferencia (Xc
, Yc
) en los del perfil (Xp
, Yp
). Para esto vamos a usar nuestra función transf_yukovski
. Recuerda que esta función recibe y da números complejos. ¿Saldrá un perfil?
Puntos_perfil = transf_yukovski(a, Xc+Yc*1j)
Xp, Yp = np.real(Puntos_perfil) , np.imag(Puntos_perfil)
#lo visulaizamos
plt.figure("perfil yukovski", figsize=(10,10))
plt.title('Perfil', {'fontsize':20})
plt.scatter(Xp, Yp)
plt.grid()
plt.gca().set_aspect(1)
#lo visulaizamos más bonito
plt.figure('perfil yukovski', figsize=(10,10))
plt.title('Perfil', {'fontsize':20})
p = plt.Polygon(list(zip(Xp, Yp)), color="#cccccc", zorder=3)
plt.gca().add_patch(p)
plt.gca().set_aspect(1)
plt.xlim(-3, 3)
plt.ylim(-0.4,1)
plt.grid()
Para visualizar ahora el flujo alrededor del cilindro recurrimos al potencial complejo de una corriente uniforme que forme un ángulo $\alpha$ con el eje $x$ en presencia de un cilindro (aplicando el teorema del círculo) y se añade un torbellino con la intensidad adecuada para que se cumpla la hipótesis de Kutta en el perfil:
\begin{equation} f(t)=U_{\infty}\text{·}\left((t-t_{0})\text{·}e^{-i\alpha}+\frac{R^{2}}{t-t_{0}}\text{·}e^{i\alpha}\right)+\frac{i\Gamma}{2\pi}\text{·}ln(t-t_{0})=\Phi+i\Psi \end{equation}donde $\Phi$ es el potencial de velocidades y $\Psi$ es la función de corriente.
$$\Gamma = 4 \pi a U (\delta + (1+\lambda) \alpha)$$$\Gamma$ es la circulación que hay que añadir al cilindro para que al transformarlo en el perfil se cumpla la condición de Kutta.
Recordando que la función de corriente toma un valor constante en las líneas de corriente, sabemos que: dibujando $\Psi=cte$ se puede visualizar el flujo.
Pintaremos estas lineas de potencial constante utilizando la función contour()
, pero antes tendremos que crear una malla con meshgrid()
. Esto será lo primero que hagamos:
Crea un parámetro N_malla
cuyo valor sea el número de puntos que va a tener la malla en cada dirección ($x$ e $y$).
Crea dos parámetros x_max
e y_max
que representen los valores extremos de la malla en $x$ e $y$ respectivamente.
Crea un array x
que vaya desde -x_max
hasta x_max
y que tenga N_malla
elementos. De manera análoga, crea el array y
.
Crea la malla XX, YY
utilizando la función meshgrid()
que necesitaremos para pintar con contour()
. Los arrays XX
e YY
son de la forma:
-x_{max} & & ... & & & +x_{max} \ & & & & & \ ... & & & & & \ & & & & & \ & & & & & \ -x_{max} & & ... & & & +x_{max} \ \end{pmatrix}
... & & & & & \ & & & & & \ & & & & & \ +y_{max} & & ... & & & +y_{max} \ \end{pmatrix}
#se crea la malla donde se va pintar la función de corriente
N_malla = 50 #nº de puntos en cada dimensión
x_max = 8
y_max = 3
x = np.linspace(-x_max, x_max, N_malla)
y = np.linspace(-y_max, y_max, N_malla)
XX, YY = np.meshgrid(x,y)
#pintamos la malla para verla
plt.figure(figsize=(10,10))
plt.scatter(XX.flatten(), YY.flatten(), marker='.')
<matplotlib.collections.PathCollection at 0x7f7c6d783208>
Crea una variable T
que tenga el valor correspondiente a la circulación $\Gamma$.
Utilizando XX
e YY
, crea un array tt
que exprese la malla de forma compleja:
\begin{pmatrix} -x_{max} + jy_{max} & & ... & & & +x_{max} + jy_{max}\\ & & & & & \\ ... & & & & & \\ & & & & & \\ & & & & & \\ -x_{max} - jy_{max} & & ... & & & +x_{max} - jy_{max} \\ \end{pmatrix}
Utilizando el array t
, el valor T
y los parámetros definidos al principio (t0, alfa, U...
) crea f
según la fórmula de arriba (no hace falta que crees una función).
Guarda la parte imaginaria de esa función (función de corriente) en una variable psi
.
# Circulación que hay que añadir al cilindro para
# que se cumpla la hipótesis de Kutta en el perfil
T = 4 * np.pi * a * U * (delta + (1+landa) * alfa)
# Malla compleja
tt = XX + YY * 1j
# Potencial complejo
f = U * ( (tt - t0) * np.exp(-alfa *1j) + R**2 / (tt - t0) * np.exp(alfa * 1j) )
f += 1j * T / (2* np.pi) * np.log(tt - t0)
# Función de corriente
psi = np.imag(f)
Como la función de corriente toma un valor constante en cada línea de corriente, podemos visualizar el flujo alrededor del cilindro pintando las lineas en las que psi
toma un valor constante. Para ello utilizaremos la función contour()
en la malla XX, YY
. Si no se ve nada prueba a cambiar el número de líneas y los valores máximo y mínimo de la función que se representan.
#lo visulaizamos
plt.figure('lineas de corriente', figsize=(10,10))
plt.contour(XX, YY, psi, np.linspace(-5,5,50))
plt.grid()
plt.gca().set_aspect(1)
#ponemos el cilindro encima
plt.figure('flujo cilindro', figsize=(10,10))
plt.contour(XX, YY, psi, np.linspace(-5,5,50), colors=['blue', 'blue'])
plt.grid()
plt.gca().set_aspect(1)
p = plt.Polygon(list(zip(Xc, Yc)), color="#cccccc", zorder=3)
plt.gca().add_patch(p)
<matplotlib.patches.Polygon at 0x7f7c6d768c18>
Bueno, lo que queríamos era hacer cosas alrededor de nuestro perfil, ¿no?
Esto lo conseguiremos pintando la función psi
(tal cual está definida) en los puntos XX_tau, YY_tau
transformados de los XX, YY
a través de la función transf_yukovski
, recuerda que la transformación que tenemos recibe y da números complejos. Como antes, debes separar parte real e imaginaria. En la siguiente celda transforma tt (donde debería estar almacenada la malla en forma compleja) para obtener XX_tau, YY_tau
.
puntos_plano_tau = transf_yukovski(a, tt)
XX_tau, YY_tau = np.real(puntos_plano_tau) , np.imag(puntos_plano_tau)
#pintamos la malla para ver como se ha transformado
plt.figure("malla", figsize=(14,6))
plt.subplot(1,2,1)
plt.scatter(XX.flatten(), YY.flatten(),\
cmap=plt.cm.Paired, c=np.linspace(0,1, N_malla**2))
plt.xlim(-x_max,x_max)
plt.ylim(-y_max,y_max)
plt.subplot(1,2,2)
plt.scatter(XX_tau.flatten(), YY_tau.flatten(),\
cmap=plt.cm.Paired, c=np.linspace(0,1, N_malla**2))
plt.xlim(-x_max,x_max)
plt.ylim(-y_max,y_max)
(-3, 3)
plt.figure("flujo perfil", figsize=(12,12))
plt.contour(XX_tau, YY_tau, psi, np.linspace(-5,5,50))
plt.xlim(-8,8)
plt.ylim(-3,3)
plt.grid()
plt.gca().set_aspect(1)
Si observas bien la figura anterior, verás como si hay una serie de líneas que parecen seguir un perfil, pero que sin embargo aparecen otras que no tienen ningún sentido. Esto está apareciendo al transformarse los puntos que están dentro de la circunferencia. Lo que vamos a hacer es eliminarlos.
Para hacer esto debes recorrer uno a uno los valores del array t
y hacer cero aquellos cuya distancia a t0
sea menor de 0.95 veces el radio R
. Pista: utiliza dos bucles y un condicional. Puedes obtener la distancia con np.abs(¿¿??)
.
# Eliminamos la singularidad que se produce al evaluar los puntos
# cercanos al origen del plano t
for i in range(np.size(tt[0,:])):
for j in range(np.size(tt[:, 0])):
if np.abs(tt[i,j] - t0) < 0.95 * R:
tt[i,j] = 0
tt[np.abs(tt - t0) < 0.95 * R] = 0
Hecho esto, transforma de nuevo, el array t, para obtener XX_tau, YY_tau
.
puntos_plano_tau = transf_yukovski(a, tt)
XX_tau, YY_tau = np.real(puntos_plano_tau) , np.imag(puntos_plano_tau)
-c:6: RuntimeWarning: divide by zero encountered in true_divide -c:6: RuntimeWarning: invalid value encountered in true_divide
#pintamos la malla para ver como se ha transformado
plt.figure("malla", figsize=(14,6))
plt.subplot(1,2,1)
plt.scatter(XX.flatten(), YY.flatten(),\
cmap=plt.cm.Paired, c=np.linspace(0,1, N_malla**2))
plt.xlim(-x_max,x_max)
plt.ylim(-y_max,y_max)
plt.subplot(1,2,2)
plt.scatter(XX_tau.flatten(), YY_tau.flatten(),\
cmap=plt.cm.Paired, c=np.linspace(0,1, N_malla**2))
plt.xlim(-x_max,x_max)
plt.ylim(-y_max,y_max)
(-3, 3)
Pintamos otra vez la función psi y....
plt.figure("flujo perfil", figsize=(12,12))
plt.contour(XX_tau, YY_tau, psi, np.linspace(-5,5,50))
plt.xlim(-8,8)
plt.ylim(-3,3)
plt.grid()
plt.gca().set_aspect(1)
#Ahora ponemos el perfil encima
plt.figure("flujo perfil", figsize=(12,12))
plt.contour(XX_tau, YY_tau, psi, np.linspace(-5, 5, 50), colors=['blue', 'blue'])
p = plt.Polygon(list(zip(Xp, Yp)), color="#cccccc", zorder=10)
plt.gca().add_patch(p)
plt.xlim(-8,8)
plt.ylim(-3,3)
plt.grid()
plt.gca().set_aspect(1)
Ahora es un buen momento para jugar con todos los parámetros del problema.
¡Prueba a cambiarlos y ejecuta el notebook entero!
Vamos a usar un interact
, ¿no?
Tenemos que crear una función que haga todas las tareas: reciba los argumentos y pinte para llamar a interact con ella. No tenemos más que cortar y pegar.
def transformacion_geometrica(a, landa, delta, N_perfil=100):
#punto del plano complejo
t0 = a * (-landa + delta * 1j)
#valor del radio de la circunferencia
R = a * np.sqrt((1 + landa)**2 + delta**2)
#se barre un ángulo de 0 a 2 pi
theta = np.linspace(0, 2*np.pi, N_perfil)
#se crean las coordenadas del los puntos
#de la circunferencia
Xc = - a * landa + R * np.cos(theta)
Yc = a * delta + R * np.sin(theta)
#se crean las coordenadas del los puntos
#del perfil
Puntos_perfil = transf_yukovski(a, Xc+Yc*1j)
Xp, Yp = np.real(Puntos_perfil) , np.imag(Puntos_perfil)
#Se pintan la cirunferencia y el perfil
fig, ax = plt.subplots(1,2)
fig.set_size_inches(15,15)
p_c = plt.Polygon(list(zip(Xc, Yc)), color="#cccccc", zorder=1)
ax[0].add_patch(p_c)
ax[0].plot(Xc,Yc)
ax[0].set_aspect(1)
ax[0].set_xlim(-3, 3)
ax[0].set_ylim(-2,2)
ax[0].grid()
p_p = plt.Polygon(list(zip(Xp, Yp)), color="#cccccc", zorder=1)
ax[1].add_patch(p_p)
ax[1].plot(Xp,Yp)
ax[1].set_aspect(1)
ax[1].set_xlim(-3, 3)
ax[1].set_ylim(-2,2)
ax[1].grid()
from IPython.html.widgets import interact
w = interact(transformacion_geometrica, landa=(-1.,1), delta=(-1.,1), a=(0,2.), N_perfil=(4, 200) )
def flujo_perfil_circunferencia(landa, delta, alfa, U=1, N_malla = 100):
N_perfil=100
a=1
#punto del plano complejo
t0 = a * (-landa + delta * 1j)
#valor del radio de la circunferencia
R = a * np.sqrt((1 + landa)**2 + delta**2)
#se barre un ángulo de 0 a 2 pi
theta = np.linspace(0, 2*np.pi, N_perfil)
#se crean las coordenadas del los puntos
#de la circunferencia
Xc = - a * landa + R * np.cos(theta)
Yc = a * delta + R * np.sin(theta)
#se crean las coordenadas del los puntos
#del perfil
Puntos_perfil = transf_yukovski(a, Xc+Yc*1j)
Xp, Yp = np.real(Puntos_perfil) , np.imag(Puntos_perfil)
#se crea la malla donde se va pintar la función de corriente
x_max = 8
y_max = 3
x = np.linspace(-x_max, x_max, N_malla)
y = np.linspace(-y_max, y_max, N_malla)
XX, YY = np.meshgrid(x,y)
tt = XX + YY * 1j
tt[np.abs(tt - t0) < 0.95 * R] = 0
alfa = np.deg2rad(alfa)
# Circulación que hay que añadir al cilindro para
# que se cumpla la hipótesis de Kutta en el perfil
T = 4 * np.pi * a * U * (delta + (1+landa) * alfa)
#Potencial complejo
f = U * ( (tt - t0) * np.exp(-alfa *1j) + R**2 / (tt - t0) * np.exp(alfa * 1j) )
f += 1j * T / (2* np.pi) * np.log(tt - t0)
#Función de corriente
psi = np.imag(f)
Puntos_plano_tau = transf_yukovski(a, tt)
XX_tau, YY_tau = np.real(Puntos_plano_tau) , np.imag(Puntos_plano_tau)
#Se pinta
fig, ax = plt.subplots(1,2)
#lineas de corriente
fig.set_size_inches(15,15)
ax[0].contour(XX, YY, psi, np.linspace(-10,10,50), colors = ['blue', 'blue'])
ax[0].grid()
ax[0].set_aspect(1)
p = plt.Polygon(list(zip(Xc, Yc)), color="#cccccc", zorder=10)
ax[0].add_patch(p)
ax[0].set_xlim(-5, 5)
ax[0].set_ylim(-2,2)
ax[1].contour(XX_tau, YY_tau, psi, np.linspace(-10,10,50), colors = ['blue', 'blue'])
ax[1].grid()
ax[1].set_aspect(1)
p = plt.Polygon(list(zip(Xp, Yp)), color="#cccccc", zorder=10)
ax[1].add_patch(p)
ax[1].set_xlim(-5, 5)
ax[1].set_ylim(-2,2)
p = interact(flujo_perfil_circunferencia, landa=(0.,1), delta=(0.,1), alfa=(0, 30), U=(0,10))
En esta clase hemos reafirmado nuestros conocimientos de NumPy, matplotlib y Python, general (funciones, bucles, condicionales...) aplicándolos a un ejemplo muy aeronáutico (¿o muy ETSIA?)
Curso AeroPython por Juan Luis Cano Rodriguez y Alejandro Sáez Mollejo se distribuye bajo una Licencia Creative Commons Atribución 4.0 Internacional.
Las siguientes celdas contienen configuración del Notebook
Para visualizar y utlizar los enlaces a Twitter el notebook debe ejecutarse como seguro
File > Trusted Notebook
%%html
<a href="https://twitter.com/Pybonacci" class="twitter-follow-button" data-show-count="false">Follow @Pybonacci</a>
<script>!function(d,s,id){var js,fjs=d.getElementsByTagName(s)[0],p=/^http:/.test(d.location)?'http':'https';if(!d.getElementById(id)){js=d.createElement(s);js.id=id;js.src=p+'://platform.twitter.com/widgets.js';fjs.parentNode.insertBefore(js,fjs);}}(document, 'script', 'twitter-wjs');</script>
# Esta celda da el estilo al notebook
from IPython.core.display import HTML
css_file = '../static/styles/style.css'
HTML(open(css_file, "r").read())