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

12 pasos a Navier-Stokes


Continuamos nuestro viaje para resolver la ecuación de Navier-Stokes con el Paso 4. ¡Pero no continues a menos que hayas completado los pasos anteriores! De hecho, el próximo paso será una combinación de los dos anteriores. Las maravillas de la reutilización de código!

Paso 4: Ecuación de Burgers


Puedes leer sobre la ecuación de Burgers en la wikipedia (inglés).

La ecuación de burgers con una dimensión espacial tiene la siguiente forma:

$$\frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} = \nu \frac{\partial ^2u}{\partial x^2}$$

Como puedes ver, es una combiniación de convección y difusión no lineal. Es sorprendente lo mucho que se puede aprender de esta pequeña y sencilla ecuación.

Podemos discretizarla usando los metodos que se han detallado de los Pasos 1 a 3. Usando diferencias hacia delante para el tiempo, hacia atrás para el espacio y nuestro método de 2º orden para la segunda derivada, se llega a:

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

Como antes, una vez dadas las condiciones iniciales, la única incógina es $u_i^{n+1}$. Despejando, haremos avanzar el tiempo de la siguiente forma:

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

Condiciones iniciales y de contorno (Initial and Boundary Conditions)

Para examinar algunas propiedades interesantes de la ecuación de Burgers, es útil usar diferentes condiciones iniciales y de contorno de las que hemos estado usando hasta ahora.

Nuestra condición inicial para este problema va a ser:

\begin{eqnarray} u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\ \phi &=& \exp \bigg(\frac{-x^2}{4 \nu} \bigg) + \exp \bigg(\frac{-(x-2 \pi)^2}{4 \nu} \bigg) \end{eqnarray}

Esto tiene una solución analítica, dada por:

\begin{eqnarray} u &=& -\frac{2 \nu}{\phi} \frac{\partial \phi}{\partial x} + 4 \\\ \phi &=& \exp \bigg(\frac{-(x-4t)^2}{4 \nu (t+1)} \bigg) + \exp \bigg(\frac{-(x-4t -2 \pi)^2}{4 \nu(t+1)} \bigg) \end{eqnarray}

Nuestra condición de frontera será:

$$u(0) = u(2\pi)$$

Esto se llama condición de frontera periódica. ¡Presta atención! Esto dará un poco de dolor de cabeza si no se anda con cuidado.

Ahorrando tiempo con SymPy

La condición inicial que estamos utilizando para la ecuación de Burgers puede ser un poco compleja de evaluar a mano. La derivada $\frac{\partial \phi}{\partial x}$ no es terriblemente difícil, pero sería fácil cambiar un signo u olvidar un factor de $x$ en alguna parte, así que vamos a utilizar SymPy para ayudarnos.

SymPy es la librería matemática de lenguaje simbólico para Python. Tiene muchas de las mismas funciones matemáticas que Mathematica con la ventaja añadida de que podemos traducir fácilmente sus resultados a nuestros cálculos en Python (también es gratuita y de código abierto).

Comenzamos cargando la librería SymPy junto con nuestra biblioteca favorita NumPy.

In [2]:
import numpy as np
import sympy

También vamos a decir a SymPy que queremos toda la salida de resultados presentados mediante $\LaTeX$. ¡Esto hará que nuestro Notebook se vea bonito!

In [3]:
from sympy import init_printing
init_printing(use_latex=True)

Comenzaremos por la creación de variables simbólicas para las tres variables en nuestra condición inicial y escribiremos la ecuación completa para $\phi$. Debemos conseguir una versión renderizada de nuestra ecuación $\phi$.

In [4]:
x, nu, t = sympy.symbols('x nu t')
phi = sympy.exp(-(x-4*t)**2/(4*nu*(t+1))) + sympy.exp(-(x-4*t-2*np.pi)**2/(4*nu*(t+1)))
phi
Out[4]:
$$e^{- \frac{\left(- 4 t + x - 6.28318530717959\right)^{2}}{4 \nu \left(t + 1\right)}} + e^{- \frac{\left(- 4 t + x\right)^{2}}{4 \nu \left(t + 1\right)}}$$

