termalización

718 days ago by pablogc

Distribución Maxwelliana de velocidades

Vamos a tomar un gas que parte con una cierta distribución inicial de velocidades.

Las colisiones entre partículas las modelizaremos con un pequeño y sencillo montecarlo que hará que la distribución inicial de velocidades se transforme poco a poco en una distribución Maxwelliana.

Ingredientes de la simulación

No nos va a interesar simular las posiciones de las partículas sino únicamente sus velocidades.

Asimismo vamos a realizar toda la simulación en dos dimensiones.

Vamos a generar una distribución inicial de velocidades con módulo de velocidad fijo:

\|\vec v\|\equiv 1 y ángulo aleatorio \theta=\mbox{random}
def genera_1(particulas): atomos=np.zeros((particulas,2)) for i in srange(particulas): theta=random()*2*pi atomos[i][0]=cos(theta) atomos[i][1]=sin(theta) return atomos 
       

Una segunda posible configuración inicial tomará valores aleatorios entre 0 y 1 para cada una de las componentes de la velocidad:
def genera_2(particulas): atomos=np.zeros((particulas,2)) for i in srange(particulas): atomos[i][0]=random() atomos[i][1]=random() return atomos 
       

Una tercera tomará un módulo de la velocidad aleatorio entre 0 y 1 y también un ángulo aleatorio:
def genera_3(particulas): atomos=np.zeros((particulas,2)) for i in srange(particulas): theta=random()*2*pi v=random() atomos[i][0]=v*cos(theta) atomos[i][1]=v*sin(theta) return atomos 
       

Dibujando la distribución de velocidades

Para poder visualizar estas distribuciones vamos a histogramarlas en un cierto intervalo (min,max) y con nd bines. Para acelerar este proceso lo más posible vamos a utilizar arrays numéricos de las librerías numpy y compilar el código en C utilizando Cython. Porteriormente, nos limitaremos a mover átomos de un bin a otro, proceso que no requerirá ya de tanta velocidad. 

 
%cython import numpy as np cimport numpy as np cdef double min cdef double max cdef int nd min=-2.5 max=2.5 nd=100 def deltabin(): return 100/(max-min) def indice(double x): #'Toma los límites y número de canales para histogramar' #'y un cierto valor x retornando la posición que ocupará' #'ese valor en el array del histograma' if x >= max: return -1 elif x <= min: return 0 else: return int(0.5+nd*(x-min)/(max-min)) def histograma(np.ndarray[double, ndim=2] atomos): #'Toma el array con las velocidades de los átomos' #'Retorna histogramados el módulo de la velocidad' #'y la componente x' h=np.zeros((nd+1,3)) d=(max-min)/nd j=0 while j <= nd: h[j][0]=min+j*d j+=1 j=0 while j < len(atomos): h[indice(sqrt(atomos[j][0]**2+atomos[j][1]**2))][1]+=1 h[indice(atomos[j][0])][2]+=1 j+=1 return h 
       

Colisiones

Dado que no hay posiciones de los átomos, lo que vamos a hacer es tirar un montecarlo que seleccione dos átomos de manera aleatoria y los haga colisionar con un parámetro de impacto b también aleatorio.
def ColAB(atomo1,atomo2,h): #'Hace colisionar los átomos dados' #'con un parámetro de impacto aleatorio' #'y retorna las velocidades de salida' b=random() h[indice(atomo1[0])][2]-=1 h[indice(atomo2[0])][2]-=1 h[indice(sqrt(atomo1[0]**2+atomo1[1]**2))][1]-=1 h[indice(sqrt(atomo2[0]**2+atomo2[1]**2))][1]-=1 vxij=atomo1[0]-atomo2[0] vyij=atomo1[1]-atomo2[1] fact=b*vxij+(1.-b**2)**0.5*vyij dvx=-b*fact dvy=-(1.-b**2)**0.5*fact atomo1[0]+=dvx atomo1[1]+=dvy atomo2[0]-=dvx atomo2[1]-=dvy h[indice(atomo1[0])][2]+=1 h[indice(atomo2[0])][2]+=1 h[indice(sqrt(atomo1[0]**2+atomo1[1]**2))][1]+=1 h[indice(sqrt(atomo2[0]**2+atomo2[1]**2))][1]+=1 reset([dvx,dvy,fact,b,vxij,vyij]) return atomo1, atomo2, h def Colisiones(ncols, atomos,h): #'Hace ncols colisiones por partícula' #'retornando la nueva población' i=0 while i < ncols*len(atomos): i1=int(random()*len(atomos)) i2=int(random()*len(atomos)) while i2==i1: i2=int(random()*len(atomos)) atomos[i1],atomos[i2],h = ColAB(atomos[i1],atomos[i2], h) i+=1 reset([i1,i2]) return atomos,h 
       

