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 para Navier-Stokes


Debes de haber completado los pasos 1 y 2 antes de continuar. Este notebook de IPython continúa la presentación de los 12 pasos para Navier-Stokes, el módulo práctico se enseña en la clase interactiva de CFD impartida por la Prof. Lorena Barba.

Paso 3: Ecuación de difusión en 1-D


La ecuación de difusión unidimensional es:

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

Lo primero que deberías notar es es que, a diferencia de las dos anteriores ecuaciones simples que hemos estudiado, esta ecuación tiene una derivada de segundo orden. ¡Primero tenemos que aprender qué hacer con ella!

Discretizando $\frac{\partial ^2 u}{\partial x^2}$

La derivada de segundo orden se puede representar geométricamente como la recta tangente a la curva dada por la primera derivada. Vamos a discretizar la derivada de segundo orden con un esquema de diferencias finitas centradas. Esto lo conseguimos mediante una combinación de diferencias hacia delante y hacia atrás de la primera derivada. Considera la expansión de Taylor de $u_{i+1}$ y $u_{i-1}$ en torno a $u_i$:

$u_{i+1} = u_i + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i + \frac{\Delta x^3}{3} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)$

$u_{i-1} = u_i - \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2} \frac{\partial ^2 u}{\partial x^2}\bigg|_i - \frac{\Delta x^3}{3} \frac{\partial ^3 u}{\partial x^3}\bigg|_i + O(\Delta x^4)$

Si sumamos estas dos expansiones, se puede ver que los términos de derivadas impares se anulan entre sí. Despreciando los términos de $O(\Delta x^4)$ o superior (y en realidad, estos son muy pequeños), podemos reordenar la suma de estas dos expansiones para resolver nuestra segunda derivada.

$u_{i+1} + u_{i-1} = 2u_i+\Delta x^2 \frac{\partial ^2 u}{\partial x^2}\bigg|_i + O(\Delta x^4)$

Reordenando para resolver $\frac{\partial ^2 u}{\partial x^2}\bigg|_i$ el resultado es:

$$\frac{\partial ^2 u}{\partial x^2}=\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$

Volviendo al paso 3

Ahora podemos escribir la versión discretizada de la ecuación de difusión en 1D:

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

Al igual que antes, nos damos cuenta de que una vez que tenemos una condición inicial, la única incógnita es $u_{i}^{n+1}$, por lo que volvemos a reorganizar la ecuación resolviendo para nuestra incógnita:

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

La ecuación discreta anterior nos permite escribir un programa para avanzar en una solución en el tiempo. Pero necesitamos una condición inicial. Vamos a continuar con nuestra favorita: la función de sombrero. Por lo tanto, en $t=0$, $u=2$ en el intervalo $0.5\le x\le 1$ y $u=1$ para todo lo demás. Estamos listos para empezar a calcular números.

In [1]:
%pylab inline
# El comando de arriba hará que figuras de este notebook se representen junto al texto

import numpy as np                 # importando nuestra librería favorita
import matplotlib.pyplot as plt    # y la útil librería para representar gráficas 2D


nx = 41
dx = 2./(nx-1)
nt = 20    # número de intervalos de tiempo que se desea calcular
nu = 0.3   # valor de la viscosity
sigma = .2 #sigma es un parametro, aprenderemos más sobre ello luego
dt = sigma*dx**2/nu #d t es definido usando el valor sigma, aprenderemos más sobre ello enseguida!


u = np.ones(nx)      # un array de numpy con nx elementos iguales a 1
u[.5/dx : 1/dx+1]=2  # estableciendo u = 2 entre 0.5 y 1 paras las condiciones iniciales (I.C.s)

un = np.ones(nx) # nuestro array marcador de posición our placeholder array, un, to advance the solution in time

for n in range(nt):  # iteración a través del tiempo
    un[:] = u[:] ## copia los valores existentes de 'u' en 'un'
    for i in range(1,nx-1):
        u[i] = un[i] + nu*dt/dx**2*(un[i+1]-2*un[i]+un[i-1])
        
plt.plot(np.linspace(0,2,nx), u)
Populating the interactive namespace from numpy and matplotlib
Out[1]:
[<matplotlib.lines.Line2D at 0x7fe9071d7450>]

Aprende más

Para una explicación paso a paso de la discretización de la ecuación de difusión y todos los Pasos del 1 al 4, puedes ver el Video Lesson 4 de la profesora Barba en YouTube.

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

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