Visualizando líneas de corriente en Python con matplotlib

Introducción

Hoy vamos a ver cómo representar diagramas de corriente en Python usando matplotlib. Este tipo de diagramas aparecen en Mecánica de Fluidos para visualizar el movimiento del fluido que estamos estudiando.
Hace unos días nos hicimos eco en Pybonacci de que se había liberado matplotlib 1.2.0, que introducía entre otras cosas soporte para Python 3 y la nueva función streamplot:


Así que vamos a estrenar las entradas con Python 3.3 y matplotlib 1.2 con un bonito ejemplo de Aerodinámica básica 🙂 El ejemplo y las gráficas están basados en la página de la Wikipedia sobre flujo potencial alrededor de un cilindro circular.
En esta entrada se han usado python 3.3.0, numpy 1.7.0b2 y matplotlib 1.2.0.

Campo de velocidades

Supongamos que tenemos un cilindro circular y que queremos visualizar el flujo a su alrededor. Lejos del cilindro este es unidireccional y uniforme, y nos vamos a centrar en el caso de fluido potencial y sin vorticidad. Al final, lo que tendremos será un campo de velocidades del tipo
$latex \displaystyle \vec{v} = \vec{f}(\vec{r})$
donde $latex \vec{r}$ es el vector de posición del punto considerado, y con la función streamplot lo visualizaremos en el plano. Además, incluiremos una representación del coeficiente de presiones $latex c_p$.
Sin entrar en mucho detalle matemático, la función potencial $latex \phi$ que satisface este problema, siendo $latex \vec{v} = \nabla{\phi}$, es, en coordenadas polares centradas en el cilindro,
$latex \displaystyle \phi(r, \theta) = U \left( r + \frac{R^2}{r} \right) \cos{\theta}$.
Por otro lado, la función de corriente $latex psi$ será
$latex \displaystyle \psi(r, \theta) = U \left( r – \frac{R^2}{r} \right) \sin{\theta}$
y el coeficiente de presiones $latex c_p$ será
$latex \displaystyle c_p = 2 \frac{R^2}{r^2} \cos(2 \theta) – \frac{R^4}{r^4}$.
El campo de velocidades se obtendrá derivando el potencial, de forma que será
$latex \displaystyle v_r = \frac{\partial \phi}{\partial r} = U\left(1-\frac{R^2}{r^2}\right)\cos\theta$
$latex \displaystyle v_{\theta} = \frac{1}{r}\frac{\partial \phi}{\partial \theta} = – U\left(1+\frac{R^2}{r^2}\right)\sin\theta$.
Vamos a ver cómo trasladamos todo esto a código Python.

Código Python

Para nuestro problema concreto vamos a asumir $latex R = U_{\infty} = 1$ y que nuestro dominio bidimensional es un cuadrado cuya diagonal va de $latex (-3, -3)$ a $latex (3, 3)$.
La función streamplot recibe como mínimo cuatro argumentos:

  • Dos arrays x e y, que representan el plano, y
  • dos arrays v_x y v_y, que son las componentes del campo de velocidades según x e y.

Como vemos, streamplot asume que vamos a dar el campo de velocidades en coordenadas rectangulares, así que tenemos que trabajar un poco. En primer lugar generamos el dominio plano utilizando la función meshgrid, que ya vimos en nuestro artículo sobre líneas de nivel en Python:
[sourcecode language=”python”]
x = np.linspace(-3, 3, 151)
y = np.linspace(-3, 3, 151)
xx, yy = np.meshgrid(x, y)
[/sourcecode]
Ahora lo tenemos que transformar a coordenadas polares, y vamos a tener la precacución de enmascarar la parte central del dominio para evitar singularidades. Ya empleamos los arrays enmascarados en nuestro artículo sobre estadística en Python con SciPy:
[sourcecode language=”python”]
rr = np.sqrt(xx ** 2 + yy ** 2)
tt = np.arctan2(yy, xx)
# Enmascaramos el centro para evitar singularidades
rr = ma.masked_less_equal(rr, R * 0.9)
[/sourcecode]
Ahora ya podemos calcular nuestras funciones potencial y de corriente, nuestro coeficiente de presiones y nuestras velocidades:
[sourcecode language=”python”]
# Función potencial
phi = U * (rr + R ** 2 / rr) * np.cos(tt)
# Función de corriente
psi = U * (rr – R ** 2 / rr) * np.sin(tt)
# Coeficiente de presiones
c_p = 2 * R ** 2 / rr ** 2 * np.cos(2 * tt) – R ** 4 / rr ** 4
# Velocidad (polares)
v_r = U * (1 – R ** 2 / rr ** 2) * np.cos(tt)
v_theta = -U * (1 + R ** 2 / rr ** 2) * np.sin(tt)
# Velocidad (rectangulares)
v_x = v_r * np.cos(tt) – v_theta * np.sin(tt)
v_y = v_r * np.sin(tt) + v_theta * np.cos(tt)
[/sourcecode]
Para la gráfica vamos a combinar la función streamplot para representar las líneas de corriente con la función contourf, que también mencionamos en nuestro artículo sobre líneas de nivel, para representar el coeficiente de presiones en nuestro dominio:
[sourcecode language=”python”]
# Creamos la figura
fig = plt.figure()
ax = fig.add_subplot(111)
# Cilindro
c = Circle((0, 0), R, color=’#bbbbbb’, linewidth=0, zorder=10)
ax.add_patch(c)
# Líneas de corriente
ax.streamplot(xx, yy, v_x, v_y, linewidth=0.8, color=’k’)
# Representación del coeficiente de presiones
cs = ax.contourf(xx, yy, c_p, np.linspace(-3, 1, 1000), cmap=cm.GnBu)
cb = fig.colorbar(cs, ticks=(-3 + np.arange(5)))
[/sourcecode]
El código completo lo puedes ver a través de nbviewer y este es el resultado:

Flujo potencial alrededor de un cilindro circular


Como ves, la única diferencia que hay en el código respecto a lo que habríamos escrito en Python 2.x es que no hace falta especificar que las cadenas contienen Unicode, porque este es ya el tipo por defecto.
Quisiera dedicar esta no demasiado lúcida entrada a todos mis compañeros aeronáuticos, a quienes se les echa mucho de menos desde la distancia 🙂
¡Un saludo!

11 pensamientos sobre “Visualizando líneas de corriente en Python con matplotlib”

  1. Pingback: Un 2012 de Python científico « Pybonacci

  2. interesante Juan
    pero para que corra el programa sin error falta el sigiente codigo
    import numpy as np
    import numpy.ma as ma
    import matplotlib.pyplot as plt
    from matplotlib.patches import Circle
    import matplotlib.cm as cm

    plt.show()

    1. Ya descubrí mi error, lo que me pasaba era que graficaba la punta de las flechas pero no las líneas de corriente. Esto era porque debía graficar las presiones antes, para que las líneas aparecieran sobre las presiones.
      Saludos! y muy buena página

Deja una respuesta

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *

− three = two