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: https://twitter.com/Pybonacci/status/266798396654882816 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: [caption id="attachment_1219" align="aligncenter" width="560"] Flujo potencial alrededor de un cilindro circular[/caption] 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!