A gas of hard spheres in 3D, checked against kinetic theory. Originally written for the PHYS389 computer modelling module and extended since.
200 particles bounce around a cubic box, colliding elastically with the walls and with each other. They all start with exactly the same speed (the rms speed for 293 K) in random directions, so the speed histogram begins as a single spike. After a few collisions each it has relaxed onto the Maxwell-Boltzmann distribution, which is the orange curve in the animation. One particle is coloured red so you can follow it.
From that one ingredient, elastic collisions between spheres, the simulation reproduces:
- the Maxwell-Boltzmann distribution, with no fitted parameters
- the pressure and collision rate of a gas of finite-sized particles
- the Carnahan-Starling equation of state up to a packing fraction of 0.21
- Boltzmann's H-theorem
- equipartition between light and heavy particles
- Brownian motion of a large sphere, from ballistic to diffusive
pip install -e .[dev]
python scripts/animate.py # animation
python scripts/animate.py --save # writes figures/animation.gif instead
python scripts/animate.py -N 500 --radius 0.02 --event-driven
python scripts/analyse.py # prints pressure and collision rate, saves a figure
pytest # tests
The other scripts each produce one of the figures below: equation_of_state.py, h_theorem.py, mixture.py and brownian.py. The first and last take a few minutes.
To use it in your own code:
from hard_sphere_gas import Simulation, EventDrivenSimulation
sim = EventDrivenSimulation(N=300, R=0.02, T=300, seed=1)
sim.run(1000)
print(sim.temperature(), sim.pressure(), sim.collision_rate())m and R can be arrays with one value per particle, and pos and vel can be passed in to set the starting state by hand. The state is in the numpy arrays sim.pos, sim.vel, sim.mass and sim.radius.
hard_sphere_gas/simulation.py- theSimulationclass, which uses a fixed timestep. Each step moves every particle byv*dt, then checks for wall hits and overlapping pairs.hard_sphere_gas/event_driven.py-EventDrivenSimulation, which has the same interface but works out exactly when each collision will happen and jumps from one to the next.hard_sphere_gas/theory.py- the kinetic theory results the simulation is compared with.scripts/- the animation and one script per figure.tests/- pytest tests for both engines and the theory functions.
Collisions are resolved by applying an impulse along the line between the two centres, which conserves both momentum and kinetic energy. In the fixed timestep engine a pair is only collided if the particles are approaching each other, and a wall only reflects a particle that is heading outwards. Without those two checks particles can get stuck to each other or to the walls, flipping direction every step.
Finding the touching pairs. For up to 200 particles the distances between all pairs are found at once with numpy. Above that a cell list is used: the box is divided into cells at least one particle diameter wide, so a particle can only be touching particles in its own cell or the 26 around it. This takes the cost per step from O(N²) to O(N), and 2000 particles run at 5 ms per step instead of 190 ms.
The event-driven engine. Between collisions the particles move in straight lines, so the time at which any two will touch is the root of a quadratic. Each particle's earliest event (a wall or another particle) is kept in a priority queue. The earliest event overall is popped, all the particles are moved forward to that time, the collision is carried out and the particles involved have their next events predicted again. Predictions are made stale by collisions that happen in the meantime, so each event records how many collisions its particles had had when it was made and is thrown away if that has changed. Spheres never overlap and there is no timestep to choose, but it is slower in pure Python, at around 3000 collisions a second.
The mass is 127 u (iodine atoms) and the box is 1 m across. The radius of 3 cm is obviously nothing like a real atom, it is chosen so that the mean free path is smaller than the box and collisions happen often enough to watch.
Left: speeds sampled after the first 1000 steps against the Maxwell-Boltzmann distribution at the temperature of the simulation. There are no fitted parameters. Right: relative change in total kinetic energy, which stays at the level of floating point rounding.
Pressure is measured by adding up the momentum given to the walls (2m|v| per hit) and dividing by time and wall area. Here it is as a ratio to the ideal gas value NkT/V, where V is the volume available to the particle centres:
| P / Pideal | Collisions per particle per second | |
|---|---|---|
Fixed timestep, dt = 5e-5 |
1.062 | 753 |
| Event-driven | 1.105 | 829 |
| Theory for an infinite system | 1.117 | 912 |
The theory values are the Carnahan-Starling equation of state and the Enskog collision rate, sqrt(2) n pi d^2 <v> g, where g is the pair distribution function at contact. (The van der Waals excluded volume correction NkT/(V - Nb) gives 1.122, the two agree at low density.) The measured mean free path is 0.27 m against 0.24 m.
Two separate effects pull the measurements below the theory:
- Timestep. With a fixed timestep two spheres overlap a little before the collision is noticed, so they behave as if they were smaller than they are. This accounts for the gap between the first two rows, and the event-driven engine removes it completely.
- Walls. The remaining gap is real physics rather than a numerical error. A particle next to a wall cannot be hit from the wall side, so it collides less often than one in the middle of the box. With 6 cm particles in a 1 m box a good fraction of the gas is within one diameter of a wall. Both the excess pressure and the collision rate come out at about 90% of the infinite-system value. Running
EventDrivenSimulation(N=1600, R=0.015), which has almost the same packing fraction but particles relatively half as close to the walls, brings both up to 95%, as expected for a surface effect. The next section gets rid of it properly.
The compressibility factor Z = PV/NkT for 250 particles at six packing fractions, using the event-driven engine. The open circles use all the wall hits and the overall density, and fall further below the Carnahan-Starling curve as the particles get bigger, because of the wall effect above.
The filled circles use only the middle of the box. The pressure on a flat wall is equal to the pressure of the bulk gas, so it is measured from hits on the central part of each wall, at least two particle diameters away from the edges. The density of the bulk gas is measured directly by counting the particles in the middle of the box, and is lower than N/V since particles collect next to the walls. Done this way the simulation agrees with Carnahan-Starling to within 1.6% at every density:
| Packing fraction (bulk) | Z measured | Carnahan-Starling |
|---|---|---|
| 0.020 | 1.084 | 1.083 |
| 0.049 | 1.227 | 1.220 |
| 0.094 | 1.484 | 1.480 |
| 0.136 | 1.761 | 1.783 |
| 0.176 | 2.148 | 2.146 |
| 0.211 | 2.574 | 2.534 |
The van der Waals excluded volume correction is only right to first order in the density, and is already well off by a packing fraction of 0.1.
Boltzmann's H is the integral of f ln f over velocity, where f is the velocity distribution. It is minus the entropy of the distribution, and the H-theorem says collisions can only make it smaller until it reaches its minimum, which is the Maxwell-Boltzmann distribution. This is the relaxation in the animation at the top, as a single number. With 2000 particles H falls smoothly to the Maxwell-Boltzmann value (-19.05 for iodine at 293 K in SI units) within about four collisions per particle, and stays there.
150 xenon-mass particles start with all of the kinetic energy and 150 argon-mass particles start at rest. Collisions share the energy out until each species has the same mean kinetic energy, 3/2 kT, after which they swap small amounts back and forth. The two species end up with different Maxwell-Boltzmann speed distributions at the same temperature.
One sphere of radius 15 cm and 20 times the mass in a bath of 500 ordinary particles. Left: its path, a random walk. Right: its mean squared displacement. For short times the sphere just keeps going and the MSD grows as t^2. After about a millisecond the bath has randomised its velocity and the MSD grows as 6Dt instead, until the sphere starts to notice the walls of the box.
The curves are not fits. The friction on the sphere is calculated from Enskog kinetic theory, which gives the diffusion coefficient D = kT / friction and the velocity relaxation time M / friction, and these go into the solution of the Langevin equation (dotted) and the same for a particle confined between walls (solid).
- The box has hard walls rather than periodic boundaries, so bulk properties pick up a surface correction of a few percent, as discussed above.
- The event-driven engine predicts each particle's collisions against every other particle, which is O(N) work per collision. Combining it with the cell list would make this O(1).
- Everything is pure Python and numpy. Compiling the inner loops with numba would make both engines much faster.
MIT, see LICENSE.