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.
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!
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)$$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.
%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
[<matplotlib.lines.Line2D at 0x7fe9071d7450>]
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.
from IPython.display import YouTubeVideo
YouTubeVideo('y2WaK7_iMRI')
from IPython.core.display import HTML
def css_styling():
styles = open("../styles/custom.css", "r").read()
return HTML(styles)
css_styling()
(La celda de arriba establece el formato de este notebook)