← Back
Computational Physics Course Project (Group)

2D N-Body Thermodynamics Simulation

AdvisorNTHU Computational Physics
PeriodSpring 2024

A group project for the NTHU Computational Physics course, simulating a 2D N-body system to verify fundamental thermodynamic principles from a microscopic perspective. The simulation encompasses three physical scenarios: ideal gas kinetic theory, adiabatic compression, and the mixing of two gas populations toward thermal equilibrium. My personal contribution formed the foundational architecture of the entire codebase — I wrote the first working version of the simulation engine, including the spatial grid partitioning (grid_N.py), all boundary reflection logic, Numba-accelerated parallel kernels, and the inter-particle viscosity force model.

  • ▪Personal Contribution — Spatial Grid Partitioning (grid_N): Designed and implemented grid_N(), which divides the 2D simulation box into N² sub-cells and builds a neighbor-grid lookup dictionary. This reduces collision-checking complexity from O(N²) to O(N × neighbors), enabling simulations of up to 10,000 particles.
  • ▪Personal Contribution — Boundary Reflection & Pressure Engine: Built the complete bound_reflection() logic handling all four walls with configurable moving velocities, elastic velocity inversion, pressure accumulation via momentum transfer per unit boundary length, and position correction to keep particles inside the domain.
  • ▪Personal Contribution — Numba Parallelization: Accelerated the inter-particle force (viscosity) and elastic collision kernels using Numba @njit(parallel=True) with prange loops for multi-threaded execution.
  • ▪Personal Contribution — Viscosity (Inter-particle Force): Implemented an exponential repulsive force model between nearby particles: F = b₁ · exp(−b₀ · r) · r̂, enabling simulation of non-ideal gas behavior and particle interaction dynamics.
  • ▪Personal Contribution — First Working Simulation: Wrote and validated the initial end-to-end version of the code, confirming that the complete architecture was physically correct and usable as the group's foundational base.
  • ▪Group Work — RK2 Integration & Analysis: The team integrated RK2 time-stepping for high accuracy, and contributed to the theoretical analysis of adiabatic compression and pressure-volume behavior.

Methodology

1. Spatial Grid Partitioning (grid_N)

To overcome the O(N²) bottleneck of all-pairs collision detection, the simulation domain is partitioned into an N×N grid of sub-cells. Each particle is assigned to a cell at every time step via the partition() function. During force and collision evaluation, only particles within the same cell or its 8 neighboring cells (as encoded by the grid_N() lookup dictionary) are considered as potential interaction pairs. This reduces the effective interaction search space by a factor of N²/9, providing a substantial speedup for large N.

2. Boundary Reflection & Pressure Calculation

The bound_reflection() function handles all four boundaries (top, bottom, left, right) generically. Each wall can be assigned an independent constant velocity (bound_vel), enabling moving-wall adiabatic compression experiments. When a particle crosses a boundary:

- Its velocity component normal to the wall is reversed and corrected for the wall's own motion: v⊥′=−v⊥+2vwallv'_\perp = -v_\perp + 2 v_{\text{wall}}

- Its position is reflected back inside the domain

- The instantaneous pressure contribution is accumulated as:

Pi=1Δt⋅L∣∑collidingΔp⊥∣P_i = \frac{1}{\Delta t \cdot L} \left| \sum_{\text{colliding}} \Delta p_\perp \right|

where LL is the length of that boundary and Δp⊥\Delta p_\perp is the normal momentum change per particle.

3. Inter-particle Viscosity Force

A soft exponential repulsive force is applied between nearby particles to model viscosity effects:

Fij=b1⋅e−b0rij⋅r^ij\mathbf{F}_{ij} = b_1 \cdot e^{-b_0 r_{ij}} \cdot \hat{\mathbf{r}}_{ij}

where b0b_0 controls the range and b1b_1 the strength of the repulsion. This force is computed only between particles in neighboring grid cells (leveraging grid_N), and the inner loop is parallelized with Numba's @njit(parallel=True) for high throughput.

4. RK2 Time Integration

Particle trajectories are advanced using a second-order Runge-Kutta (RK2) scheme:

vn+1=vn+a1+a22Δt\mathbf{v}_{n+1} = \mathbf{v}_n + \frac{\mathbf{a}_1 + \mathbf{a}_2}{2} \Delta t
xn+1=xn+v1+v22Δt\mathbf{x}_{n+1} = \mathbf{x}_n + \frac{\mathbf{v}_1 + \mathbf{v}_2}{2} \Delta t

where a1\mathbf{a}_1, a2\mathbf{a}_2 are accelerations evaluated at the current and predicted-next positions respectively.

5. Elastic Inter-Particle Collisions

A collision is triggered when ∣rij∣<2r|\mathbf{r}_{ij}| < 2r and Δv⋅rij<0\Delta\mathbf{v} \cdot \mathbf{r}_{ij} < 0. Velocities are exchanged via the center-of-mass frame:

vi′=2vcm−vi,vj′=2vcm−vj\mathbf{v}'_i = 2\mathbf{v}_{\text{cm}} - \mathbf{v}_i, \quad \mathbf{v}'_j = 2\mathbf{v}_{\text{cm}} - \mathbf{v}_j

ensuring exact conservation of both momentum and kinetic energy.

Results & Analysis

Particle Trajectory Comparison

The trajectory plots (0–25 particles) clearly show the qualitative effect of each physical mechanism: straight-line bouncing for ideal gas, chaotic scattering with inter-particle collisions, smooth oscillatory paths with viscosity only, and highly complex trajectories when both collision and viscosity are active simultaneously.

Adiabatic Compression

Pressure curves matched theoretical adiabatic predictions most closely when the boundary velocity was low (≈ 2.5 code units), maintaining a quasi-static condition. At higher boundary velocities (v=15), the system deviated significantly from the theoretical curve — consistent with non-equilibrium thermodynamics breaking the quasi-static assumption.

Mixing of Two Gas Populations

Two groups of particles initialized with different mean velocities were allowed to interact. The total kinetic energy remained approximately constant (confirming elastic collision energy conservation), while the average velocities of both groups converged to a common equilibrium value within ~15 time units, demonstrating the equipartition theorem from statistical mechanics.

Figures

Figure 1: Particle trajectories (0–25) with both inter-particle collisions and viscosity — showing full non-ideal gas complexity

Figure 1: Particle trajectories (0–25) with both inter-particle collisions and viscosity — showing full non-ideal gas complexity

Figure 2: Total kinetic energy vs time — near-constant curve confirms elastic collision energy conservation

Figure 2: Total kinetic energy vs time — near-constant curve confirms elastic collision energy conservation

Figure 3: Average velocity of two particle groups (v1, v2) converging to thermal equilibrium within ~15 time units

Figure 3: Average velocity of two particle groups (v1, v2) converging to thermal equilibrium within ~15 time units

Figure 4: Adiabatic compression pressure vs time (boundary velocity = 2.5) — simulation closely tracks theoretical prediction

Figure 4: Adiabatic compression pressure vs time (boundary velocity = 2.5) — simulation closely tracks theoretical prediction

Simulation Videos

Video 1: Mixing of two particle groups with different initial temperatures — velocity distribution converges to thermal equilibrium

Video 2: Adiabatic compression with a moving boundary (v=20) — particles accelerate as the boundary compresses the domain