Distribución de Energía Cinética de Bolas en una Caja 2D

Me preguntaba si sería posible simular la distribución de energía cinética de un montón de bolas que se están moviendo en una caja.

Así que me puse a la tarea e implementé estos pequeños programas en Python, como sencillos ficheros de texto que se ejecutan desde un terminal en Linux.

Las bolas son similares a las bolas de billar, porque chocan entre sí y contra las paredes, pero no tienen ningún otro efecto adicional. No hay perdidas de energía, y el número de bolas es constante durante toda la simulación.

Hay tres programas que hay que ejecutar:

  • 1createPositions.py, que crea un fichero data.dat con la información sobre las bolas.
  • 2makePlotsAll.py, que lee el fichero anterior data.dat y crea una gráfica para cada instante de tiempo, dentro del directorio tmp.
  • 3createAnimation.sh, que junta todas las gráficas anteriores para crear un pequeño vídeo animation.mp4.

Los tres programas se ejecutan uno detrás del otro, de este modo:

./1createPositions.py
./2makePlotsAll.py
./3createAnimation.sh

Para hacer que estos programas sean ejecutables, tendrá el lector que utilizar los comandos (en Linux):

chmod u+x 1createPositions.py
chmod u+x 2makePlotsAll.py
chmod u+x 3createAnimation.sh

El primer programa crea un fichero data.dat con toda la informacion de las bolas en cada instante de tiempo.

Las condiciones iniciales son bastante excepcionales. Las bolas se encuentran distribuidas aleatoriamente en posicion, pero su velocidad es muy particular. Solo hay una componente horizontal, hacia la derecha, igual para todas las bolas. De esta forma, todas las bolas tienen la misma energia al principio, y no chocan entre si hasta que impactan con la pared derecha.

El segundo programa de Python crea gráficas como la siguiente. El panel superior muestra la posición de las bolas en la caja, el panel inferior la distribución de energía cinética en ese instante de tiempo (1/2*masa*velocidad², pero como todas las bolas son iguales, ignorando el 1/2 y la masa, quedándonos sólo con el cuadrado de la velocidad de la bola).

figure_00298

1createPositions.py


#!/usr/bin/env python

###############################################################################
# Initialisation
###############################################################################

# Imports libraries
import random
import math

###############################################################################
# Constants
###############################################################################

# Number of steps
Nsteps=300

# Number of balls
Nballs=1500

# Size of balls
r=0.002

# Size of field (X and Y)
x0=0.0
x1=1.0
y0=0.0
y1=1.0

# Initial velocity (the same for all balls)
v0=0.005

###############################################################################
# Displays general information on screen
###############################################################################

print
print 'This scripts creates a set of testing data'
print 'Number of time-steps:',Nsteps
print 'Number of balls:',Nballs
print

###############################################################################
# Initialise data
###############################################################################

x=[0.0]*Nballs
y=[0.0]*Nballs
vx=[0.0]*Nballs
vy=[0.0]*Nballs

for nb in range(Nballs):
    x[nb]=x0+(x1-x0)*random.random()
    y[nb]=y0+(y1-y0)*random.random()
    ang=0.0
    vx[nb]=v0*math.cos(ang)
    vy[nb]=v0*math.sin(ang)

###############################################################################
# Open output file
###############################################################################

file=open('data.dat','w')

# Write header
file.write('#Time nBall      xPos      yPos       xVel       yVel\n')

###############################################################################
# Main part
###############################################################################

