De Maxwell-Boltzmann snelheidsverdeling

Doelen¶

We hebben al gezien in het vorige werkblad dat de snelheid van de deeltjes en de temperatuur aan elkaar gerelateerd zijn. In dit werkblad gaan we dat verder bekijken.

Eerst nemen we alle delen over van de code die we opnieuw moeten gebruiken:

class voor particles
functies voor een lijst aan deeltjes

Daarna voegen we de code toe voor het bekijken van de dynamiek van de deeltjes:

We kijken naar de verschillende richtingen
We bepalen de kansverdeling van de snelheden van de deeltjes

Zie voor de verdere inhoudelijke achtergrond de Feynman lectures on Physics.

Laden van eerdere code¶

Eerst nemen we de pakketten van Python en de constanten voor de simulatie over. Voer je eigen getallen in, die je ook in de vorige werkbladen hebt gebruikt.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.optimize import curve_fit

BOX_SIZE_0 = 10 #m           # Hoogte en breedte startvolume (in vul je gekozen eenheid in)
N = 40                          # Aantal deeltjes 
V_0 = 1 #m/s                  # Startsnelheid van deeltjes (vul je gekozen eenheid in)
RADIUS = 0.3 #m              # Straal van moleculen (vul je gekozen eenheid in)
DT = 0.1 * RADIUS / V_0 #s                   # Tijdstap om geen botsing te missen

#your code/answer
K_b = 1.38e-23 # J/K

Dan maken we weer gebruik van de klasse voor het deeltje:

class ParticleClass:
    def __init__(self, m, v, r, R):
        """ maakt een deeltje (constructor) """
        self.m = m                         
        self.v = np.array(v, dtype=float)  
        self.r = np.array(r, dtype=float)  
        self.R = R

    def update_position(self):
        """ verandert positie voor één tijdstap """
        self.r += self.v * DT 
            
    @property
    def momentum(self):
        return self.m * self.v
    
    @property
    def kin_energy(self):
        return 1/2 * self.m * np.dot(self.v, self.v)
    
def collide_detection(p1: ParticleClass, p2: ParticleClass) -> bool:
    """ Geeft TRUE als de deeltjes overlappen """
    return np.linalg.norm(p1.r - p2.r) < (p1.R + p2.R)


def particle_collision(p1: ParticleClass, p2: ParticleClass):
    """ past snelheden aan uitgaande van overlap """
    m1, m2 = p1.m, p2.m
    delta_r = p1.r - p2.r
    delta_v = p1.v - p2.v
    dot_product = np.dot(delta_r, delta_v)
    # Als deeltjes van elkaar weg bewegen dan geen botsing
    if dot_product >= 0: # '='-teken voorkomt ook problemen als delta_r == \vec{0}
        return
    distance_squared = np.dot(delta_r, delta_r) 
    # Botsing oplossen volgens elastische botsing in 2D
    p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
    p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r

En we laten de deeltjes met de wanden botsen, zodat er sprake is van een druk en een temperatuur in een gesloten volume.

def box_collision(particle: ParticleClass):
    ''' botsing met harde wanden '''
    if abs(particle.r[0]) + particle.R > BOX_SIZE_0/2: 
        particle.v[0] = -particle.v[0]                                        # Omdraaien van de snelheid
        particle.r[0] = np.sign(particle.r[0]) * (BOX_SIZE_0/2 - particle.R)  # Zet terug net binnen box                 
    if abs(particle.r[1]) + particle.R > BOX_SIZE_0/2: 
        particle.v[1] = -particle.v[1]     
        particle.r[1] = np.sign(particle.r[1]) * (BOX_SIZE_0/2 - particle.R)

Functies aan een lijst van deeltjes¶

Waarbij we al deze functies uitvoeren en samennemen over een lijst met deeltjes:

def create_particles(particles):
    """ Leegmaken en opnieuw aanmaken van deeltjes  in lijst """
    particles.clear()
    for i in range(N):
        vx = np.random.uniform(-V_0, V_0)
        vy = np.random.choice([-1, 1]) * np.sqrt(V_0**2 - vx**2)        
        pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2)
        particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=RADIUS))
        
