import numpy as np
import matplotlib.pyplot as pltRecap en uitbreiding¶
Je hebt al gezien hoe je een botsend deeltjesmodel maakt. Je startte toen met 1 of 2 deeltjes en hebt dit uitgebreid naar ~100 deeltjes. Daarbij was het logisch om gebruik te maken van classes.
Voor het maken van een simulatie van gasmoleculen in een gesloten volume zal je logischerwijs de volgende stappen doorlopen:
bepalen van initiële condities
definiëren van de deeltjes
definiëren van volume met randvoorwaarden
simulaties bestaande uit:
bepaling van positie van deeltjes
controleren op onderlinge botsingen
controleren op botsingen met wanden
aanpassen van eigenschappen zoals totaal volume of temperatuur
opslaan van gegevens
De initiële condities¶
In onze simulatie gaan we uit van een aantal aannames. Bij deze aannames wil je dat de simulatie een voorspellende waarde heeft voor een reëel systeem, maar ook dat de benodigde rekenkracht voor het uitvoeren van de simulatie op te brengen is. We kiezen daarom voor de volgende condities:
Ons model heeft een beperkt en constant aantal deeltjes.
We maken gebruik van een 2D simulatie.
De deeltjes voelen geen onderlinge kracht en ondergaan alleen elastische botsingen.
(tijdens simuleren) De tijdstap is voldoende klein. (d.w.z. de afgelegde weg is klein in vergelijking met de diameter van de deeltjes, zodat we geen botsingen missen)
Hieronder geven we een aantal constanten voor onze code. Om de simulatie straks te kunnen interpreteren, vergelijken we deze eerst met een realistische situatie:
De simulatie die we gaan maken bevat maar twee dimensies. We kunnen het oppervlak per deeltje in deze simulatie daarom kiezen als het zojuist berekende .
Voor de gemiddelde snelheid, , van een molecuul geldt dat
Hierbij is de massa van het deeltje, het aantal dimensies, de constante van Boltzmann en de temperatuur.
Onze simulatie modelleert dus een klein systeem zonder rekening te houden met de kracht tussen de luchtmoleculen. Je zult zien dat we toch al behoorlijk wat van de processor vragen om deze berekeningen te doen. Om een meer realistische modelsimulatie te maken is daarom een serieuze uitdaging, die veel wiskunde en programmeerervaring vraagt. Op dit ogenblik wordt dit bijvoorbeeld gevraagd bij een van de eerste opdrachten van het vak ‘Computational Physics’ van de master Applied Physics.
# Berekening van het gemiddelde volume per molecuul bij 10 bar
k_B = 1.380649e-23 # [J/K] Boltzmannconstante
T = 300 # [K] kamertemperatuur
P = 10 * 1e5 # [Pa] druk van 10 bar
#V_molecuul = (k_B * T) / P # [m^3] gemiddeld volume per molecuul
V_molecuul = 0.01
print(f"geschatte Gemiddeld volume per molecuul: {V_molecuul:.3e} m^3")
geschatte Gemiddeld volume per molecuul: 1.000e-02 m^3
#opdracht 2:
# Constanten
k_B = 1.380649e-23 # Boltzmann-constante (J/K)
T = 293 # Temperatuur (K)
m = 4.8e-26 # Massa van N2-molecuul (kg)
f = 2 # aantal dimensies (2D)
# Berekening <v^2> voor 2D
v_squared = (f * k_B * T) / m
# Gemiddelde snelheid
v = np.sqrt(v_squared)
print(f"Gemiddelde snelheid van een gasmolecuul in 2D bij 293 K: {v:.2f} m/s")
Gemiddelde snelheid van een gasmolecuul in 2D bij 293 K: 410.55 m/s
#4
BOX_SIZE_0 = 1.0 # 1 meter (of 1 eenheid)
RADIUS = 0.02 # 2 cm bijvoorbeeld
V_0 = 410 # snelheid 1 m/s
N = 40 # aantal deeltjes
# beweging per frame ~ 0.01 meter
# TIJDSTAP ZODAT v*DT << diameter
# Doel: verplaatsing = v * DT < (2R / 100)
#DT = (2 * RADIUS) / (100 * V_0) # [s] tijdstap klein genoeg voor botsingsdetectie
DT = 0.00001
print("Gekozen DT=",DT,"s")Gekozen DT= 1e-05 s
Particle class¶
Onze ParticleClass definieert van een deeltje de massa , snelheid , positie , en straal . De positie moet per tijdstap bepaald worden.
Nieuw aan ParticleClass zijn twee properties, nl. momentum en energy. Een property is eigenlijk een variabele die samenhangt met de waarde van andere variabelen binnen de klasse. Om de interne consistentie te bewaren moet je de juiste waarde van deze variabele dus telkens afleiden en dat doe je met een functie. Zo geeft kin_energy de kinetische energie terug van het deeltje door die af te leiden van de variabelen m en v.
Buiten ParticleClass voegen we twee functies toe, die je ook al hebt gezien. De eerste heet collide_detection en detecteert of twee deeltjes onderling overlappen en dus botsen. De tweede heet particle_collision en berekent de nieuwe snelheden van de twee deeltjes ten gevolge van hun botsing.
# Maken van de class
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 """
""" Geeft TRUE als de deeltjes overlappen """
dx = p1.r[0]-p2.r[0]
dy = p1.r[1]-p2.r[1]
rr = p1.R + p2.R
return dx**2+dy**2 < rr**2
#Note: np.linalg.norm(p1.r - p2.r) < (p1.R + p2.R) had ook gekund, maar is veel langzamer. Controleer zelf!
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_rVolume en randvoorwaarden¶
Voor het volume maken we gebruik van het oppervlak van een pyplot, net zoals we dat in Q1 gedaan hebben. Daar is op dit ogenblik niets anders voor nodig dan een functie die de deeltjes laat botsen met de wanden.
Je hebt hiervoor al een voorbeeld van gezien, maar hieronder gebruiken we een andere vorm. Als je code ontwikkelt is het zaak om regelmatig en structureel te controleren of je code klopt. Het maakt daarbij niet uit of de code door jezelf is gemaakt, door een ander (zoals de docent) is aangeboden of van een generatieve AI komt. Hoe langer je code toevoegt zonder een controle uit te voeren, hoe groter de kans is dat je fout op fout stapelt en niet meer kan traceren wat er precies mis is.
We maken gebruik van een truc die in veel simulaties wordt gebruikt: we maken gebruik van periodieke randvoorwaarden. In dit geval wil dit zeggen dat een deeltje dat botst met de wand aan de ene kant van het volume, aan de andere kant verschijnt met exact dezelfde snelheid. Op deze manier verandert de snelheid van een deeltje niet bij het botsen met de wand en blijven alle behoudswetten gelden. Die behoudswetten kunnen we straks goed controleren.
def box_collision(particle: ParticleClass):
''' botsing met wanden geeft periodieke randvoorwaarden '''
for i in range(2):
if particle.r[i] > BOX_SIZE_0 / 2:
particle.r[i] -= BOX_SIZE_0
elif particle.r[i] < -BOX_SIZE_0 / 2:
particle.r[i] += BOX_SIZE_0Maken van de simulatie¶
Nu dat we alle benodigde functies hebben, kunnen we ze in de juiste structuur zetten om de simulatie uit te voeren.
Bepaling positie deeltjes¶
Eerst moeten we alle deeltjes aanmaken en de startwaardes van hun variabelen invoeren:
# De code in de functie create_particles zorgt ervoor dat:
# 1. Alle deeltjes dezelfde snelheid hebben:
# - vx wordt willekeurig gekozen tussen -V_0 en V_0
# - vy wordt berekend uit vy = ±√(V_0^2 - vx^2) zodat de totale snelheid
# v = sqrt(vx^2 + vy^2) = V_0 constant blijft voor elk deeltje
# - Hierdoor hebben alle deeltjes exact dezelfde snelheid maar in een
# willekeurige richting
#
# 2. De deeltjes uniform verdeeld zijn over de box:
# - pos = np.random.uniform(-BOX_SIZE_0/2 + RADIUS, BOX_SIZE_0/2 - RADIUS, 2)
# kiest een willekeurige x- en y-positie voor elk deeltje binnen de box
# - De straal van het deeltje (RADIUS) wordt in rekening gebracht zodat
# de deeltjes niet op elkaar liggendef create_particles(particles):
"""Leegmaken en opnieuw aanmaken van deeltjes"""
particles.clear()
for i in range(N):
# genereer vx willekeurig in [-V_0, V_0] en bereken vy zo dat |v| ~= V_0
vx = np.random.uniform(-V_0, V_0)
# bescherm tegen kleine afrondfouten met max(..., 0)
vy = np.random.choice([-1, 1]) * np.sqrt(max(V_0**2 - vx**2, 0.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))
Controleren op onderlinge botsingen¶
Om de botsingen van alle deeltjes goed te verwerken moeten we controleren welke paren deeltjes met elkaar moeten botsen en deze botsingen één keer uitvoeren.
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])Controleren op botsingen met wanden¶
Voor het bepalen of deeltjes botsen met de wanden maken we gebruik van de functie die we net hebben aangepast.
def handle_walls(particles):
''' botsing met wanden controleren voor alle deeltjes in volume '''
for p in particles:
box_collision(p)Integreren in tijdstap¶
Nu moeten we alle functies op de juiste manier combineren om de simulatie een tijdstap te laten maken:
def take_time_step(particles):
for p in particles:
p.update_position()
handle_collisions(particles)
handle_walls(particles)
Uitvoeren simulatie en tonen resultaat¶
Zoals aangegeven, kunnen we een animatie maken van de positie en snelheid als functie van de tijd, maar we kunnen ook het eindresultaat tonen en interpreteren:
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_aspect('equal')
dot, = ax.plot([], [], 'ro', ms=13);
def init():
dot.set_data([], [])
return dot,
def update(frame):
take_time_step(particles)
dot.set_data([p.r[0] for p in particles], [p.r[1] for p in particles])
return dot,
ani = FuncAnimation(fig, update, frames=100, init_func=init, blit=True, interval=50)
HTML(ani.to_jshtml())
Toevoegen animatie¶
Zoals gezegd: het kost veel meer rekenkracht om een animatie te maken over de tijd. Soms kan het wel helpen om te zien of de simulatie aan alle verwachtingen voldoet.
from matplotlib.animation import FuncAnimation
from IPython.display import HTML
particles = []
create_particles(particles)
fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_ylim(-BOX_SIZE_0/2, BOX_SIZE_0/2)
ax.set_aspect('equal')
scat = ax.scatter([], [], c='k', s=30)
def init():
scat.set_offsets(np.empty((0,2)))
return scat,
def update(frame):
take_time_step(particles)
coords = np.array([p.r for p in particles])
scat.set_offsets(coords)
return scat,
ani = FuncAnimation(fig, update, frames=300, init_func=init, blit=True, interval=50)
plt.close(fig)
HTML(ani.to_jshtml())
Controle code¶
Om te controleren of onze code werkt passen we nu de code aan.
# Lijst om totale kinetische energie per tijdstap op te slaan
total_energy = []
# Simulatie met bestaande particles van de animatie
for i in range(100):
take_time_step(particles) # update posities en snelheden
# Tel totale kinetische energie
total_energy.append(sum(p.kin_energy for p in particles))
# Plotten
plt.figure(figsize=(8,5))
plt.plot(total_energy, marker='o')
plt.xlabel("Tijdstap")
plt.ylabel("Totale kinetische energie")
plt.title("Totale kinetische energie vs Tijdstap")
plt.grid(True)
plt.show()
#10
# Een andere grootheid die behouden moet blijven is het totale momentum van het systeem.
# Dit geldt omdat we alleen elastische botsingen uitvoeren en er geen externe krachten zijn.
# Kleine variaties kunnen optreden door:
# - Numerieke afrondingsfouten bij berekening van nieuwe snelheden
# - Discrete tijdstap (DT), waardoor botsingen niet exact continu worden berekend
# Code om totale momentum per tijdstap te berekenen en te plotten
total_momentum = []
for i in range(100):
take_time_step(particles)
# Totale momentum vector optellen (2D)
P_tot = np.sum([p.momentum for p in particles], axis=0)
total_momentum.append(np.linalg.norm(P_tot)) # grootte van het momentum vector
# Plot van totale momentum
plt.figure(figsize=(8,5))
plt.plot(total_momentum, marker='o')
plt.xlabel("Tijdstap")
plt.ylabel("Totale momentum")
plt.title("Totale momentum vs Tijdstap")
plt.grid(True)
plt.show()
#11
# De functie handle_collisions doet het volgende:
# 1. Bepaalt het aantal deeltjes in de lijst 'particles'.
# 2. Loopt over alle paren van deeltjes (i,j) waarbij j>i.
# Hierdoor wordt elk paar precies één keer gecontroleerd (geen dubbele botsing).
# 3. Controleert met collide_detection of de twee deeltjes elkaar overlappen.
# 4. Als dat zo is, wordt particle_collision aangeroepen om de snelheden aan te passen volgens een elastische botsing.
#
# Belangrijk:
# - Als drie deeltjes tegelijk zouden botsen, behandelt de code ze paargewijs per tijdstap.
# - Dit betekent dat in een tijdstap nooit exact drie-deeltjes-botsing wordt berekend, alleen paren.
# - Dit is acceptabel bij veel deeltjes omdat triple-botsingen zeldzaam zijn en de meeste interacties pairwise plaatsvinden.
# Lijst van deeltjes aanmaken voor test
particles = []
# Drie deeltjes, twee bewegen naar het derde
particles.append(ParticleClass(m=1.0, v=[1, 0], r=[-1, 0], R=0.1))
particles.append(ParticleClass(m=1.0, v=[-1, 0], r=[1, 0], R=0.1))
particles.append(ParticleClass(m=1.0, v=[0, 0], r=[0, 0], R=0.1))
# Lijst om het aantal botsingen per tijdstap bij te houden
collision_count = []
# Simulatie over een aantal tijdstappen
for t in range(10):
# Controleer en verwerk botsingen
handle_collisions(particles)
# Tel het aantal overlappende paren in deze tijdstap
count = 0
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]):
count += 1
collision_count.append(count)
# Update posities
for p in particles:
p.update_position()
# Plot van aantal botsingen per tijdstap
plt.figure(figsize=(7,5))
plt.plot(collision_count, marker='o')
plt.xlabel("Tijdstap")
plt.ylabel("Aantal botsingen")
plt.title("Aantal botsingen per tijdstap")
plt.grid(True)
plt.show()
import numpy as np
import matplotlib.pyplot as plt
# --- Instellingen ---
dt = 0.01 # tijdstap
T = 5 # totale simulatie tijd
radius = 0.1 # straal van deeltjes voor botsing
# --- Initialisatie van deeltjes ---
# Posities: x, y
particles = np.array([
[0.0, 0.0], # Deeltje A
[1.0, 0.0], # Deeltje B
[0.5, 0.5] # Deeltje C (doel)
])
# Snelheden: vx, vy
v = 0.2
velocities = np.array([
[v, v], # Deeltje A beweegt diagonaal op C af
[-v, v], # Deeltje B beweegt diagonaal op C af
[0.0, 0.0] # Deeltje C staat stil
])
# Opslaan voor plotting
positions = [particles.copy()]
# --- Simulatie ---
collisions = 0
for t in np.arange(0, T, dt):
# Update posities
particles += velocities * dt
# Controleer botsingen
# Alleen paargewijze afstand meten
for i in range(len(particles)):
for j in range(i+1, len(particles)):
dist = np.linalg.norm(particles[i] - particles[j])
if dist < radius*2: # botsing
collisions += 1
# eenvoudige elastische botsing (omkeren van snelheden)
velocities[i], velocities[j] = velocities[j], velocities[i]
positions.append(particles.copy())
positions = np.array(positions)
# --- Plot trajecten ---
plt.figure(figsize=(6,6))
for i in range(len(particles)):
plt.plot(positions[:,i,0], positions[:,i,1], label=f'Deeltje {i+1}')
plt.scatter(particles[:,0], particles[:,1], color='red') # laatste posities
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.title('Trajecten van deeltjes')
plt.axis('equal')
plt.show()
print(f"Totaal aantal botsingen: {collisions}")
Totaal aantal botsingen: 642
# Opslaan botsingen per tijdstap
collision_timeline = np.zeros(len(np.arange(0, T, dt)))
collisions = 0
for k, t in enumerate(np.arange(0, T, dt)):
particles += velocities * dt
step_collisions = 0
for i in range(len(particles)):
for j in range(i+1, len(particles)):
dist = np.linalg.norm(particles[i] - particles[j])
if dist < radius*2:
step_collisions += 1
velocities[i], velocities[j] = velocities[j], velocities[i]
collision_timeline[k] = step_collisions
plt.figure()
plt.plot(np.arange(0, T, dt), collision_timeline)
plt.xlabel('Tijd')
plt.ylabel('Aantal botsingen per stap')
plt.title('Botsingen per tijdseenheid')
plt.show()