Es quizás un poco pequeña, pero se puede ver. Ahora, evaluar nuestro derivada parcial $\frac{\partial \phi}{\partial x}$ es algo trivial.

In [5]:
phiprime = phi.diff(x)
phiprime
Out[5]:
$$- \frac{e^{- \frac{\left(- 4 t + x\right)^{2}}{4 \nu \left(t + 1\right)}}}{4 \nu \left(t + 1\right)} \left(- 8 t + 2 x\right) - \frac{1}{4 \nu \left(t + 1\right)} \left(- 8 t + 2 x - 12.5663706143592\right) e^{- \frac{\left(- 4 t + x - 6.28318530717959\right)^{2}}{4 \nu \left(t + 1\right)}}$$

Si quieres ver la versión sin renderizar, sólo tienes que utilizar el comando de impresión (print) de Python.

In [16]:
print phiprime
-(-8*t + 2*x)*exp(-(-4*t + x)**2/(4*nu*(t + 1)))/(4*nu*(t + 1)) - (-8*t + 2*x - 12.5663706143592)*exp(-(-4*t + x - 6.28318530717959)**2/(4*nu*(t + 1)))/(4*nu*(t + 1))

¿Ahora qué?

Ahora que ya tenemos la versión "Pythonic" de nuestra derivada, podemos terminar de escribir la ecuación de condición inicial y luego traducirlo en una expresión Python utilizable. Para ello, vamos a utilizar la función lambdify, que a partir una ecuación simbólica SymPy la convierte en una función que se puede ejecutar.

In [6]:
from sympy.utilities.lambdify import lambdify

u = -2*nu*(phiprime/phi)+4
u
Out[6]:
$$- \frac{2 \nu}{e^{- \frac{\left(- 4 t + x - 6.28318530717959\right)^{2}}{4 \nu \left(t + 1\right)}} + e^{- \frac{\left(- 4 t + x\right)^{2}}{4 \nu \left(t + 1\right)}}} \left(- \frac{e^{- \frac{\left(- 4 t + x\right)^{2}}{4 \nu \left(t + 1\right)}}}{4 \nu \left(t + 1\right)} \left(- 8 t + 2 x\right) - \frac{1}{4 \nu \left(t + 1\right)} \left(- 8 t + 2 x - 12.5663706143592\right) e^{- \frac{\left(- 4 t + x - 6.28318530717959\right)^{2}}{4 \nu \left(t + 1\right)}}\right) + 4$$

Lambdify

Para transfomar con lambdify la expresión en una función que se pueda usar, decimos a lambdify qué variables son las de entradas (t,x,nu) y las de salida (u).

In [7]:
ufunc = lambdify((t, x, nu), u)
ufunc(1,4,3)
Out[7]:
$$3.4917066420644494$$

De vuelta a la ecuación de Burgers

Ahora que tenemos las condiciones iniciales establecidas, podemos continuar y finalizar la solución del problema. Podemos representar la condición inicial utilizando nuestro lambdify

In [8]:
import matplotlib.pyplot as plt

### Declaración de variables
nx = 101
nt = 100
dx = 2*np.pi/(nx-1)
nu = .07
dt = dx*nu

x = np.linspace(0, 2*np.pi, nx)
#u = np.empty(nx)
un = np.empty(nx)
t = 0