for ns in range(Nsteps):
    print 'Time-step',ns,'out of',Nsteps-1
    for nb in range(Nballs):
        x[nb]+=vx[nb]
        y[nb]+=vy[nb]
        # Bounce on other particles
        for nbb in range(nb):
            dist=math.sqrt((x[nb]-x[nbb])**2.0+(y[nb]-y[nbb])**2.0)+2e-9
            if dist<=(2.0*r): # Bounce # Second ball is quiet, the first impacts on it xCM=(x[nb]+x[nbb])/2.0 yCM=(y[nb]+y[nbb])/2.0 vxRef=vx[nbb] vyRef=vy[nbb] ang=math.atan2(y[nb]-y[nbb],x[nb]-x[nbb]) # radial and tangential components vr=+(vx[nb]-vxRef)*math.cos(ang)+(vy[nb]-vyRef)*math.sin(ang) vt=-(vx[nb]-vxRef)*math.sin(ang)+(vy[nb]-vyRef)*math.cos(ang) vx[nb]=-vt*math.sin(ang)+vxRef vy[nb]=+vt*math.cos(ang)+vyRef vx[nbb]=vr*math.cos(ang)+vxRef vy[nbb]=vr*math.sin(ang)+vyRef # Make sure that both balls are not inside each other x[nb]=xCM+r*math.cos(ang) y[nb]=yCM+r*math.sin(ang) x[nbb]=xCM-r*math.cos(ang) y[nbb]=yCM-r*math.sin(ang) # Bounce on box limits if x[nb]>(x1-r):
            x[nb]=x1-r
            vx[nb]=-abs(vx[nb])
        if x[nb]<(x0+r): x[nb]=x0+r vx[nb]=+abs(vx[nb]) if y[nb]>(y1-r):
            y[nb]=y1-r
            vy[nb]=-abs(vy[nb])
        if y[nb]<(y0+r):
            y[nb]=y0+r
            vy[nb]=+abs(vy[nb])
        # Save to file
        line='{0:5d} {1:5d} {2:.7f} {3:.7f} {4:+.7f} {5:+.7f}\n'.format(ns,nb,x[nb],y[nb],vx[nb],vy[nb])
        file.write(line)

###############################################################################
# Close output file
###############################################################################

file.close()

print
print 'File data.dat is created'
print

###############################################################################
# The End
###############################################################################

2makePlotsAll.py


#!/usr/bin/env python

# Import libraries
import matplotlib.pyplot as plt
import numpy as np
import gc
import pylab as pl
import os

# Displays some information on screen
print
print 'This script reads the file data.dat and creates'
print 'plots in folder ./tmp'
print

# Load data
matrix=np.loadtxt('data.dat') # Load a huge matrix containing all data
time=matrix[:,0]  # Time-step, first=0
nBall=matrix[:,1] # Number of ball, first=0
x=matrix[:,2]     # Position in X axis
y=matrix[:,3]     # Position in Y axis
vx=matrix[:,4]    # Velocity in X axis
vy=matrix[:,5]    # Velocity in Y axis

# Free memory
del matrix
gc.collect()

# Number of time-steps
Nsteps=int(max(time)-min(time))

# Number of balls
Nballs=int(max(nBall)-min(nBall)+1)

# Define energy array (kinetic energy = 1/2*m*v^2)
vxSub=np.array([0.0]*Nballs)
vySub=np.array([0.0]*Nballs)
v2=np.array([0.0]*Nballs)

# Creates temporal folder
if not os.path.exists('./tmp'):
    os.makedirs('./tmp')
    print './tmp is created'

# Prepare plot
fig = plt.figure()

# Main part
for nStep in range(Nsteps):
    fileOutput='tmp/figure_{0:05d}.png'.format(nStep)
    print fileOutput
    n0=nStep*Nballs
    # Calculate energy and normalise
    vxSub=vx[n0:n0+Nballs]    # Subset of the vx array for this time-step
    vySub=vy[n0:n0+Nballs]    # Subset of tye vy array for this time-step
    v2=vxSub*vxSub+vySub*vySub
    v2=v2*1e6        # [micro-units]
    # Make plot
    ax = fig.add_subplot(211)
    ax.plot( x[n0:n0+Nballs], y[n0:n0+Nballs], 'o' )
    plt.xlim([0.0,1.0])
    plt.ylim([0.0,1.0])
    # Make histogram
    ax = fig.add_subplot(212)
    ax.hist(v2,50,range=[0.0,140.0],normed=0)
    ax.set_xlabel('Average Energy, v^2 [micro-units]')
    ax.set_ylabel('Number of Cases')
    # Save figure to file
    pl.savefig(fileOutput)
    # Clear plot
    plt.clf()

# The End

3createAnimation.sh

#!/bin/sh

echo
echo &amp;quot;This script combines all plots from ./tmp and&amp;quot;
echo &amp;quot;creates the file animation.mp4&amp;quot;
echo

# 24 frames/second
avconv -r 24 -i tmp/figure_%5d.png -s 1600x1200 animation.mp4

Autor: willyfog

Turista laboral por la Unión Europea. Por favor que dure. Lo que veo, leo o me cuentan no lo suelo encontrar en español, así que me gusta escribirlo por aquí.

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión /  Cambiar )

Google photo

Estás comentando usando tu cuenta de Google. Cerrar sesión /  Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión /  Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión /  Cambiar )

Conectando a %s