Texto y código sujeto bajo Creative Commons Attribution license, CC-BY-SA. (c) Original por Lorena A. Barba, 2013, traducido libremente por F.J. Navarro-Brull para CAChemE.org
Version 0.1 (Julio 2013)

12 pasos para Navier-Stokes


¡Hola! Bienvenidos a 12 pasos para Navier-Stokes . Este es un módulo práctico que se utiliza en el principio de un curso interactivo de Dinámica de Fluidos Computacional (en inglés Computational Fluid Dynamics, CFD) y se imparte por la Prof. Lorena Barba desde el 2009 en la Universidad de Boston. El curso supone sólo el conocimiento básico de programación (en cualquier lenguaje) y, por supuesto, algún fundamento en las ecuaciones en derivadas parciales y mecánica de fluidos. El módulo práctico se inspiró en las ideas del Dr. Rio Yokota, que era un post-doctorado en el laboratorio de Barba, y ha sido depurado por la profesora Barba y sus alumnos durante varios semestres de enseñanza del curso. El curso se imparte íntegramente utilizando Python y los estudiantes que no saben Python simplemente aprenderán a medida que se trabaje a través del módulo.

Este IPython notebook será la primera etapa de programación de tu propio resolvedor de Navier-Stokes en Python. Vamos a ello ahora mismo. No te preocupes si no entiendes todo lo que está pasando en primer lugar, cubriremos en detalle cada aspecto a medida que se avance y mediante los videos de apoyo (en inglés) (EN) Clases de la Prof. Barba en YouTube.

Para obtener los mejores resultados, una vez visto este notebook (cuaderno), escribe tu propio código para el Paso 1, ya sea como un script en Python o en un notebook IPython vacío.

Para ejecutar este cuaderno, se supone que has arrancado el servidor de notebook usando: ipython notebook --pylab inline en tu línea de comandos.

Paso 1: 1-D Convección lineal (Linear Convection)


La ecuación en 1-D de convección lineal (1-D Linear Convection) es la ecuación más simple que permite el modelado más sencillo y permite aprender bastante sobre CFD (Computational fluid dynamics). ¡Es sorprendente como esta pequeña ecuación puede enseñarnos tanto! aquí está:

$$\frac{\partial u}{\partial t} + c \frac{\partial u}{\partial x} = 0$$

Con unas condiciones iniciales dadas (entendido como onda o wave), la ecuación representa la propagación de esta onda inicial con velocdad $c$, sin cambio de su forma. Siendo la condición inicial $u(x,0)=u_0(x)$. Entonces, la solución exacta de la ecuación es $u(x,t)=u_0(x-ct)$.

Para discretizar esta ecuación en tiempo y espacio, usamos el método de diferencias finitas, hacia delante para la derivada de tiempo y hacia atrás para el espacio. Así, se discretiza la coordenada espacial $x$ en puntos con un índice de $i=0$ a $N$, y el tiempo en intervalos de un tamaño específico $\Delta t$.

De la definición de una derivada (simplemente quitando el límite), sabemos que:

$$\frac{\partial u}{\partial x}\approx \frac{u(x+\Delta x)-u(x)}{\Delta x}$$

Nuestra ecuación discreta queda como:

$$\frac{u_i^{n+1}-u_i^n}{\Delta t} + c \frac{u_i^n - u_{i-1}^n}{\Delta x} = 0 $$

Donde $n$ y $n+1$ son dos puntos consecutivos en el tiempo, mientras que, $i-1$ y $i$ son dos puntos vecinos de la coordenada $x$ discretizada. Si se dan las condiciones iniciales, entonces la única incógnita en esta discretización es $u_i^{n+1}$. Podemos resolver la incógnita y obtener una ecuación que nos permite avanzar en el tiempo, de la siguiente manera:

$$u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)$$

Ahora vamos a tratar de implementar esto en Python.

Vamos a empezar importando algunas librerías para ayudarnos.

  • numpy es una biblioteca que proporciona un grupo de operaciones matriciales de forma similar a MATLAB
  • matplotlib es una biblioteca de gráficas 2D que vamos a utilizar para representar los resultados
  • time y sys proporcionan funciones básicas de tiempo que vamos a utilizar para reducir la velocidad de visualización de animaciones
In [1]:
# Recuerda: Los comentarios en python se indican con el signo de almohadilla

%pylab inline
# El comando de arriba hará que figuras de este notebook se representen junto al texto

import numpy as np                 # aquí importamos numpy, llamándolo "np" de ahora en adelante
import matplotlib.pyplot as plt    # aquí importamos matplotlib, llamándolo "plt"
import time, sys                   # e importamos algunas funciones útiles
from IPython.core.display import clear_output # utilizado para la animación que se verá adelante
Populating the interactive namespace from numpy and matplotlib

Ahora vamos a definir algunas variables, queremos crear una red de puntos espaciados uniformemente dentro de un dominio espacial que tenga 2 unidades de longitud de ancho, es decir, $x_i\in(0,2)$. Vamos a definir una variable nx, que será el número de puntos de la cuadrícula que queremos y dx será la distancia entre cualquier par de puntos de las cuadrículas adyacentes.

In [2]:
nx = 41  # prueba a cambiar este número de 41 a 81 y ejecutar todo ... ¿qué pasa?
dx = 2./(nx-1)
nt = 25    # nt es el número de pasos de tiempo que queremos calcular
dt = .025  # dt es la cantidad de tiempo que cada momento de tiempo abarca (delta t)
c = 1.     # asumir la velocidad de la onda de c = 1