u = np.asarray([ufunc(t, x0, nu) for x0 in x])
u
Out[8]:
array([ 4.        ,  4.06283185,  4.12566371,  4.18849556,  4.25132741,
        4.31415927,  4.37699112,  4.43982297,  4.50265482,  4.56548668,
        4.62831853,  4.69115038,  4.75398224,  4.81681409,  4.87964594,
        4.9424778 ,  5.00530965,  5.0681415 ,  5.13097336,  5.19380521,
        5.25663706,  5.31946891,  5.38230077,  5.44513262,  5.50796447,
        5.57079633,  5.63362818,  5.69646003,  5.75929189,  5.82212374,
        5.88495559,  5.94778745,  6.0106193 ,  6.07345115,  6.136283  ,
        6.19911486,  6.26194671,  6.32477856,  6.38761042,  6.45044227,
        6.51327412,  6.57610598,  6.63893783,  6.70176967,  6.76460125,
        6.82742866,  6.89018589,  6.95176632,  6.99367964,  6.72527549,
        4.        ,  1.27472451,  1.00632036,  1.04823368,  1.10981411,
        1.17257134,  1.23539875,  1.29823033,  1.36106217,  1.42389402,
        1.48672588,  1.54955773,  1.61238958,  1.67522144,  1.73805329,
        1.80088514,  1.863717  ,  1.92654885,  1.9893807 ,  2.05221255,
        2.11504441,  2.17787626,  2.24070811,  2.30353997,  2.36637182,
        2.42920367,  2.49203553,  2.55486738,  2.61769923,  2.68053109,
        2.74336294,  2.80619479,  2.86902664,  2.9318585 ,  2.99469035,
        3.0575222 ,  3.12035406,  3.18318591,  3.24601776,  3.30884962,
        3.37168147,  3.43451332,  3.49734518,  3.56017703,  3.62300888,
        3.68584073,  3.74867259,  3.81150444,  3.87433629,  3.93716815,  4.        ])
In [9]:
plt.figure(figsize=(11,7), dpi=100)
plt.plot(x,u, marker='o', lw=2)
plt.xlim([0,2*np.pi])
plt.ylim([0,10])
Out[9]:
$$\begin{pmatrix}0, & 10\end{pmatrix}$$

Esta definitivamente no es la función del sombrero que hemos estado tratando hasta ahora. Nosotros lo llamamos una "función de diente de sierra" (saw-tooth function). Vamos a continuar y ver qué pasa.

Nota: Para ejecutar este cuaderno, se supone que tienes en ejecución el servidor de notebook usando: ipython notebook --pylab inline. Para ello, basta con poner el comando anterior en la ventana de comandos (cmd.exe o terminal si estas en windows, linux o mac).

Condiciones de contorno periódicas

Una de las grandes diferencias entre el paso 4 y las lecciones anteriores es el uso de las condiciones de contorno periódicas. Si haces experimentos con los Pasos 1 y 2 y realizas la simulación durante más tiempo (aumentando nt) te darás cuenta de que la onda seguirá moviéndose hacia la derecha hasta que ya ni siquiera aparece en la malla.

Con condiciones de contorno periódicas, cuando un punto llega al lado derecho del cuadro, que retorna de vuelta a la parte izquierda del mismo.

Recordemos la discretización que nos salió a principios de este notebook:

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

Qué significado tiene $u_{i+1}^n$ cuando $i$ llega al final del borde?

Piensa en ello durante un minuto antes de continuar.

In [10]:
for n in range(nt):
    un[:]=u[:]
    for i in range(nx-1):
        u[i] = un[i] - un[i] * dt/dx *(un[i] - un[i-1]) + nu*dt/dx**2*\
                (un[i+1]-2*un[i]+un[i-1])
    u[-1] = un[-1] - un[-1] * dt/dx * (un[-1] - un[-2]) + nu*dt/dx**2*\
                (un[0]-2*un[-1]+un[-2])
        
u_analytical = np.asarray([ufunc(nt*dt, xi, nu) for xi in x])
In [11]:
plt.figure(figsize=(11,7), dpi=100)
plt.plot(x,u, marker='o', lw=2, label='Computational')
plt.plot(x, u_analytical, label='Analytical')
plt.xlim([0,2*np.pi])
plt.ylim([0,10])
plt.legend()
Out[11]:
<matplotlib.legend.Legend at 0x7fdc0bbef910>

¿Y ahora qué?

Los pasos siguientes, del 5 al 12, serán en dos dimensiones. Pero es fácil de ampliar las fórmulas de diferencias finitas de 1D a las derivadas parciales en 2D o 3D. Basta con aplicar la definición -una derivada parcial con respecto a $x$ es la variación en la dirección $x$ manteniendo $y$ constante.

Antes de pasar al Paso 5, asegúrate de que has completado tu propio código de los Pasos 1 al 4 y has jugado con los parámetros pensando sobre lo que está sucediendo. Además, te recomendamos tomar un ligero descanso para aprender sobre operaciones matriciales (arrays) con NumPy.

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

(La celda de arriba establece el formato de este notebook)