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,\text{init}} + m_2 v_{2,\text{init}} = m_1 v_{1,\text{final}} + m_2 v_{2,\text{final}} \]

Kinetic Energy Conservation#

\[ \frac{1}{2} m_1 v_{1,\text{init}}^2 + \frac{1}{2} m_2 v_{2,\text{init}}^2 = \frac{1}{2} m_1 v_{1,\text{final}}^2 + \frac{1}{2} m_2 v_{2,\text{final}}^2 \]

where:

\(v_{1,\text{init}}\), \(v_{2,\text{init}}\) are the initial speeds of particles 1 and 2.
\(v_{1,\text{final}}\), \(v_{2,\text{final}}\) are the final speeds after the collision.

2. Collision Setup#

For two colliding particles, we define:

Positions: \( \mathbf{r}_1, \mathbf{r}_2 \)
Velocities (before collision): \( \mathbf{v}_{1,\text{init}}, \mathbf{v}_{2,\text{init}} \)
Masses: \( m_1, m_2 \)

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

\[ \mathbf{v}_{\text{rel}} = \mathbf{v}_{1,\text{init}} - \mathbf{v}_{2,\text{init}} \]

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

\[ \hat{\mathbf{r}} = \frac{\mathbf{r}_1 - \mathbf{r}_2}{|\mathbf{r}_1 - \mathbf{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:

\[ \mathbf{v}_{1,\text{final}} = \mathbf{v}_{1,\text{init}} - \frac{2 m_2}{m_1 + m_2} (\mathbf{v}_{\text{rel}} \cdot \hat{\mathbf{r}}) \hat{\mathbf{r}} \]

\[ \mathbf{v}_{2,\text{final}} = \mathbf{v}_{2,\text{init}} + \frac{2 m_1}{m_1 + m_2} (\mathbf{v}_{\text{rel}} \cdot \hat{\mathbf{r}}) \hat{\mathbf{r}} \]