A gas in a box: pressure from collisions

Script by Ordinate 2D ✓ · Rotation and rolling

Forty atoms bounce in a closed box; the wall's pressure against N⟨KE⟩/A.

Use in Ordinate 2DOpens in the Script tab: it runs when you press Run. Free, no account needed.

"""A gas in a box: pressure from collisions, and P·A = N·⟨KE⟩.

Sixty small disks bounce elastically in a closed box with no gravity. Each time one hits the right wall, the momentum
it gives the wall is added up; the total over time is a force, and force per wall length is the 2D pressure. For an ideal
gas in two dimensions P·A = N·k·T = N·⟨kinetic energy⟩ — the script writes the pressure so far next to it every 2 s.

Runs in the studio (Script tab → Run) and locally:  python3 gas_in_box.py
"""
import math
import random
import studio as s

N, R, MASS = 60, 0.07, 0.05
W, H = 6.0, 4.0                                  # the box, m
rng = random.Random(1)

s.clear()
s.gravity = (0, 0)
walls = {
    "left": s.box(-W / 2 - 0.1, H / 2, 0.2, H + 0.4, static=True, name="Left wall", restitution=1, friction=0, color="#8b939b"),
    "right": s.box(W / 2 + 0.1, H / 2, 0.2, H + 0.4, static=True, name="Right wall", restitution=1, friction=0, color="#e5484d"),
    "floor": s.box(0, -0.1, W + 0.4, 0.2, static=True, name="Floor", restitution=1, friction=0, color="#8b939b"),
    "roof": s.box(0, H + 0.1, W + 0.4, 0.2, static=True, name="Roof", restitution=1, friction=0, color="#8b939b"),
}
gas = []
for i in range(N):
    a = rng.uniform(0, 2 * math.pi)
    speed = rng.gauss(3.0, 0.8)
    gas.append(s.circle(rng.uniform(-W / 2 + 0.3, W / 2 - 0.3), rng.uniform(0.3, H - 0.3), R, mass=MASS, restitution=1, friction=0,
                        name=f"Atom {i + 1}", color="#4f8ff7", velocity=(speed * math.cos(a), speed * math.sin(a))))
last_vx = [b.vx for b in gas]
impulse = 0.0


@s.on_step
def measure(t, dt):
    global impulse
    for i, b in enumerate(gas):
        if b.x > W / 2 - R - 0.05 and last_vx[i] > 0 > b.vx:     # just bounced off the right wall
            impulse += MASS * (last_vx[i] - b.vx)
        last_vx[i] = b.vx
    if round(t / dt) % 480 == 0 and t > 0:           # every 2 s
        pressure = impulse / t / H
        ke = sum(0.5 * MASS * b.speed ** 2 for b in gas) / N
        s.record(t=t, pressure=pressure, ideal_gas=N * ke / (W * H), mean_KE=ke)
        print(f"t = {t:4.1f} s: P = {pressure:.3f} N/m, N⟨KE⟩/A = {N * ke / (W * H):.3f} N/m")


s.run(20)

table collisions thermal