Interfaz gráfica

A continuación vamos a desarrollar una interfaz sencilla que permita elegir los parámetros de la simulación y generar visualizaciones:
global atomos, h, lista, accion @interact def _(inicial=['v=1, theta=rand','vx=rand, vy=rand','v=rand, theta=rand'], particulas = slider(0,100000,1000,5000,label='Partículas'), colisiones = slider(0,20,1,1,label='Colisiones (por partícula)'), nfoto = slider(1,10,1,1,label='Fotogramas por Colisión'), output=['plot modulo', 'plot componente', 'anima modulo', 'anima componente'], accion=['Stop', 'Calcula']): if accion == 'Stop': html('Seleccione los parámetros de la simulación y pulse <strong>Calcular</strong>') elif accion == 'Calcula': if inicial == 'v=1, theta=rand': html('Condición inicial: $\|\overline v\|=1\;,\;\Theta=\mbox{random}$') generador=genera_1 teorico=plot(2*x*exp(-x**2),0,2.5, rgbcolor=(0,1,0)) elif inicial == 'vx=rand, vy=rand': html('Condición inicial: $v_x=\mbox{random}\; ,\; v_y=\mbox{random}$') generador=genera_2 teorico=plot(3*x*exp(-2/3*x**2),0,2.5, rgbcolor=(0,1,0)) elif inicial == 'v=rand, theta=rand': html('Condición inicial: $\|\overline v\|=\mbox{random}\;,\;\Theta=\mbox{random}$') generador=genera_3 teorico=plot(6*x*exp(-3*x**2),0,2.5, rgbcolor=(0,1,0)) atomos=generador(particulas) h=histograma(atomos) try: lista.append(h*1) except: lista=[] lista.append(h*1) atm=atomos*1 for k in srange(colisiones*nfoto): atomos,h=Colisiones(0.5/nfoto,atm,h) lista.append(h*1) if output == 'plot modulo': i=list_plot(zip(lista[0][:,0],lista[0][:,1]),'l',rgbcolor=(1,0.5,0.5)) i.show(dpi=75,xmin=0) html('Ploteando $\|\overline v\|$ al cabo de %s colisiones por partícula'%(colisiones)) o=list_plot(zip(lista[-1][:,0],lista[-1][:,1]*deltabin()/particulas),'l')+teorico o.show(dpi=75,xmin=0) elif output == 'plot componente': i=list_plot(zip(lista[0][:,0],lista[0][:,2]),'l',rgbcolor=(1,0.5,0.5)) i.show(dpi=75) html('Ploteando $v_x$ al cabo de %s colisiones por partícula'%(colisiones)) o=list_plot(zip(lista[-1][:,0],lista[-1][:,2]*deltabin()/particulas),'l') o.show(dpi=75) elif output == 'anima modulo': html('Animando $\|\overline v\|$ al cabo de %s colisiones por partícula con %s fotogramas por colision'%(colisiones,nfoto)) o = animate(list_plot(zip(lista[i][:,0],lista[i][:,1]*deltabin()/particulas),'l',xmin=0)+teorico for i in srange(len(lista))) o.save('output.gif', show_path=True) elif output == 'anima componente': html('Animando $v_x$ al cabo de %s colisiones por partícula con %s fotogramas por colision'%(colisiones,nfoto)) o= animate(list_plot(zip(lista[i][:,0],lista[i][:,2]*deltabin()/particulas),'l')+teorico for i in srange(len(lista))) o.save('output.gif', show_path=True) accion='Stop' lista=[] 
       

Click to the left again to hide and once more to show the dynamic interactive window