También tenemos que establecer nuestras condiciones iniciales. La velocidad inicial $u_0$ viene dada por $u = 2$ en el intervalo $0.5 \leq x \leq 1$ y $u = 1$ en el resto de $(0,2)$ (ésto es la función de sombrero).

Aquí, utilizamos la función ones() definiendo un 'array' de numpy que tiene nx elementos de largo con cada valor igual a 1.

In [12]:
u = np.ones(nx)      # numpy función ones() creando un array de 1 con nx elmentos
u[.5/dx : 1/dx+1]=2  # estableciendo u = 2 entre 0.5 y 1 como nuestras C.I.s (condiciones de iniciales)
print u
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  2.  2.  2.  2.  2.  2.  2.  2.
  2.  2.  2.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.  1.  1.  1.]

Ahora echemos un vistazo a las condiciones iniciales representándolo gráficamente con Matplotlib. Hemos importado matplotlib como plt y, dentro de ella, la función representación de datos plot, por lo que llamamos a plt.plot. Para saber más sobre las múltiples posibilidades de Matplotlib, échale un vistazo a la Galería de ejemplos de gráficas.

Aquí, utilizamos la sintaxis de una simple gráfica 2D: plot (x, y), donde los valores de x se distribuyen uniformemente en puntos de la cuadrícula:

In [13]:
plt.plot(np.linspace(0,2,nx), u)
Out[13]:
[<matplotlib.lines.Line2D at 0x7fc77856e190>]

¿Por qué la función del sombrero (hat function) tiene las líneas perfectamente rectas? Piénsalo un poco.

Ahora es el momento de poner en práctica la discretización de la ecuación de convección usando el método de diferencias finitas.

Para cada elemento de nuestra matriz u, necesitamos realizar la operación $u_i^{n+1} = u_i^n - c \frac{\Delta t}{\Delta x}(u_i^n-u_{i-1}^n)$

Vamos a guardar el resultado en un nuevo array (temporal), un, que será la solución de $u$ para el siguiente punto en el tiempo. Vamos a repetir esta operación para todos los intervalos de tiempo que especifiquemos y, entonces, podemos ver hasta qué punto la onda experimenta convección.

Primero inicializamos nuestro array marcador de posición, un, y mantenemos los valores que calculados para el $n+1$ intervalos de tiempo, utilizando una vez más la función de NumPy ones().

Entonces, podemos pensar que tenemos dos operaciones iterativas: una en el espacio y una en el tiempo (aprenderemos una manera diferente de hacerlo después), por lo que vamos a empezar por la inclusión o anidadamiento de un bucle dentro de otro. Ten en cuenta el uso ingenioso de la función range(). Cuando escribimos: for i in range (1, nx) vamos a recorrer el array u, pero vamos a saltar al primer elemento (el elemento cero). ¿Por qué?

In [14]:
un = np.ones(nx) # crea un array temporal de unos con nx elementos

for n in range(nt):  # bucle para los valores de n desde 0 a nt-1, por lo que se ejecutará nt veces
    un[:] = u[:] # copia los valores existentes de u en un
    for i in range(1,nx): ## puedes convertir esta línea en un comentario...
    #for i in range(nx): ## ... e intercambiarla por esta otra para ver que pasa!
        u[i] = un[i]-c*dt/dx*(un[i]-un[i-1])
        
        

Nota: Aprenderemos más tarde que el código escrito arriba es bastante ineficiente y hay mejores maneras de escribir esto al estilo Python. Pero vamos a continuar de momento así.

Ahora vamos a intentar representar nuestro array u después de que el tiempo haya avanzado.

In [15]:
plt.plot(np.linspace(0,2,nx),u)
Out[15]:
[<matplotlib.lines.Line2D at 0x7fc778540c10>]

¡OK! Nuestra función sombrero (de copa) se ha movido definitivamente hacia la derecha, pero ya no es un sombrero "recto". ¿Qué está pasando?

Aprende más


Para una explicación más detallada del método de diferencias finitas, incluyendo temas como el error de truncamiento, el orden de convergencia y otros detalles, puedes ver en ingés las Video Lessons 2 and 3 de la profesora Barba en YouTube.

In [7]:
from IPython.display import YouTubeVideo
YouTubeVideo('iz22_37mMkk')
Out[7]:
In [8]:
YouTubeVideo('xq9YTcv-fQg')
Out[8]:

Para una explicación detallada (paso por paso) de la discretización de la ecuación de convección lineal con diferencias finitas (y también los siguientes pasos de estos notebooks hasta el Paso 4), mira el vídeo Lesson 4 de la profesora Barba en YouTube.

In [9]:
YouTubeVideo('y2WaK7_iMRI')
Out[9]:

Por último pero no por ello menos importante (Last but not least)

Recuerda que deberías reescribir el código de este notebook desde cero como script en Python o experimentar con el cambiando los valores de discretización mediante tu instalación de IPython notebook. Una vez hecho esto, estás listo para el Paso 2.


In [2]:
from IPython.core.display import HTML
def css_styling():
    styles = open("../styles/custom.css", "r").read()
    return HTML(styles)
css_styling()
Out[2]:

(La celda de arriba establece el estilo de este notebook. Se modificó el estilo encontrado en GitHub por CamDavidsonPilon, @Cmrn_DP.)