Sistema de discos rígidos. Teoría Cinética. Evolución hacia el equilibrio, reversibilidad e irreversibilidad.

• Teoría: Mecánica clásica, Teoría Cinética, Mecánica Estadística.
• Técnica numérica: Simulaciones de N-cuerpos. Colisiones elásticas. Algoritmos de Ordenación.
• Programación: Fortran, C, introducción a X-windows, gnuplot.

Modelo:

Se trata de estudiar el comportamiento de un sistema de N discos rígidos de masa unidad colocados en una mesa cuadrada de lado L. Los discos de radio a se mueven rectilineamente hasta que colisionan elásticamente entre si y/o contra las paredes que les rodean. Dadas unas posiciones y velocidades iniciales, r_i y v_i i = 1,..,N respectivamente, la evolución del sistema es determinista y reversible temporalmente. En particular la energía cinética total es una magnitud conservada en el tiempo.

Problema:

• Obligatorio: Diseñar y escribir el programa tal que dada cualquier condición inicial se pueda conocer la posición y velocidad de cada partícula después de un tiempo t.

• Voluntario 1: Sabemos (desde Clausius) que la temperatura del sistema es proporcional a la energía cinética total de las partículas dividida por el número de partículas. Puesto que ésta es una constante del movimiento en nuestro caso, dada una condición inicial cualquiera, estamos fijando la temperatura del sistema. La presión de las partículas sobre las paredes del recipiente, P, se calcula obteniendo la variación total del momento lineal de las partículas (la componente perpendicular a la pared) cuando colisionan con las paredes, Dp, durante un intervalo temporal Dt. Asi la presión viene dada por P = Dp/4LDt . Se pregunta:

  Dada una temperatura inicial, dejar que el sistema evolucione para que relaje a su estado de equilibrio termodinámico. Calcular la presión. Realizar el mismo experimento para distintas temperaturas iniciales. Obtener la relación entre presión y temperatura.

• Voluntario 2: Supongamos que inicialmente todas las partículas están confinadas en una de las esquinas de la caja. Dejar evolucionar al sistema un tiempo t. Invertir las velocidades de las partículas. Dejar evolucionar al sistema un tiempo t (igual al anterior). Puesto que el sistema es reversible temporalmente las partículas deberán volver a su estado inicial (partículas concentradas en una esquina). Se pregunta:
  • (1) Para que intervalo temporal el sistema es numéricamente reversible?.
  • (2) Elegir un tiempo t para el que el sistema es reversible. En el momento que invertimos las velocidades, elegir una partícula y cambiarla ligeramente de posición. Comentar lo que ocurre si dejamos evolucionar un tiempo t al sistema.

• Voluntario 3: Una vez que el sistema se encuentra en equilibrio, calcular la distribución de velocidades para varias temperaturas i.e. el número de partículas con un módulo de la velocidad entre v y v+dv. Obtener la distribución para varias temperaturas. Comparar los resultados numéricos con la predicción teórica.

Teoría:

Hemos de hallar las condiciones para que dos partículas colisionen y, si colisionan, obtener el cambio de velocidades correspondiente a un choque elástico. Igualmente hemos de hallar las condiciones de colisión de una partícula con las paredes y su cambio de velocidad si colisionan. Veamos primero las condiciones de colisión entre dos partículas.

Sean dos partículas con posiciones y velocidades (r₁,v₁) y (r₂,v₂) respectivamente. Al cabo de un tiempo t las partículas se encontrarán en

Las condiciones de colisión con las paredes son mucho mas simples de hallar. Nuestro sistema se compone de cuatro paredes: A: y = 0, B: x = L, C: y = L y D: x = 0. Sea una partícula con posición (x,y) y velocidad (v_x,v_y). Entonces las condiciones de colisión son:

