Tutorial: Numpy simultions of 2D gas

Tutorial: Numpy simultions of 2D gas#

import numpy as np
import matplotlib.pyplot as plt
from ipywidgets import interactive

import matplotlib.animation as animation
from IPython.display import HTML

class GasSimulation:

    def __init__(self, N=100, L=10.0, dt=0.01, radius=0.2, temperature=1.0, mass=1.0, steps=500):
        """Initialize gas particle simulation."""

        self.N, self.L, self.dt, self.radius, self.temperature, self.mass, self.steps = N, L, dt, radius, temperature, mass, steps
        
        # Initialize positions (left side) and velocities (Maxwell-Boltzmann distribution)
        self.positions  = np.column_stack((np.random.uniform(0, L/2, N), np.random.uniform(0, L, N)))
        self.velocities = np.random.normal(0, np.sqrt(temperature / mass), (N, 2))

        # Pre-allocate arrays for performance
        self.positions_history  = np.zeros((steps, N, 2), dtype=np.float32)
        self.velocities_history = np.zeros((steps, N, 2), dtype=np.float32)

    def update_positions(self):
        """Update positions and handle wall collisions."""

        self.positions += self.velocities * self.dt
        collisions = (self.positions - self.radius < 0) | (self.positions + self.radius > self.L) # Grab particles that collided with the the walls
        self.velocities[collisions] *= -1                                                         # Reverse the velocity upon collision
        self.positions = np.clip(self.positions, self.radius, self.L - self.radius)               # Clip positions and make sure no particle is exceeding the limits

    def run_simulation(self):
        """Run simulation and store trajectory history."""

        for step in range(self.steps):
            self.update_positions()
            self.positions_history[step], self.velocities_history[step] = self.positions.copy(), self.velocities.copy()

# Create and run simulation
sim = GasSimulation(N=100, steps=500)
sim.run_simulation()

print(sim.positions_history.shape) # Needs some tools to summarize the information

(500, 100, 2)

Post-processing simulation#

# Create and run simulation
sim = GasSimulation(N=100, steps=500)
sim.run_simulation()

visualize(sim)

animate(sim)

Project 1: Plot velocity and particle distributions#

#plot_speed_distribution(sim.velocities_history)
#plot_position_heatmap(sim.positions_history)

Project 2: Implement elastic collisions between gas particles#

1. Conservation Laws#

When two particles collide elastically, both momentum and kinetic energy are conserved.

Momentum Conservation#

m_{1} v_{1, init} + m_{2} v_{2, init} = m_{1} v_{1, final} + m_{2} v_{2, final}

Kinetic Energy Conservation#

\frac{1}{2} m_{1} v_{1, init}^{2} + \frac{1}{2} m_{2} v_{2, init}^{2} = \frac{1}{2} m_{1} v_{1, final}^{2} + \frac{1}{2} m_{2} v_{2, final}^{2}

where:

$v_{1, init}$ , $v_{2, init}$ are the initial speeds of particles 1 and 2.
$v_{1, final}$ , $v_{2, final}$ are the final speeds after the collision.

2. Collision Setup#

For two colliding particles, we define:

Positions: $r_{1}, r_{2}$
Velocities (before collision): $v_{1, init}, v_{2, init}$
Masses: $m_{1}, m_{2}$

The relative velocity between the two particles before the collision is:

v_{rel} = v_{1, init} - v_{2, init}

The unit vector along the collision axis (the direction along the line joining the particle centers) is:

\hat{r} = \frac{r_{1} - r_{2}}{| r_{1} - r_{2} |}

3. Velocities After Collision#

After the collision, the velocity component along the collision axis changes, while the perpendicular components remain unchanged.

Using momentum and energy conservation, the updated velocities are:

v_{1, final} = v_{1, init} - \frac{2 m_{2}}{m_{1} + m_{2}} (v_{rel} \cdot \hat{r}) \hat{r}

v_{2, final} = v_{2, init} + \frac{2 m_{1}}{m_{1} + m_{2}} (v_{rel} \cdot \hat{r}) \hat{r}