def handle_collisions(particles):
    """ alle onderlinge botsingen afhandelen voor deeltjes in lijst """
    num_particles = len(particles)
    for i in range(num_particles):
        for j in range(i+1, num_particles):
            if collide_detection(particles[i], particles[j]):
                particle_collision(particles[i], particles[j])

def handle_walls(particles):
    """ botsing met wanden controleren voor alle deeltjes in lijst """
    for p in particles:
        box_collision(p)

def take_time_step(particles):
    """ zet tijdstap voor een lijst deeltjes en verwerk alle botsingen onderling en met wanden """
    for p in particles:
        p.update_position()
    handle_collisions(particles)
    handle_walls(particles)

Eerste simulatie ter controle¶

Zoals we inmiddels gewend zijn draaien we eerst een korte simulatie om te verifiëren of alle code correct is overgenomen:

particles = []
create_particles(particles)
for i in range(100):
    take_time_step(particles)

plt.figure()
plt.xlabel('x')
plt.ylabel('y')
plt.gca().set_aspect('equal')
plt.xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
plt.ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)

for p in particles:
    plt.plot(p.r[0], p.r[1], 'k.', ms=10)
    plt.arrow(p.r[0], p.r[1], p.v[0], p.v[1], 
              head_width=0.05, head_length=0.1, color='red')
plt.show()

Equipartitiebeginsel¶

Dit beginsel is heel belangrijk voor de thermodynamica en stelt dat de energie in thermodynamisch evenwicht gelijk wordt verdeeld over de toegankelijke vrijheidsgraden. Laten we dit eerst verifiëren in onze simulatie.

# jouw antwoord
def create_uniform_particles(particles):
    """ Leegmaken en opnieuw aanmaken van deeltjes met uniforme snelheid in lijst """
    particles.clear()
#your code/answer
for i in range(N):
        vx = 0.0
        vy = np.random.choice([-1, 1]) * V_0 
        pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2)
        particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=RADIUS))

class ParticleClass:
    def __init__(self, m, v, r, R):
        self.m = m
        self.v = np.array(v, dtype=float)
        self.r = np.array(r, dtype=float)
        self.R = R

    def update_position(self):
        self.r += self.v * DT

def collide_detection(p1, p2):
    return np.linalg.norm(p1.r - p2.r) < (p1.R + p2.R)

def particle_collision(p1, p2):
    m1, m2 = p1.m, p2.m
    delta_r = p1.r - p2.r
    delta_v = p1.v - p2.v
    dot_product = np.dot(delta_r, delta_v)
    if dot_product >= 0: return
    distance_squared = np.dot(delta_r, delta_r)
    p1.v -= 2 * m2 / (m1 + m2) * dot_product / distance_squared * delta_r
    p2.v += 2 * m1 / (m1 + m2) * dot_product / distance_squared * delta_r

def box_collision(particle):
    if abs(particle.r[0]) + particle.R > BOX_SIZE_0/2:
        particle.v[0] = -particle.v[0]
        particle.r[0] = np.sign(particle.r[0]) * (BOX_SIZE_0/2 - particle.R)
    if abs(particle.r[1]) + particle.R > BOX_SIZE_0/2:
        particle.v[1] = -particle.v[1]
        particle.r[1] = np.sign(particle.r[1]) * (BOX_SIZE_0/2 - particle.R)

def handle_collisions(particles):
    num_particles = len(particles)
    for i in range(num_particles):
        for j in range(i+1, num_particles):
            if collide_detection(particles[i], particles[j]):
                particle_collision(particles[i], particles[j])

def handle_walls(particles):
    for p in particles:
        box_collision(p)

def take_time_step(particles):
    for p in particles:
        p.update_position()
    handle_collisions(particles)
    handle_walls(particles)

# Functie om deeltjes te maken met snelheid alleen in y-richting
def create_uniform_particles(particles):
    particles.clear()
    for i in range(N):
        vx = 0.0
        vy = np.random.choice([-1, 1]) * V_0
        # Willekeurige posities
        pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2)
        particles.append(ParticleClass(m=1.0, v=[vx, vy], r=pos, R=RADIUS))