• Pared A: Colisionará si se acerca a la pared, i.e. v_y < 0. El punto de contacto con la pared será x_col = x+v_x(a- y)/v_y y la colisión ocurrirá en t_col = (a-y)/v_y. Notar que la pared A se extiende desde x = 0 hasta x = L por lo tanto, en nuestro caso solo se hace efectiva la colisión si a < x_col < L-a.
• Pared B: Colisionará si v_x > 0 y si a < y_col < L-a donde y_col = y+v_y(L-a-x)/v_x. En ese caso el tiempo para que se produzca la colisión es: t_col = (L-a-x)/v_x.
• Pared C: Colisionará si v_y > 0 y si a < x_col < L-a donde x_col = x+v_x(L-a-y)/v_y. En ese caso el tiempo para que se produzca la colisión es: t_col = (L-a-y)/v_y.
• Pared D: Colisionará si v_x < 0 y si a < y_col < L-a donde y_col = y+v_y(a-x)/v_x. En ese caso el tiempo para que se produzca la colisión es: t_col = (a-x)/v_x.

Una vez que calculamos todos los tiempos de las colisiones posibles, se ha de elegir el menor de ellos y se ha de mover el sistema ese tiempo hasta que se realiza la colisión. La colisión que se produce puede ser entre partículas o entre una partícula y una pared. En cualquier caso, las colisiones son elásticas de forma que la energía cinética y el vector momento lineal antes y después de la colisión se conservan. Veamos como se modifican las velocidades de la/s partícula/s envuelta/s en cada caso.

• Colisión entre partículas: Supongamos que las partículas que están colisionando tiene posiciones r₁, r₂ y velocidades v₁, v₂. Sea r^pa = r₂-r₁ = (x,y) el vector que conecta los centros de los discos colisionando (x²+y² = 4a²). Sea r^pe = (-y,x) un vector perpendicular a r^pa. Las reglas de colisión son:

  v¢1pe = v1pe    v¢2pe = v2pe     v¢1pa = v2pa     v¢2pa = v1pa

  (7)

  donde v^pe es la proyección de v sobre la dirección dada por el vector r^pe, i.e. v^pe = v·r^pe/4a². Igualmente v^pa el la correspondiente proyección del vector sobre la dirección r^pa, i.e. v^pa = v·r^pa/4a². Los argumentos sin primas son las velocidades antes de la colisión y los que tienen primas corresponden a velocidades después de la misma. Las anteriores reglas de colisión se pueden expresar en coordenadas cartesianas:

  v¢1x = v1x+ x 4a2 [x(v2x-v1x)+y(v2y-v1y) ]

  (8)

  v¢1y = v1y+ y 4a2 [x(v2x-v1x)+y(v2y-v1y) ]

  (9)

  v¢2x = v2x- x 4a2 [x(v2x-v1x)+y(v2y-v1y) ]

  (10)

  v¢2y = v2y- y 4a2 [x(v2x-v1x)+y(v2y-v1y) ]

  (11)

• Colisión de una partícula contra la pared A (y = 0):v¢^x = v^x, v¢^y = -v^y.
• Colisión de una partícula contra la pared B (x = L):v¢^x = -v^x, v¢^y = v^y.
• Colisión de una partícula contra la pared C (y = L):v¢^x = v^x, v¢^y = -v^y.
• Colisión de una partícula contra la pared D (x = 0):v¢^x = -v^x, v¢^y = v^y.

Algoritmo numérico:

En la anterior sección teórica prácticamente hemos definido el algoritmo de evolución pero precisémoslo más.

• (0) Fijar un conjunto de posiciones y velocidades iniciales para las N partículas (NOTA: comprobar que los discos no solapen inicialmente). Las posiciones pueden variar x Î [0,L] y Î [0,L] y las velocidades pueden ser elegidas como se quieran, por ejemplo que tengan módulo unidad y un ángulo elegido al azar entre 0 y 2p.
• (1) Calcular de todas las posibles parejas de colisiones, la que tardará en producirse antes:

  • t_abs = 10⁹, n₁ = 0, n₂ = 0
  • do i=1,N
  • do j=i+1,N+4
  • Si no colisiona la partícula i con la j (j = N+1,N+2,N+3,N+4 son la paredes) ir a (a)
  • calcular t_col. Si t_col < t_abs entonces n₁ = i, n₂ = j y t_abs = t_col.
  • (a)enddo
  • enddo

• (2) Mover todas las partículas con sus respectivas velocidades un tiempo t_abs.
• (3) Aplicar la regla de colisión a la pareja n₁, n₂. Ir a (1)