12 pasos para Navier-Stokes¶
¿Hicistes cambios en los Pasos 1 y 2 usando diferentes valores de parámetros? Si es así, es probable que te encontraras con un comportamiento inesperado. ¿Tu solución nunca se descontroló? (En mi experiencia, los estudiantes de CFD aman hacer que las cosas fallen).
Probablemente te estés preguntando porqué el cambio de los parámetros de discretización afecta a la solución de una manera tan drástica. Este notebook complementa nuestras lecciones interactivas de CFD al discutir la condición CFL. Y aprender más observando las clases en YouTube de la profesora Barba (enlaces más abajo).
Convergencia y la condición CFL¶
Para los primeros pasos, hemos estado usando las mismas condiciones iniciales y de contorno generales. Con los parámetros que inicialmente sugerimos, la cuadrícula cuenta con 41 puntos y el incremento de tiempo es de 0,025 segundos. Ahora, vamos a experimentar con el aumento del tamaño de nuestra cuadrícula. El siguiente código es idéntico al código que se utilizó en el Paso 1 pero aquí se ha incluido en una función para que podamos examinar fácilmente lo que sucede a medida que se ajusta una sola variable: el tamaño de la cuadrícula.
%pylab inline
# El comando de arriba hará que figuras de este notebook se representen junto al texto
import numpy as np # numpy es una libreria que realiza operaciones matriciales estilo MATLAB
import matplotlib.pyplot as plt # matplotlib es una librería para dibujar gráficas en 2D
# Listos para crear la función linearconv()
def linearconv(nx):
dx = 2./(nx-1)
nt = 20 # nt es el número de intervalos de tiempo que se desea calcular
dt = .025 # dt es la cantidad de tiempo que cada incremento de tiempo comprende (delta t)
c = 1
u = np.ones(nx) # definiendo 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) #nicializar el marcador de posición de conjunto de 'un', para almacenar la solución de tiempo instante por instante (n+1)
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):
u[i] = un[i]-c*dt/dx*(un[i]-un[i-1])
plt.plot(np.linspace(0,2,nx),u)
Populating the interactive namespace from numpy and matplotlib
Ahora vamos a examinar los resultados de nuestro problema de convección lineal con una malla (cuadrícula) mas fina cada vez.
linearconv(41) # convecction usando 41 puntos en la malla (mesh)
Este es el mismo resultado que en nuestro de cálculo del Paso 1, reproducido arriba para servir como referencia.
linearconv(61)
Aquí, aún está presente cierta difusión numérica (numerical diffusion), pero es menos grave.
linearconv(71)
Aquí se representa el mismo patrón, la onda es más cuadradada que en previos resultados.
linearconv(85)
Esto no se parece en nada a nuestra función sombrero original.
¿Qué ha ocurrido?¶
Para responder a esta pregunta, tenemos que pensar un poco acerca de lo que en realidad estamos implementando en el código.
En cada iteración de nuestro bucle que incrementa el tiempo, se utilizan los datos existentes acerca de nuestra onda para estimar la velocidad de la onda en el intervalo de tiempo siguiente. Inicialmente, el aumento en el número de puntos de la cuadrícula nos devolvió respuestas más precisas. Hubo menos difusión numérica y la onda cuadrada se veía mucho más como tal con respecto a lo que lo hizo en el primer ejemplo.
Cada iteración de nuestro bucle de tiempo abarca un incremento de longitud $\Delta t$, que hemos estado definiendo como 0.025
Durante esta iteración, se evalúa la velocidad de la onda en cada uno de los puntos $x$ que hemos creado. En la última ocasión, algo ha ido claramente mal.
Lo que ha pasado es que a lo largo del período de tiempo $\Delta t$, la onda se desplaza una distancia que es mayor que dx. La longitud dx de cada celda de la malla (cuadrícula) está relacionado con el número de puntos totales nx, por lo que la estabilidad puede asegurarse si el incremento de tiempo $\Delta t$ se calcula con respecto al tamaño de dx.
donde $u$ es la velocidad de la onda, $\sigma$ se conoce como el número de Courant y el valor de $\sigma_{max}$ que garantice la estabilidad depende de la discretización utilizada.
En una nueva versión de nuestro código, vamos a utilizar el número CFL para calcular un incremento de tiempo dt apropiado en función del tamaño de dx.
import numpy as np
import matplotlib.pyplot as plt
def linearconv(nx):
dx = 2./(nx-1)
nt = 20 # nt es el número de intervalos de tiempo que se desea calcular
c = 1
sigma = .5
dt = sigma*dx
u = np.ones(nx)
u[.5/dx : 1/dx+1]=2
un = np.ones(nx)
for n in range(nt): # teración a través del tiempo
un[:] = u[:] ## copia los valores existentes de 'u' en 'un'
for i in range(1,nx):
u[i] = un[i]-c*dt/dx*(un[i]-un[i-1])
plt.plot(np.linspace(0,2,nx),u)
linearconv(41)
linearconv(61)
linearconv(81)
linearconv(101)
linearconv(121)
Date cuenta que como el número de puntos, nx, aumenta, la onda experimenta la convección cada vez de forma más corta (menor "tamaño"). Si bien el número de iteraciones de tiempo que hemos hecho avanzar la solución se mantiene constante a nt = 20, dependiendo del valor de nx y los valores correspondientes de dx y dt se examina un intervalo de tiempo más corto en general, lo que provoca el efecto mencionado de acortamiento.
Aprende más¶
Es posible, en algunos casos, hacer un análisis riguroso de la estabilidad de los sistemas numéricos. Échale un vistazo a la presentación en inglés de la Profesora Barba sobre este tema en YouTube.
from IPython.display import YouTubeVideo
YouTubeVideo('Yw1YPBupZxU')
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.)