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

12 pasos para Navier-Stokes


Ves dónde va esto... haremos la difusión en 2D ahora y en el próximo vamos a combinar los pasos 6 y 7 para resolver la ecuación de Burgers. Así que asegurate de que tus pasos anteriores funcionan bien antes de continuar.

Paso 7: Difusion en 2D


Y aqui está la ecuación de difusión en 2D:

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

Recordarás que se nos ocurrió un método para discretizar derivadas de segundo orden en el paso 3, concretamente en la investigación de la difusión en 1D. Vamos a utilizar el mismo esquema aquí, con nuestras diferencias hacia delante en el tiempo y dos derivadas de segundo orden.

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

Una vez más, reorganizamos la ecuación de datos discretos y resolvemos para $u_{i,j}^{n+1}$

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

In [1]:
%pylab inline 
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D ##library for 3d projection plots
from matplotlib import cm ##cm = "colormap" for changing the 3d plot color palette

### Declaración de variables
nx = 31
ny = 31
nt = 17
nu=.05
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
sigma = .25
dt = sigma*dx*dy/nu

x = np.linspace(0,2,nx)
y = np.linspace(0,2,ny)

u = np.ones((ny,nx)) ## crea un vector 1xn de unos
un = np.ones((ny,nx)) ##

### Asignar variables inciales

u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2 ## establece función de sobrero como C.I. : u(.5<=x<=1 && .5<=y<=1 ) is 2

fig = plt.figure()
ax = fig.gca(projection='3d')
X,Y = np.meshgrid(x,y)
surf = ax.plot_surface(X,Y,u[:], rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=False)
plt.show()
ax.set_xlim(1,2)
ax.set_ylim(1,2)
ax.set_zlim(1,2.5)
#ax.zaxis.set_major_locator(LinearLocator(5))
Populating the interactive namespace from numpy and matplotlib
Out[1]:
(1, 2.5)

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

In [2]:
### Función que ejecuta la solución hasta un tiempo nt
def diffuse(nt):
    u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2
    
    for n in range(nt+1): 
        un[:] = u[:]
        u[1:-1,1:-1]=un[1:-1,1:-1]+nu*dt/dx**2*(un[2:,1:-1]-2*un[1:-1,1:-1]+un[0:-2,1:-1])+nu*dt/dy**2*(un[1:-1,2:]-2*un[1:-1,1:-1]+un[1:-1,0:-2])
        u[0,:]=1
        u[-1,:]=1
        u[:,0]=1
        u[:,-1]=1

    
    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X,Y,u[:], rstride=1, cstride=1, cmap=cm.coolwarm,
        linewidth=0, antialiased=True)
    ax.set_zlim(1,2.5)
    plt.show()
    
In [3]:
diffuse(10)
In [4]:
diffuse(14)
In [5]:
diffuse(50)

Aprender más

La lección del vídeo que te lleva a través de los detalles de los pasos 5 al 8 es la Video Lesson 6 en You Tube:

In [1]:
from IPython.display import YouTubeVideo
YouTubeVideo('tUg_dE3NXoY')
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)