Convección en una dimensión



In [1]:

    
%matplotlib inline



In [2]:

    
import scipy as np
from matplotlib import pyplot as plt

Parámetros



In [3]:

    
L  = 1.0           # longitud del sistema 1D
nx = 42            # nodos espaciales
dx = L / (nx-2)    # sí, quitamos dos nodos ...
x = np.linspace( 0 , L , num=nx )

T= 0.1            # tiempo total
nt = 100            # pasos temporales
dt = T / nt

c = 1              # velocidad de la onda

¡Número de Courant !



In [4]:

    
Co = c * dt / dx
Co









    Out[4]:





0.04

Condiciones iniciales



In [5]:

    
u0 = 1 * np.ones(nx)                    # todo uno
x1 = L/4 ; n1 = int(x1 / dx)
x2 = L/2 ; n2 = int(x2 / dx)
u0[ n1 : n2 ] = 2



In [6]:

    
plt.plot( x , u0 )









    Out[6]:





[<matplotlib.lines.Line2D at 0x7f643a62c190>]

Un paso en el tiempo

Recordemos que queremos implementar $u_i^{n+1} = u_i^n - \mathrm{Co}/2 (u_{i+1}^n-u_{i-1}^n)$



In [7]:

    
u = u0.copy()



In [8]:

    
un = u.copy()         # distribución actual

i = 1
u[i] = un[i] - (Co / 2.0) * (un[i+1] - _valor_izdo )



for i in range( 2 , nx - 2 ):    # Ahora queda claro por qué hemos quitado los extremos !!
   u[i] = un[i] - (Co / 2.0) * (un[i+1] - un[i-1])



In [9]:

    
plt.plot(x,u)









    Out[9]:





[<matplotlib.lines.Line2D at 0x7f643a591050>]

Tiempo completo



In [93]:

    
u = u0.copy()



In [94]:

    
for n in range(nt):
    un = u.copy()
    for i in range( 1 , nx - 1 ): 
        u[i] = un[i] - (Co / 2.0) * (un[i+1] - un[i-1])



In [95]:

    
plt.plot(x , u , x , u0 , 'r')









    Out[95]:





[<matplotlib.lines.Line2D at 0x7fb46945bb90>,
 <matplotlib.lines.Line2D at 0x7fb46945bd90>]

Formulaciones alternativas

¿Qué pasa si probamos el algoritmo "peor"? $u_i^{n+1} = u_i^n - \mathrm{Co} (u_{i}^n-u_{i-1}^n)$



In [96]:

    
u = u0.copy()



In [97]:

    
for n in range(nt):
    un = u.copy()
    for i in range( 1 , nx ): 
        u[i] = un[i] - Co  * (un[i] - un[i-1])



In [98]:

    
plt.plot(x , u , x , u0 , 'r')









    Out[98]:





[<matplotlib.lines.Line2D at 0x7fb469386790>,
 <matplotlib.lines.Line2D at 0x7fb469386990>]

Reflexiones

¿Qué hemos hecho en los bordes? (condiciones de contorno)
¿Por qué los algoritmos difieren tantísimo?
¿Cuántos parámetros hay realmente? (Pista: 1)