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


Hasta ahora, todo nuestro trabajo ha sido en una dimensión espacial "x" (Pasos 1 al 4). Podemos aprender mucho con sólo 1D, pero vamos a crecer hasta flatland: el mundo de las dos dimensiones.

En los siguientes ejercicios, se extienden los primeros pasos descritos (del 1 al 4) a 2D. Para ampliar las fórmulas de diferencias finitas de derivadas parciales de 1D a 2D o 3D, sólo se aplica la definición: una derivada parcial con respecto a $x$ es la variación en la dirección $x$ manteniendo constante $y$.

En el espacio 2D, una malla rectangular (uniforme) se define por los puntos con coordenadas:

$$x_i = x_0 +i \Delta x$$

$$y_i = y_0 +i \Delta y$$

Ahora, define $u_{i,j} = u(x_i,y_j)$ y aplica las fórmulas de diferencias finitas en cada variable $x,y$ actuando de forma separada en los indices $i$ y $j$. Todas las derivadas están basadas en la expansión de Taylor EN 2D de un punto de la malla con el valor aproximado de $u_{i,j}$.

Por lo tanto, para una derivada parcial de primer orden en la dirección $x$, la fórmula de diferencias finitas es:

$$ \frac{\partial u}{\partial x}\biggr\rvert_{i,j} = \frac{u_{i+1,j}-u_{i,j}}{\Delta x}+\mathcal{O}(\Delta x)$$

y de forma similar en la dirección $y$. Así, podemos escribir la forma diferencial con diferencias finitas hacia atrás, hacia delante o centradas para los siguientes pasos. ¡Vamos a ello!

Paso 5: Convección lineal en 2D.


La ecuación diferencial parcial (EDP) que gobierna la convección en 2D es:

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

Esta es exáctamente la misma fórmula que con la de Convección lineal en 1D, salvo por el hecho de que ahora tenemos dos dimensiones espaciales ($x$ e $y$) a tener en cuenta a la hora de avanzar en nuestra solución en el tiempo.

Una vez más, el paso de tiempo se puede discretizar como una diferencia hacia adelante y los incrementos espaciales se pueden discretizar como diferencias hacia atrás.

Con las implementaciones en 1D, utilizamos subíndices $i$ para denotar movimiento en el espacio (por ejemplo, $u_{i}^n-u_{i-1}^n$). Ahora que tenemos dos dimensiones para tener en cuenta, hay que añadir un segundo subíndice, $j$, para contabilizar toda la información en el espacio.

Aquí, vamos a volver a utilizar $i$ como el índice de nuestros valores $x$, y vamos a añadir los subíndices $j$ para nuestros valores de $y$.

Con esto en mente, nuestra discretización de la EDP será obtenida de forma directa.

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

Como antes, despejando para la única incógnita:

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

Resolveremos esta ecuación con las siguientes condiciones iniciales

$$u(x) = \begin{cases} \begin{matrix} 2\ \text{para} & 0.5 \leq x \leq 1 \cr 1\ \text{para} & \text{cualquier otro sitio}\end{matrix}\end{cases}$$

y condiciones de contorno:

$$u = 1\ \text{para } \begin{cases} \begin{matrix} x = 0,\ 2 \cr y = 0,\ 2 \end{matrix}\end{cases}$$

In [1]:
#Representa las figuras en este mismo notebook:
%pylab inline

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D    ## Nueva librería requerida para proyecciones en 3D

import numpy as np

### Declaración de variables
nx = 81
ny = 81
nt = 100
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
un = np.ones((ny,nx)) ##

### Asignación de las condiciones iniciales
# Sse establece una función de sombrero como tal
# u(.5<=x<=1 && .5<=y<=1 ) is 2

u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2

### Representar condiciones iniciales
fig = plt.figure(figsize=(11,7), dpi=100)
##Los parametros figsize pueden se usados para cambiar el tamño y resolución

ax = fig.gca(projection='3d')                      
X, Y = np.meshgrid(x,y)                            
surf = ax.plot_surface(X,Y,u[:])
Populating the interactive namespace from numpy and matplotlib

Notas de representación de gráficos en 3D

Para representar un resultado en 3D, asegúrate de que has agregado 'Axes3D' proveniente de la librería matplotlib.

from mpl_toolkits.mplot3d import Axes3D

Los comandos para representar 3D son un poco más complicados que con 2d.

fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X,Y,u[:])

Con la primera línea se está creando una ventana con la figura. Los figsize y dpi son comandos opcionales y simplemente especifican el tamaño y la resolución de la figura que se reproduce. Se puede omitir, pero todavía se necesita el

fig = plt.figure()

La siguiente línea asigna a la figura la etiqueta ejes 'ax' y también le especifica que será una gráfica de proyección en 3D. La última línea utiliza el comando

plot_surface()

que es equivalente al comando normal de representar figuras, pero se necesita una malla de valores de X e Y para las posiciones de los puntos de datos.

Nota

Los valores de X e Y que representas con plot_surface no son vectores de x e y en 1D. Para usar las funciones 3D de matplotlib, necesitas crear una malla de valores x, y que corresponden con cada punto de la base. Las coordenadas de esta malla se generan usando la función de numpy meshgrid.

X, Y = np.meshgrid(x,y)

Iterando en dos dimensiones

Para evaluar la onda en dos dimensiones se necesita el uso de varios bucles for anidados para cubrir todos los indices i y j. Dado que Python no es un lenguaje compilado, puede haber notables ralentizamientos en la ejecución de código con múltiples bucles for. Primero, trata de evaluar el código de convección en 2D y observa qué resultados produce.

In [2]:
u = np.ones((ny,nx))
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2

for n in range(nt+1): ## Bucle para los incrementos de tiempo
    un[:] = u[:]
    for i in range(1, len(u)):
        for j in range(1, len(u)):
            u[i,j] = un[i, j] - (c*dt/dx*(un[i,j] - un[i-1,j]))-(c*dt/dy*(un[i,j]-un[i,j-1]))
            u[0,:] = 1
            u[-1,:] = 1
            u[:,0] = 1
            u[:,-1] = 1

fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X,Y,u[:]) 

Operaciones vectorizadas (Array Operations)

Aquí el mismo código de convección en 2D se ha implementado, pero en lugar de usar bucles for anidados, los mismos cálculos fueron evaluados usando operaciones vectorizadas.

In [16]:
u = np.ones((ny,nx))
u[.5/dy:1/dy+1,.5/dx:1/dx+1]=2

for n in range(nt+1): ## Bucle para los incrementos de tiempo
    un[:] = u[:]
    u[1:,1:]=un[1:,1:]-(c*dt/dx*(un[1:,1:]-un[0:-1,1:]))-(c*dt/dy*(un[1:,1:]-un[1:,0:-1]))
    u[0,:] = 1
    u[-1,:] = 1
    u[:,0] = 1
    u[:,-1] = 1

fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
surf2 = ax.plot_surface(X,Y,u[:])
plt.show()

Aprende más

El siguiente video te enseñará un desarrollo detallado de este Paso 5 (en dirección hasta el Paso 8). Video Lesson 6 de YouTube:

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

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