# Simulatie uitvoeren
particles = []
create_uniform_particles(particles)

# Bereken aantal stappen
total_distance = 2 * BOX_SIZE_0
total_time = total_distance / V_0
num_steps = int(total_time / DT)

# Data opslaan voor plot
time_values = []
E_kin_x = []
E_kin_y = []

for step in range(num_steps):
    take_time_step(particles)
    
    # Kinetische energie in x- en y-richting berekenen
    Ex = sum(0.5 * p.m * p.v[0]**2 for p in particles)
    Ey = sum(0.5 * p.m * p.v[1]**2 for p in particles)
    
    time_values.append(step * DT)
    E_kin_x.append(Ex)
    E_kin_y.append(Ey)

# Plotten
plt.figure()
plt.plot(time_values, E_kin_x, label='$E_{kin,x}$')
plt.plot(time_values, E_kin_y, label='$E_{kin,y}$')
plt.xlabel('Tijd (s)')
plt.ylabel('Kinetische Energie (J)')
plt.title('Kinetische energie in x- en y-richting als functie van tijd')
plt.legend()
plt.show()

Aan het resultaat van je simulatie kan je zien dat de deeltjes zich inderdaad gedragen volgens equipartitiebeginsel. Afwijkingen van het gemiddelde zijn het gevolg van de statistiek en worden voor grotere aantallen deeltjes relatief kleiner.

De snelheidsverdeling¶

In het vorige werkblad hebben we al gezien dat er een relatie is tussen de snelheid van de deeltjes en de temperatuur. Je kan hierdoor al vermoeden dat de snelheid van de deeltjes wordt gegeven door een verdeling die van de temperatuur afhangt. De vraag is dus of we deze functie kunnen vinden.

Om de snelheidsverdeling van de deeltjes te bepalen kan je gebruik maken van de functie histogram van python. Laten we daarom een simulatie draaien waarin we vanuit een willekeurige beginsituatie 100 tijdstappen zetten en daarna de snelheidsverdeling plotten:

particles = []

create_uniform_particles(particles)
for i in range(num_steps):
    take_time_step(particles)

speeds = [np.linalg.norm(p.v) for p in particles]
counts, bins = np.histogram(speeds, bins=10, density='True')

plt.figure()
plt.xlabel('Speed')
plt.ylabel('Count')
plt.stairs(counts, bins, fill='True')
plt.show()

Voldoende statistiek¶

Als je de simulatie hierboven een aantal keer uitvoert, dan zal je zien dat de statistiek onvoldoende is om een reproduceerbaar antwoord te krijgen. We kunnen het aantal deeltjes toe laten nemen om meer statistiek te krijgen, maar dat kost heel veel rekenkracht. Het is een goedkopere oplossing om de statistiek te bepalen op verschillende momenten in de tijd en deze statistische resultaten te middelen. De onderstaand code laat de deeltjes 5000 tijdstappen zetten en noteert de snelheid van alle deeltjes op elke 250e tijdstap. Het histogram wordt dan bepaald door de snelheden bij elke 250e tijdstap samen te nemen.

Door de simulatie een aantal keer te draaien zie je dat de verdeling al een stuk stabieler wordt.

particles = []
num_bins = 10
tot_counts = np.zeros(num_bins, dtype=float)

plt.figure()
plt.xlabel('Speed')
plt.ylabel('Count')

create_uniform_particles(particles)

for i in range(5000):
    take_time_step(particles)
    if i % 250 == 0:
        speeds = [np.linalg.norm(p.v) for p in particles]
        partial_counts, bins = np.histogram(speeds, bins=num_bins, density='True', range=[0,3*V_0])
        tot_counts += partial_counts

norm_counts = tot_counts / 20
plt.stairs(norm_counts, bins, fill='True')
plt.show()

De mathematische vorm van de snelheidsverdeling¶

