Por un momento, recordemos las ecuaciones de Navier-Stokes para un fluido incompresible, donde $\vec{v}$ representa el campo de velocidades:
\begin{eqnarray} \nabla \cdot\vec{v} &=& 0\\\ \frac{\partial \vec{v}}{\partial t}+(\vec{v}\cdot\nabla)\vec{v} &=& -\frac{1}{\rho}\nabla p + \nu \nabla^2\vec{v} \end{eqnarray}La primera ecuación representa la conservación de la masa manteniendo la densidad constante. La segunda ecuación es la conservación del momento. Sin embargo, aparece un problema: la ecuación de continuidad para el flujo incompresible no tiene una variable dominante y no hay manera obvia para acoplar la velocidad y la presión. En el caso de flujo compresible, por contra, la continuidad de la masa proporcionaría una ecuación con la evolución de la densidad $\rho$, que está acoplada con una ecuación de estado que relaciona $\rho$ y $p$.
En flujo incompresible, la ecuación de continuidad $\nabla \cdot\vec{v}=0$ proporciona una restricción cinemática que requiere que el campo de presión evolucione de manera que la tasa de expansión $\nabla \cdot\vec{v}$ debería desaparecer en todas partes. Una forma de salir de esta dificultad es construir un campo de presión que garantice que la continuidad se satisface, tal relación se puede obtener mediante la adopción de la divergencia de la ecuación de momento. En este proceso, ¡una ecuación de Poisson para la presión aparece!
La ecuación de Poisson se obtiene añadiendo un término fuente (source term) a la derecha de la ecuación de Laplace:
$$\frac{\partial ^2 p}{\partial x^2} + \frac{\partial ^2 p}{\partial y^2} = b$$Así, a diferencia de la ecuación de Laplace, hay algo de valor finito en el interior del campo que afecta a la solución. La ecuación de Poisson actúa para "relajar" las fuentes iniciales en el campo.
En forma de datos discretos, esto se ve casi igual que en el Paso 9, a excepción del término de fuente:
$$\frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2 p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2}=b_{i,j}^{n}$$Al igual que antes, lo reorganizamos para obterner una ecuación para $p$ en el punto $i,j$. Por lo tanto, obtenemos:
$$p_{i,j}^{n}=\frac{(p_{i+1,j}^{n}+p_{i-1,j}^{n})\Delta y^2+(p_{i,j+1}^{n}+p_{i,j-1}^{n})\Delta x^2-b_{i,j}^{n}\Delta x^2\Delta y^2}{2(\Delta x^2+\Delta y^2)}$$Vamos a resolver esta ecuación, suponiendo un estado inicial de $p=0$ en todas partes, y aplicando las siguientes condiciones de contorno:
$p=0$ en $x=0, \ 2$ y $y=0, \ 1$
y el término fuente se compone de dos picos iniciales en el interior del dominio, de la siguiente manera:
$b_{i,j}=100$ en $i=nx/4, j=ny/4$
$b_{i,j}=-100$ en $i=nx*3/4, j=3/4 ny$
$b_{i,j}=0$ en el resto de puntos.
Las iteraciones avanzarán en pseudo-tiempo para relajar los picos iniciales. La relajación bajo la ecuación de Poisson se vuelve más y más lenta a medida que evolucionan. ¿Por qué?
Echemos un vistazo a una posible manera de escribir el código de la ecuación de Poisson. Como siempre, cargamos nuestras librerías Python favoritas. También queremos hacer algunas figuras preciosas en 3D. Obtengamos nuestros parámetros que hemos definidos y la inicialización. ¿Qué notas del planteamiento propuesto a continuación?
%pylab inline
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt
import numpy as np
# Parámetros
nx = 50
ny = 50
nt = 100
xmin = 0.
xmax = 2.
ymin = 0.
ymax = 1.
dx = (xmax-xmin)/(nx-1)
dy = (ymax-ymin)/(ny-1)
# Inicialización
p = np.zeros((nx,ny))
pd = np.zeros((nx,ny))
b = np.zeros((nx,ny))
x = np.linspace(xmin,xmax,nx)
y = np.linspace(xmin,xmax,ny)
# Fuente
b[nx/4][ny/4] = 100
b[3*nx/4][3*ny/4] = -100
Populating the interactive namespace from numpy and matplotlib
Con esto, estamos listos para avanzar en la estimación inicial de pseudo-tiempo. ¿Cómo es el código de abajo diferente de la función que se utiliza en el Paso 9 para resolver la ecuación de Laplace?
for it in range(nt):
pd[:][:]=p[:][:]
p[1:nx-1,1:ny-1] = ( dy**2/(2*(dx**2+dy**2))*(pd[2:nx,1:ny-1]+pd[0:nx-2,1:ny-1]) +
dx**2/(2*(dx**2+dy**2))*(pd[1:nx-1,2:ny]+pd[1:nx-1,0:ny-2]) -
b[1:nx-1,1:ny-1]*dx**2*dy**2/(2*(dx**2+dy**2)) )
p[0,:] = p[nx-1,:] = p[:,0] = p[:,ny-1] = 0.0
Tal vez podríamos reutilizar el la funcion de representar del Paso 9, ¿no te parece?
def plot2D(x, y, p):
fig = plt.figure(figsize=(11,7), dpi=100)
ax = fig.gca(projection='3d')
X,Y = np.meshgrid(x,y)
surf = ax.plot_surface( X,Y,p[:], rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False )
ax.set_xlim(0,2)
ax.set_ylim(0,1)
ax.view_init(30,225)
plot2D(x, y, p)
¡Ah! Las maravillas de la reutilización del código! Ahora, probablemente pienses: "Bueno, si he escrito esta pequeña función tan estupenda que hace algo tan útil, pero quiero usarla una y otra vez ¿Cómo puedo hacer esto sin tener que copiar y pegar todo el rato? -Si tienes mucha curiosidad acerca de esto, tendrás que aprender acerca de packaging. Pero esto va más allá del alcance de nuestras lecciones de CFD. Sólo tienes que buscarlo en Google si realmente quieres saber cómo ;)
Para obtener más información sobre la función de la ecuación de Poisson en CFD, puedes ver el Video Lesson 11 en YouTube:
from IPython.display import YouTubeVideo
YouTubeVideo('ZjfxA3qq2Lg')
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)