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


Deberías de haber completado tu propio código hasta el Paso 5 antes de continuar con esta lección. Al igual que con los pasos 1 a 4, vamos a avanzar gradualmente, por lo que es importante completar el paso anterior!

Seguimos...

Paso 6: Convección 2D


Ahora resolveremos un problema de convección 2D, representado por el par de ecuaciones diferenciales parciales acopladas que se presentan a continuación:

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

Discretizando las ecuaciones con los métodos que hemos aplicado anteriormente se obtiene:

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

Reorganización de las dos ecuaciones, se resuelve por $u_{i,j}^{n+1}$ and $v_{i,j}^{n+1}$, respectivamente. Ten en cuenta que estas ecuaciones también están acopladas.

Condiciones iniciales

Las condiciones iniciales son las mismos que se utilizaron para la convección en 1D, aplicadas tanto en la dirección x como en la dirección y.

$$u,\ v\ = \begin{cases}\begin{matrix} 2 & \text{para } x,y \in (0.5, 1)\times(0.5,1) \cr 1 & \text{en el resto de puntos} \end{matrix}\end{cases}$$

Condiciones de contorno

Las condiciones de contorno son u y v igual a 1 a lo largo de los límites de la malla (cuadrícula).

$$u = 1,\ v = 1 \text{ para } \begin{cases} \begin{matrix}x=0,2\cr y=0,2 \end{matrix}\end{cases}$$
In [1]:
%pylab inline 
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

### Declaración de variables
nx = 101
ny = 101
nt = 80
c = 1
dx = 2.0/(nx-1)
dy = 2.0/(ny-1)
sigma = .2
dt = sigma*dx

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

u = np.ones((ny,nx)) ## crea un vector 1xn de unos
v = np.ones((ny,nx))
un = np.ones((ny,nx))
vn = 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
v[.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

for n in range(nt+1): ## Bucle para los incrementos de tiempo
    un[:] = u[:]
    vn[:] = v[:]

    u[1:,1:]=un[1:,1:]-(un[1:,1:]*dt/dx*(un[1:,1:]-un[0:-1,1:]))-vn[1:,1:]*dt/dy*(un[1:,1:]-un[1:,0:-1]) 
    v[1:,1:]=vn[1:,1:]-(un[1:,1:]*dt/dx*(vn[1:,1:]-vn[0:-1,1:]))-vn[1:,1:]*dt/dy*(vn[1:,1:]-vn[1:,0:-1])
    
    u[0,:] = 1
    u[-1,:] = 1
    u[:,0] = 1
    u[:,-1] = 1
    
    v[0,:] = 1
    v[-1,:] = 1
    v[:,0] = 1
    v[:,-1] = 1
Populating the interactive namespace from numpy and matplotlib
In [2]:
from matplotlib import cm ##cm = "colormap" para cambiar la paleta de color de la figura 3d 
fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
X,Y = np.meshgrid(x,y)

ax.plot_surface(X,Y,u, cmap=cm.coolwarm)
Out[2]:
<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7f41c1cac710>

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 YouTube:

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)