De meest algemene vorm van de snelheidsverdeling heeft de vorm $f_{2D}(\vec{v})$ . Dat wil zeggen dat er bij elke unieke combinatie van $x$ en $y$ -component van de snelheid een specifieke kans hoort. We weten echter al dat dit niet het geval kan zijn. De natuur heeft helemaal geen voorkeur voor richting en wij kunnen zelf kiezen hoe ons assenstelsel georiënteerd is. De kans is dus alleen afhankelijk van de modulus van de snelheid en onafhankelijk van de richting. We kunnen de verdeling daarom weergeven als $f_{2D}(v)$ (zonder pijl voor de vector).

We kunnen de verdeling nog scherper definiëren door te stellen dat de componenten van de snelheid onderling onafhankelijk zijn. De functie $f_{2D}$ is daarom te splitsen in aparte functies voor de $x$ en $y$ -richting. Combineren we dit met onze conclusie van de vorige paragraaf dan moet dus gelden dat we de functie kunnen splitsen in twee functies voor $x$ en $y$ die onderling precies hetzelfde zijn en ook hetzelfde als $f(v)$ :

f_{2D}(v)=f(v_x)f(v_y)

(1)

Laten we nu de snelheidsverdeling beschouwen langs een willekeurige richting $r$ , die een lineaire combinatie van de $x$ - en $y$ -richting is. Omdat dit een deelverzameling is van de twee-dimensionale snelheidsverdeling moet bovenstaande relatie hiervoor nog steeds gelden. Om de vorm van deze functie te vinden stellen we nu een differentiaalvergelijking op door te differentiëren naar $v_x$ :

\frac{d}{dv_x}f(v_r)=\frac{d}{dv_x}f(v_x)f(v_y)

(2)

Om de linkerkant uit te werken, passen we de kettingregel toe:

\left(\frac{\partial f(v)}{\partial v_r}\right)\left(\frac{\partial v_r}{\partial v_x}\right)=f(v_y)\frac{df(v_x)}{dv_x}

(3)

Met de Stelling van Pythagoras wordt dit:

\frac{v_x}{v_r}\left(\frac{\partial f(v_r)}{\partial v_r}\right)=f(v_y) \frac{df(v_x)}{dv_x}

(4)

Om deze differentiaalvergelijking op te lossen willen we een scheiding van variabelen toepassen. Daarvoor willen we eerst af van de variabele $v_y$ . Dat kan door de vergelijking te delen door $f(v_r)=f(v_x)f(v_y)$ en de $v_x$ naar de andere kant te brengen:

\frac{1}{v_r f(v_r)}\left(\frac{\partial f(v_r)}{\partial v_r}\right)=\frac{1}{v_x f(v_x)} \frac{df(v_x)}{dv_x}

(5)

Je kan precies dezelfde redenering ook opzetten vanuit de $y$ -coördinaat in plaats van de $x$ -coördinaat. De termen die je hier gevonden hebt zijn daarmee functies van verschillende en onderling onafhankelijke variabelen die toch hetzelfde antwoord geven. Ze moeten daarom constant zijn. Die constante noemen we $-2\alpha$ , omdat dit de formules verderop vereenvoudigt:

\frac{1}{v_r f(v_r)}\left(\frac{\partial f(v_r)}{\partial v_r}\right)=\frac{1}{v_x f(v_x)} \frac{df(v_x)}{dv_x}=\frac{1}{v_y f(v_y)} \frac{df(v_y)}{dv_y} = -2\alpha

(6)

Misschien herken je hier al de vorm die $f(v_x)$ moet hebben om hieraan te voldoen, maar om dit makkelijker te maken, kunnen we de vergelijking een beetje herschrijven:

\frac{df(v_x)}{f(v_x)}=-2\alpha v_x dv_x

(7)

We kunnen beide zijden nu rustig integreren, zodat de oplossing wordt gegeven door:

f(v_x)=A \exp\left(-\alpha v_x^2\right)

(8)

Bij het vak Multivariabele Analyse zal je deze integraal tegenkomen aan het einde van dit blok. Daar wordt bewezen dat de integraal onder deze functie precies “1” is als (we zeggen ook wel: de functie is genormaliseerd):

A = \sqrt{\frac{\alpha}{\pi}}

(9)

Om ook de waarde van $\alpha$ te bepalen, kunnen we nu eisen dat deze formule de juiste waarde moet geven voor het gemiddelde van het kwadraat van de snelheid in de $x$ -richting:

\left< v_x^2 \right> = \frac{kT}{m} = \sqrt{\frac{\alpha}{\pi}} \int_{-\infty}^{\infty} v_x^2 \exp\left(-\alpha v_x^2\right)dv_x

(10)

Partieel integreren levert op:

\alpha = \frac{m}{2kT},

(11)

zodat de snelheidsverdeling voor twee dimensies uiteindelijk de vorm heeft:

f_{2D}(v_x,v_y)=f(v_x)f(v_y)=\frac{m}{2\pi kT} \exp\left( -\frac{m(v_x^2+v_y^2)}{2kT} \right)

(12)

Om hieruit de snelheidsverdeling $f_{2D}(v)$ te bepalen, moeten we nog een extra stap nemen. Het aantal combinaties van $v_x$ en $v_y$ dat overeenkomt met de snelheid $v$ is namelijk afhankelijk van $v$ . Dit wordt gegeven door de cirkelomtrek met middelpunt $(v_x=0,v_y=0)$ en straal $v$ . Zodoende is de snelheidsverdeling voor de modus van de snelheid gegeven door:

f_{2D}(v)=\frac{mv}{kT} \exp\left( -\frac{mv^2}{2kT} \right)

(13)

particles = []
create_particles(particles)

all_speeds = []

# Parameters
total_simulation_steps = 3000
equilibration_steps = 500  
sample_interval = 20       

# Loop
for step in range(total_simulation_steps):
    take_time_step(particles)
    
    if step > equilibration_steps and step % sample_interval == 0:
        
        current_speeds = [np.linalg.norm(p.v) for p in particles]
        all_speeds.extend(current_speeds)

# Histogram maken
counts, bin_edges = np.histogram(all_speeds, bins=25, density=True)
bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2

# De fit functie
def mb_2d_dist(v, a):
    return a * v * np.exp(-a * v**2 / 2)

# Fitten
popt, pcov = curve_fit(mb_2d_dist, bin_centers, counts, p0=[1.0])
a_fit = popt[0]

# Temperatuur berekenen
m_particle = 1.0
T_calculated = m_particle / a_fit

print(f"Gefitte parameter a: {a_fit:.3f}")
print(f"Berekende Temperatuur (in simulatie-eenheden met k_B=1): {T_calculated:.3f}")

# Lijn plot
v_axis = np.linspace(0, max(all_speeds), 100)
fit_curve = mb_2d_dist(v_axis, a_fit)

# Plotten
plt.figure(figsize=(10, 6))

# Histogram van de data
plt.hist(all_speeds, bins=25, density=True, alpha=0.6, color='skyblue', label='Simulatie data', edgecolor='black')

# De fit lijn
plt.plot(v_axis, fit_curve, 'r-', linewidth=2.5, label=f'Maxwell-Boltzmann Fit (T={T_calculated:.2f})')

plt.title('2D Maxwell-Boltzmann Snelheidsverdeling')
plt.xlabel('Snelheid $v$ (m/s)')
plt.ylabel('Kansdichtheid')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Gefitte parameter a: 1.902
Berekende Temperatuur (in simulatie-eenheden met k_B=1): 0.526

# Definitie constanten
m = 1.0       # Massa van een deeltje
V_0 = 1.0     # Initiële snelheid
k_B_sim = 1.0 # Boltzmann constante in simulatie-eenheden

T_theorie = (0.5 * m * V_0**2) / k_B_sim

# Vergelijk met de gefitte waarde uit Opdracht 3
print(f"Ingestelde snelheid V_0:      {V_0:.2f}")
print(f"Theoretische Temperatuur:     {T_theorie:.3f}")
print(f"Gefitte Temperatuur (uit Q3): {T_calculated:.3f}")

# Bereken de afwijking
afwijking = abs(T_theorie - T_calculated) / T_theorie * 100
print(f"Afwijking:                    {afwijking:.2f}%")

Ingestelde snelheid V_0:      1.00
Theoretische Temperatuur:     0.500
Gefitte Temperatuur (uit Q3): 0.526
Afwijking:                    5.16%