Doelen¶
Nu we de microscopische grootheden van de moleculen hebben verbonden aan de macroscopische grootheden van het gas kunnen we de thermodynamica van het gas echt bestuderen met onze simulatie. In dit werkblad kijken we hoe de temperatuur en de druk veranderen onder invloed van een zuiger die het volume verandert.
Eerst herhalen we de delen van de code die we nodig hebben:
klasse voor het deeltje met bijbehorende functies
variabelen en randcondities van de controle volume
functies voor (een lijst) deeltjes
Daarna voegen we code toe voor de dynamiek van de zuiger:
zuiger implementeren in volume en dynamische formules
En vervolgens:
bestuderen van temperatuur en druk als functie van volume
onderzoeken of we terug kunnen keren naar startcondities
In onderstaande animatie laten we het proces zien dat je gaat programmeren.
Laden van eerdere code¶
We beginnen weer met de noodzakelijke pakketten en de constanten. Daar voegen we nu een constante aan toe: de startsnelheid van de zuiger.
# ruimte voor uitwerking
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.optimize import curve_fit
BOX_SIZE_0 = 10 # Hoogte en lengte startvolume
N = 40 # Aantal deeltjes
V_0 = 1 # Startsnelheid van deeltjes
RADIUS = 0.3 # Straal van moleculen
DT = 0.1 * RADIUS / V_0 # Tijdstap om geen botsing te missen
V_PISTON_0 = -0.1 * V_0 # Startsnelheid van zuiger [nm/ps]
# (negatief betekent zowel links als rechts naar binnen gericht)
#your code/answer
Zoals altijd laden we de klasse voor de gasmoleculen en de functies voor hun onderlinge interactie:
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 """
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
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)
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_rHet volume en de randvoorwaarden zullen we moeten aanpassen aan onze simulatie met bewegende zuiger: Het volume zal nu niet meer altijd een vierkant zijn.
De simulatie bestaat uit een volume met links en rechts een bewegende wand: de zuiger.
Laten we aannemen dat de zuiger altijd in de horizontale richting verplaatst en het volume symmetrisch houdt ten opzichte van de oorsprong, d.w.z. er is een zuiger aan de linker wand die een tegengestelde verplaatsing heeft aan die in de rechter wand.
We maken eerst een aantal variabelen aan die bij het volume horen:
box_height = BOX_SIZE_0
box_length = BOX_SIZE_0
impulse_outward = 0.0
pressure = 0.0
v_piston = V_PISTON_0
work = 0.0De functies die bij het volume en de randvoorwaarden horen moeten we een klein beetje aanpassen, zodat we niet langer uitgaan van de constante waarde van de lengte en hoogte. Om de variabelen zoals box_height en box_length die we hierboven gedefinieerd hebben, later in functies te gebruiken, moeten we ze telkens oproepen met het keyword global. Dit is hieronder uitgewerkt.
def top_down_collision(particle: ParticleClass):
""" botsingen met wanden onder en boven controleren en totale stoot bepalen """
global impulse_outward, box_height
if abs(particle.r[1]) + particle.R > box_height / 2:
particle.r[1] = np.sign(particle.r[1]) * (box_height/2 - particle.R)
impulse_outward += abs(particle.momentum[1]) * 2
particle.v[1] *= -1
def left_right_collision(particle: ParticleClass):
""" verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
global box_length, v_piston, impulse_outward, work
if abs(particle.r[0]) + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = np.sign(particle.r[0]) * v_piston
relative_velocity = particle.v[0] - piston_velocity # stelsel zuiger
particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
impulse_outward += 2 * particle.m * abs(relative_velocity)
work += 2 * particle.m * relative_velocity * piston_velocityEn dan laden we ook alle functies die over de gehele lijst met deeltjes werken, waarbij een paar kleine aanpassingen nodig zijn vanwege de splitsing in het botsen met de wanden. Ook hier roepen we een aantal variabelen met global aan.
def create_particles(particles):
""" Leegmaken en opnieuw aanmaken van deeltjes in lijst """
global box_length, box_height
particles.clear()
for _ in range(N):
vx = np.random.uniform(-V_0, V_0)
vy = np.random.choice([-1, 1]) * np.sqrt(V_0**2 - vx**2)
x = np.random.uniform(-box_length/2 + RADIUS, box_length/2 - RADIUS)
y = np.random.uniform(-box_height/2 + RADIUS, box_height/2 - RADIUS)
particles.append(ParticleClass(m=1.0, v=[vx, vy], r=[x, y], R=RADIUS))
def temperature(particles) -> float:
""" Berekent temperatuur uit gemiddelde kinetische energie """
total_kin_energy = sum(p.kin_energy for p in particles)
avg_kin_energy = total_kin_energy / len(particles)
# Voor 2D: <Ekin> = kB * T, dus T = <Ekin> / kB
# We nemen kB = 1 in onze eenheden
return avg_kin_energy
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 en gemiddeld bepaling druk """
global pressure, impulse_outward, box_height, box_length
impulse_outward = 0.0
for p in particles:
left_right_collision(p)
top_down_collision(p)
pressure = 0.95 * pressure + 0.05 * impulse_outward / ((2 * box_length + 2 * box_height) * DT)
def take_time_step(particles):
""" zet tijdstap voor een lijst deeltjes en verwerk alle botsingen onderling en met wanden """
global box_length, v_piston
box_length += 2 * v_piston * DT # zowel links als rechts zuiger
for p in particles:
p.update_position()
handle_walls(particles)
handle_collisions(particles)
Implementeren (symmetrische) zuiger¶
Voordat we nog meer veranderingen aan de code doorvoeren, moeten we eerst controleren of alles nog werkt. Onderstaande functie is een beetje aangepast ten opzichte van vorige werkbladen, omdat we merkten dat de vorm van de pijlen wel eens de mist in ging bij een heel andere keuze voor eenheden in de constanten.
print("Test simulatie met bewegende zuiger...")
particles = []
create_particles(particles)
for i in range(100):
take_time_step(particles)
plt.figure(figsize=(8, 8))
plt.xlabel('x [nm]')
plt.ylabel('y [nm]')
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=25)
plt.arrow(p.r[0], p.r[1], p.v[0]*DT*30, p.v[1]*DT*30, width=.2*RADIUS,
head_width=RADIUS, head_length=RADIUS, color='red')
plt.title('Deeltjes na 100 tijdstappen')
plt.show()Test simulatie met bewegende zuiger...
Nu implementeren we de zuiger door het toegestane gebied voor de gasdeeltjes bij elke tijdstap te verkleinen met een stap . De factor 2 is opgenomen omdat zowel de linker en de rechter wand een zuigerwand zijn.
def take_time_step(particles):
""" zet tijdstap voor een lijst deeltjes en verwerk alle botsingen onderling en met wanden """
global box_length, v_piston
box_length += 2 * v_piston * DT # zowel links als rechts zuiger
for p in particles:
p.update_position()
handle_walls(particles)
handle_collisions(particles)Hieronder draaien we een kleine simulatie om te kijken of we de box kleiner zien worden en vervolgens te bestuderen hoe de temperatuur zich gedraagt als functie van het oppervlak/volume.
# Deel van de animated simulatie
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=14)
def init():
dot.set_data([], [])
return dot,
# we kiezen het aantal datapunten zodat het volume tot 1/3 van het begin volume reduceert
num_steps = round(4/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
box_length = BOX_SIZE_0 # zetten zuiger terug
v_piston = V_PISTON_0
create_particles(particles) # resetten deeltjes
def update(frame):
take_time_step(particles)
dot.set_data([p.r[0] for p in particles], [p.r[1] for p in particles])
volumes[frame] = box_length * box_height # gebruik 'frame' in plaats van 'i'
temperatures[frame] = temperature(particles) # gebruik 'frame' in plaats van 'i'
return dot,
ani = FuncAnimation(fig, update, frames=int(num_steps/2), init_func=init, blit=True, interval=50)
HTML(ani.to_jshtml())Animation size has reached 20974053 bytes, exceeding the limit of 20971520.0. If you're sure you want a larger animation embedded, set the animation.embed_limit rc parameter to a larger value (in MB). This and further frames will be dropped.
Loading...
# we kiezen het aantal datapunten zodat het volume tot 1/3 van het begin volume reduceert
num_steps = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
box_length = BOX_SIZE_0 # zetten zuiger terug
v_piston = V_PISTON_0
create_particles(particles) # resetten deeltjes
for i in range(num_steps):
take_time_step(particles)
volumes[i] = box_length * box_height
temperatures[i] = temperature(particles)
temperatures = np.asarray(temperatures)
plt.figure()
plt.xlabel('Volume')
plt.ylabel('Temperature')
plt.plot(volumes, temperatures, '-r')
plt.show()
Dit kan niet kloppen. Hier zien we dat de temperatuur vrijwel constant is (let op de vermenigvuldigingsfactor vermeld aan de bovenkant van de verticale as). Maar de zuiger voert arbeid uit, op baiss van de wet van behoud van energie betekent zou de temperatuur moet veranderen!
Om het model kloppend te maken moeten we kijken naar de botsing van de deeltjes met de wand. In de vorige werkbladen stonden de wanden stil en veranderde de snelheid van de deeltjes alleen van teken in de component loodrecht op de wand. Nu dat de wanden een zuiger zijn en snelheid hebben, moeten we daarvoor corrigeren.
De snelheid van de deeltjes klapt nog steeds om van teken in het referentiestelsel van de wand, maar omdat de wand beweegt ten opzichte van het volume met snelheid , wordt de juiste functie:
def left_right_collision(particle: ParticleClass):
""" verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
global box_length, v_piston, impulse_outward
if abs(particle.r[0]) + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = np.sign(particle.r[0]) * v_piston
relative_velocity = particle.v[0] - piston_velocity # stelsel zuiger
particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
impulse_outward += 2 * particle.m * abs(relative_velocity) # stoot gevoeld door zuigerNu kunnen we de simulatie opnieuw uitvoeren:
num_steps = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
box_length = BOX_SIZE_0
v_piston = V_PISTON_0
pressure = 0.0
work = 0.0
create_particles(particles)
for i in range(num_steps):
take_time_step(particles)
volumes[i] = box_length * box_height
temperatures[i] = temperature(particles)
plt.figure(figsize=(10, 6))
plt.xlabel('Volume [nm²]')
plt.ylabel('Temperatuur [u·nm²/ps²]')
plt.title('Temperatuur als functie van Volume')
plt.plot(volumes, temperatures, '-r', linewidth=2)
plt.grid(True, alpha=0.3)
plt.show()
We zien nu een heel duidelijke afhankelijkheid van de temperatuur op het volume, zoals we ook verwachten vanwege de wet van behoud van energie. Een volgende logische stap is of deze grafiek ook daadwerkelijk overeenkomt met onze verwachting, maar daarvoor is het waardevol om ook informatie te halen uit de andere grootheden.
num_steps = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))
particles = []
volumes = np.zeros(num_steps, dtype=float)
pressures = np.zeros(num_steps, dtype=float)
pressure = 0.0
box_length = BOX_SIZE_0
v_piston = V_PISTON_0
work = 0.0
create_particles(particles)
for i in range(num_steps):
take_time_step(particles)
volumes[i] = box_length * box_height
pressures[i] = pressure
plt.figure(figsize=(10, 6))
plt.xlabel('Volume [nm²]')
plt.ylabel('Druk [u/(nm·ps²)]')
plt.title('Druk als functie van Volume')
plt.plot(volumes, pressures, '-b', linewidth=2)
plt.grid(True, alpha=0.3)
plt.show()
Als je simulatie klopt heeft deze grafiek de vorm van een machtsfunctie met een negatieve exponent.
def power_law(vol, a, n):
""" de fitfunctie voor het P,V-diagram """
return a * (vol)**n
# Fit uitvoeren met goede beginwaardes
p0 = [pressures[0] * volumes[0]**1.4, -3] # Verwachting: n ≈ -γ ≈ -1.4 voor 2D
bounds = ((0, -5), (np.inf, -1))
popt, pcov = curve_fit(power_law, volumes, pressures, p0=p0, bounds=bounds)
a_fit, n_fit = popt
print(f"\nFit resultaten:")
print(f"a = {a_fit:.3f}")
print(f"n = {n_fit:.3f}")
print(f"Verwachte waarde n voor adiabatisch proces in 2D: ≈ -1.4 tot -1.5")
plt.figure(figsize=(10, 6))
plt.xlabel('Volume [nm²]')
plt.ylabel('Druk [u/(nm·ps²)]')
plt.title(f'Druk vs Volume met fit: P = {a_fit:.2f} × V^{n_fit:.3f}')
plt.plot(volumes, pressures, 'b.', label='Simulatie data', alpha=0.6)
plt.plot(volumes, power_law(volumes, a_fit, n_fit), 'r-',
label=f'Fit: n = {n_fit:.3f}', linewidth=2)
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
Fit resultaten:
a = 2033022.323
n = -3.534
Verwachte waarde n voor adiabatisch proces in 2D: ≈ -1.4 tot -1.5
Als je de simulatie een aantal keer uitvoert, zie je dat er een structureel verschil zit tussen de waarde die je verwacht en de waarde die je uit de fit krijgt.
Eerste hoofdwet¶
Een goede controle voor ons simulatie is te onderzoeken of deze voldoet aan de eerste hoofdwet van de thermodynamica.
Omdat er (nog) geen warmte wordt uitgewisseld, moeten we alleen de hoeveelheid arbeid bepalen. De arbeid wordt geleverd door de zuiger(s) op de deeltjes. Dus in de functie waar we de botsingen van de deeltjes met de wand behandelen, kunnen we ook de verrichte arbeid bepalen. De arbeid wordt gegeven\ door:
Hier staan de 1 voor de begintoestand en de 2 voor de eindtoestand van het proces. In de differentiaalvorm wordt dit geschreven als:
Zoals dat in het boek gebeurt, kiezen we voor de notatie met in plaats van , om aan te geven dat deze integraal afhankelijk is van het gekozen proces en niet alleen van het begin- en eindpunt. Voor ons experiment, waar het volume niet geleidelijk maar in stapjes verandert, kan je de hoeveelheid arbeid per tijdstip dus herschrijven tot:
waarin de hoogte is van de zuiger.
Dit betekent dat we de arbeid per tijdstap in onze code kunnen bepalen op hetzelfde moment dat we de druk bepalen met behulp van de gezamenlijke stoot van de moleculen. Daarvoor moeten we de code voor de functie left_right_collision nogmaals aanpassen:
def left_right_collision(particle: ParticleClass):
""" verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
global box_length, v_piston, impulse_outward, work
if abs(particle.r[0]) + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = np.sign(particle.r[0]) * v_piston
relative_velocity = particle.v[0] - piston_velocity # stelsel zuiger
particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
impulse_outward += 2 * particle.m * abs(relative_velocity)
work += 2 * particle.m * relative_velocity * piston_velocityVoordat we een simulatie draaien is het goed even stil te staan en na te denken. Net als bij een experiment in het lab doorlopen we bewust de onderzoekscyclus. We moeten daarvoor eerst een hypothese opstellen en daarna pas de simulatie uitvoeren.
num_steps = 1000
particles = []
pressure = 0.0
work_array = np.zeros(num_steps, dtype=float)
energy_array = np.zeros(num_steps, dtype=float)
work = 0.0
box_length = BOX_SIZE_0
v_piston = V_PISTON_0
create_particles(particles)
initial_energy = sum(p.kin_energy for p in particles)
for i in range(num_steps):
take_time_step(particles)
work_array[i] = work
current_energy = sum(p.kin_energy for p in particles)
energy_array[i] = current_energy - initial_energy
plt.figure(figsize=(10, 6))
plt.xlabel('Stap')
plt.ylabel('Energie [u·nm²/ps²]')
plt.title('Eerste Hoofdwet: ΔE_kin vs Arbeid (W)')
plt.plot(work_array, label='Arbeid door zuiger', linewidth=2)
plt.plot(energy_array, label='Verandering kinetische energie',
linewidth=2, linestyle='--')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
print(f"\nControle Eerste Hoofdwet:")
print(f"Totale arbeid: {work_array[-1]:.3f}")
print(f"Verandering kinetische energie: {energy_array[-1]:.3f}")
print(f"Verschil: {abs(work_array[-1] - energy_array[-1]):.3f}")
Controle Eerste Hoofdwet:
Totale arbeid: -71.525
Verandering kinetische energie: 71.525
Verschil: 143.050
# Validatie van de Eerste Hoofdwet: ΔU = W (voor adiabatisch proces)
import numpy as np
import matplotlib.pyplot as plt
# Constante voor Boltzmann (in SI eenheden)
k_B = 1.380649e-23 # J/K
num_steps = 1000
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
works = np.zeros(num_steps, dtype=float)
internal_energies = np.zeros(num_steps, dtype=float)
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0
v_piston = V_PISTON_0
create_particles(particles)
# Initiële interne energie
T_0 = temperature(particles)
U_0 = N * k_B * T_0 # Voor 2D ideaal gas: U = N*k_B*T
for i in range(num_steps):
take_time_step(particles)
volumes[i] = box_length * box_height
temperatures[i] = temperature(particles)
works[i] = work
# Interne energie relatief t.o.v. begin
internal_energies[i] = N * k_B * temperatures[i] - U_0
# Plot resultaten
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Plot 1: ΔU vs W (hoofdvalidatie)
ax1 = axes[0, 0]
ax1.plot(works, internal_energies, 'b-', linewidth=2, label='Simulatie data')
ax1.plot([works[0], works[-1]], [works[0], works[-1]], 'r--',
linewidth=2, label='ΔU = W (ideaal)')
ax1.set_xlabel('Arbeid W [J]', fontsize=11)
ax1.set_ylabel('Verandering interne energie ΔU [J]', fontsize=11)
ax1.set_title('Validatie Eerste Hoofdwet (adiabatisch proces)', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Bereken en toon afwijking
final_deviation = abs(internal_energies[-1] - works[-1])
relative_error = (final_deviation / abs(works[-1])) * 100 if works[-1] != 0 else 0
ax1.text(0.05, 0.95, f'Relatieve fout: {relative_error:.2f}%',
transform=ax1.transAxes, fontsize=10, verticalalignment='top',
bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
# Plot 2: ΔU en W als functie van tijd
ax2 = axes[0, 1]
time_steps = np.arange(num_steps)
ax2.plot(time_steps, internal_energies, 'b-', linewidth=2, label='ΔU')
ax2.plot(time_steps, works, 'r--', linewidth=2, label='W')
ax2.set_xlabel('Tijdstap', fontsize=11)
ax2.set_ylabel('Energie [J]', fontsize=11)
ax2.set_title('ΔU en W in de tijd', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
# Plot 3: Temperatuur vs Volume
ax3 = axes[1, 0]
ax3.plot(volumes, temperatures, 'g-', linewidth=2)
ax3.set_xlabel('Volume [m²]', fontsize=11)
ax3.set_ylabel('Temperatuur [K]', fontsize=11)
ax3.set_title('T-V diagram', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
# Plot 4: Percentage afwijking in de tijd
ax4 = axes[1, 1]
# Vermijd deling door nul
valid_indices = works != 0
percentage_diff = np.zeros_like(works)
percentage_diff[valid_indices] = ((internal_energies[valid_indices] - works[valid_indices]) /
np.abs(works[valid_indices])) * 100
ax4.plot(time_steps, percentage_diff, 'purple', linewidth=2)
ax4.axhline(y=0, color='k', linestyle='--', alpha=0.5)
ax4.set_xlabel('Tijdstap', fontsize=11)
ax4.set_ylabel('Relatieve afwijking (%)', fontsize=11)
ax4.set_title('(ΔU - W) / |W| × 100%', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Numerieke samenvatting
print("VALIDATIE EERSTE HOOFDWET")
print(f"Begintoestand:")
print(f" Volume: {BOX_SIZE_0**2:.6e} m²")
print(f" Temperatuur: {T_0:.2f} K")
print(f" Interne energie: {U_0:.6e} J")
print(f"\nEindtoestand:")
print(f" Volume: {volumes[-1]:.6e} m²")
print(f" Temperatuur: {temperatures[-1]:.2f} K")
print(f" ΔU: {internal_energies[-1]:.6e} J")
print(f" W: {works[-1]:.6e} J")
print(f"\nValidatie:")
print(f" |ΔU - W|: {final_deviation:.6e} J")
print(f" Relatieve fout: {relative_error:.3f}%")
if relative_error < 5:
print(f" ✓ Simulatie voldoet goed aan de eerste hoofdwet!")
elif relative_error < 10:
print(f" ~ Simulatie voldoet redelijk aan de eerste hoofdwet")
else:
print(f" ✗ Grote afwijking - controleer simulatieparameters")
plt.figure()
plt.plot(works[0:999], temperatures[0:999], "r")
plt.plot(works[1000:1999], temperatures[1000:1999], "g")
plt.show()
VALIDATIE EERSTE HOOFDWET
Begintoestand:
Volume: 1.000000e+02 m²
Temperatuur: 0.50 K
Interne energie: 2.761298e-22 J
Eindtoestand:
Volume: 4.000000e+01 m²
Temperatuur: 2.27 K
ΔU: 9.769010e-22 J
W: -7.075665e+01 J
Validatie:
|ΔU - W|: 7.075665e+01 J
Relatieve fout: 100.000%
✗ Grote afwijking - controleer simulatieparameters
Reversibiliteit¶
Als laatste onderdeel van dit werkblad maken we een eenvoudige cyclus. We comprimeren het volume gedurende 1000 tijdstappen en keren daarna in 1000 gelijke tijdstappen terug naar het startvolume.
# Conversie constanten
U_TO_KG = 1.66054e-27 # atomaire massa naar kg
NM_TO_M = 1e-9 # nanometer naar meter
PS_TO_S = 1e-12 # picoseconde naar seconde
K_B = 1.380649e-23 # Boltzmann constante [J/K]
def convert_to_SI(work_sim, temp_sim, N):
"""
Converteer simulatie-eenheden naar SI-eenheden
Parameters:
- work_sim: arbeid in u·nm²/ps²
- temp_sim: temperatuur in u·nm²/ps²
- N: aantal deeltjes
Returns:
- work_J: arbeid in Joule
- temp_K: temperatuur in Kelvin
"""
# Arbeid: u·nm²/ps² → J
# 1 u·nm²/ps² = (1.66054e-27 kg)·(1e-9 m)²/(1e-12 s)² = 1.66054e-3 J
work_J = work_sim * U_TO_KG * (NM_TO_M**2) / (PS_TO_S**2)
# Temperatuur: u·nm²/ps² → K
# T_sim = (2/3) * E_kin_gemiddeld = (2/3) * (1/2) * m * v²
# In SI: E_kin [J] = T_sim [u·nm²/ps²] * conversie
# T [K] = E_kin / k_B
temp_J = temp_sim * U_TO_KG * (NM_TO_M**2) / (PS_TO_S**2)
temp_K = temp_J / K_B
return work_J, temp_Knum_steps = 1000
particles = []
work_cycle = np.zeros(2*num_steps) # Arbeid opslaan
temps_cycle = np.zeros(2*num_steps) # Temperatuur opslaan
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0
v_piston = V_PISTON_0
create_particles(particles)
# Compressie
for i in range(num_steps):
take_time_step(particles)
work_cycle[i] = work # ← Arbeid opslaan (niet volume!)
temps_cycle[i] = temperature(particles) # ← Temperatuur opslaan (niet druk!)
# Omkeren zuiger
v_piston = -V_PISTON_0
# Decompressie
for i in range(num_steps, 2*num_steps):
take_time_step(particles)
work_cycle[i] = work # ← Arbeid opslaan
temps_cycle[i] = temperature(particles) # ← Temperatuur opslaan
# Converteer naar SI
work_cycle_J, temps_cycle_K = convert_to_SI(work_cycle, temps_cycle, N=40)
# Plot cyclus 1
plt.figure(figsize=(10, 8))
plt.plot(temps_cycle_K[:num_steps], work_cycle_J[:num_steps] * 1e3,
'b-', label='Compressie', linewidth=2)
plt.plot(temps_cycle_K[num_steps:], work_cycle_J[num_steps:] * 1e3,
'r-', label='Decompressie', linewidth=2)
plt.xlabel('Temperatuur [K]', fontsize=12)
plt.ylabel('Arbeid [mJ]', fontsize=12)
plt.title('W-T diagram: Cyclus (langzame zuiger)', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
num_steps_fast = 200
particles = []
work_fast = np.zeros(2*num_steps_fast) # Arbeid opslaan
temps_fast = np.zeros(2*num_steps_fast) # Temperatuur opslaan
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0
v_piston = 5 * V_PISTON_0
create_particles(particles)
# Compressie
for i in range(num_steps_fast):
take_time_step(particles)
work_fast[i] = work # ← Arbeid opslaan
temps_fast[i] = temperature(particles) # ← Temperatuur opslaan
# Omkeren zuiger
v_piston = -5 * V_PISTON_0
# Decompressie
for i in range(num_steps_fast, 2*num_steps_fast):
take_time_step(particles)
work_fast[i] = work # ← Arbeid opslaan
temps_fast[i] = temperature(particles) # ← Temperatuur opslaan
# Converteer naar SI
work_fast_J, temps_fast_K = convert_to_SI(work_fast, temps_fast, N=40)
# Plot cyclus 2
plt.figure(figsize=(10, 8))
plt.plot(temps_fast_K[:num_steps_fast], work_fast_J[:num_steps_fast] * 1e3,
'b-', label='Compressie', linewidth=2)
plt.plot(temps_fast_K[num_steps_fast:], work_fast_J[num_steps_fast:] * 1e3,
'r-', label='Decompressie', linewidth=2)
plt.xlabel('Temperatuur [K]', fontsize=12)
plt.ylabel('Arbeid [mJ]', fontsize=12)
plt.title('W-T diagram: Cyclus (snelle zuiger, 5x)', fontsize=14)
plt.legend(fontsize=11)
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\nZuigersnelheid in simulatie: {abs(5 * V_PISTON_0):.1f} nm/ps")
print(f"Zuigersnelheid in SI: {abs(5 * V_PISTON_0) * 1000:.1f} m/s")
print("\nBij snelle zuiger (5x) zie je dat het systeem niet terugkeert")
print("naar de begintoestand - dit toont non-equilibrium gedrag.")
Zuigersnelheid in simulatie: 0.5 nm/ps
Zuigersnelheid in SI: 500.0 m/s
Bij snelle zuiger (5x) zie je dat het systeem niet terugkeert
naar de begintoestand - dit toont non-equilibrium gedrag.
In deze laatste simulatie zie je (bij een correcte code) dat nog altijd netjes wordt voldaan aan de eerste hoofdwet. Ondanks dat, zie je dat het systeem niet terugkeert naar de begintoestand. In het boek wordt dit omschreven in het deel over ‘quasi-equilibrium’: Processen moeten voldoende traag plaatsvinden zodat er zich een evenwicht kan vormen binnen het gas. Als processen te snel plaatsvinden is er geen sprake van equilibrium en kun je de macroscopische thermodynamische formules niet langer gebruiken.
"""
Uitbreiding: Parallel bewegende zuigers
Beide zuigers bewegen in dezelfde richting → volume blijft constant
Onderzoek: effect op temperatuur
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
# Nieuwe functie: parallelle zuigercollisie
def left_right_collision_parallel(particle: ParticleClass):
"""Botsingen met wanden die parallel bewegen (zelfde snelheid)"""
global box_length, v_piston_parallel, impulse_outward, work, box_center_x
# Beide wanden bewegen met dezelfde snelheid, dus centrum verschuift
if abs(particle.r[0] - box_center_x) + particle.R > box_length / 2:
# Corrigeer positie relatief t.o.v. bewegend centrum
particle.r[0] = box_center_x + np.sign(particle.r[0] - box_center_x) * (box_length/2 - particle.R)
# Botsing berekenen
piston_velocity = v_piston_parallel # Beide zuigers zelfde richting
relative_velocity = particle.v[0] - piston_velocity
particle.v[0] = -relative_velocity + piston_velocity
impulse_outward += 2 * particle.m * abs(relative_velocity)
# Netto arbeid is nul (beide zuigers doen werk, maar tegengesteld effect)
# Wel is er energieoverdracht via snelheidsverandering moleculen
def take_time_step_parallel(particles):
"""Tijdstap met parallelle zuigers"""
global box_center_x, v_piston_parallel
# Centrum verschuift, volume blijft constant
box_center_x += v_piston_parallel * DT
for p in particles:
p.update_position()
handle_walls_parallel(particles)
handle_collisions(particles)
def handle_walls_parallel(particles):
"""Wandbotsingen met parallelle zuigers"""
global pressure, impulse_outward, box_height, box_length
impulse_outward = 0.0
for p in particles:
left_right_collision_parallel(p)
top_down_collision(p)
pressure = 0.95 * pressure + 0.05 * impulse_outward / ((2 * box_length + 2 * box_height) * DT)
# Simulatie setup
def run_parallel_piston_simulation(piston_speed, initial_temp_factor=1.0, num_steps=2000):
"""
Voer simulatie uit met parallelle zuigers
Parameters:
- piston_speed: snelheid van beide zuigers (zelfde richting)
- initial_temp_factor: factor voor starttemperatuur
- num_steps: aantal tijdstappen
"""
global box_length, box_center_x, v_piston_parallel, pressure, V_0
# Reset variabelen
particles = []
box_length = BOX_SIZE_0
box_center_x = 0.0
v_piston_parallel = piston_speed
pressure = 0.0
V_0_original = V_0
V_0 = V_0_original * np.sqrt(initial_temp_factor) # T ∝ v²
create_particles(particles)
V_0 = V_0_original # Reset voor volgende runs
# Data arrays
temperatures = np.zeros(num_steps)
center_positions = np.zeros(num_steps)
avg_velocities_x = np.zeros(num_steps)
volumes = np.zeros(num_steps)
for i in range(num_steps):
take_time_step_parallel(particles)
temperatures[i] = temperature(particles)
center_positions[i] = box_center_x
avg_velocities_x[i] = np.mean([p.v[0] for p in particles])
volumes[i] = box_length * box_height
return temperatures, center_positions, avg_velocities_x, volumes
# Experiment 1: Effect van zuigersnelheid
print("Experiment 1: Variatie zuigersnelheid")
fig1, axes1 = plt.subplots(2, 2, figsize=(14, 10))
speeds = [0.1 * V_0, 0.5 * V_0, 1.0 * V_0, 2.0 * V_0]
colors = ['blue', 'green', 'orange', 'red']
for speed, color in zip(speeds, colors):
T, center, v_avg, vol = run_parallel_piston_simulation(speed, num_steps=3000)
displacement = center / BOX_SIZE_0 # Genormaliseerd
# Plot 1: Temperatuur vs verplaatsing
axes1[0, 0].plot(displacement, T, color=color, linewidth=2,
label=f'v_piston = {speed/V_0:.1f}×v₀')
# Plot 2: Gemiddelde snelheid vs verplaatsing
axes1[0, 1].plot(displacement, v_avg, color=color, linewidth=2)
print(f"v_piston = {speed/V_0:.1f}×v₀:")
print(f" ΔT = {T[-1]-T[0]:.2f} K ({((T[-1]-T[0])/T[0]*100):.1f}%)")
print(f" Eind <v_x> = {v_avg[-1]:.3f} m/s")
axes1[0, 0].set_xlabel('Verplaatsing centrum (×BOX_SIZE)', fontsize=11)
axes1[0, 0].set_ylabel('Temperatuur [K]', fontsize=11)
axes1[0, 0].set_title('Temperatuur vs Verplaatsing', fontsize=12, fontweight='bold')
axes1[0, 0].legend()
axes1[0, 0].grid(True, alpha=0.3)
axes1[0, 1].set_xlabel('Verplaatsing centrum (×BOX_SIZE)', fontsize=11)
axes1[0, 1].set_ylabel('Gemiddelde v_x [m/s]', fontsize=11)
axes1[0, 1].set_title('Drift snelheid vs Verplaatsing', fontsize=12, fontweight='bold')
axes1[0, 1].grid(True, alpha=0.3)
# Experiment 2: Effect van starttemperatuur
print("\nExperiment 2: Variatie starttemperatuur")
temp_factors = [0.5, 1.0, 2.0, 4.0]
v_piston_test = 0.5 * V_0
for temp_factor, color in zip(temp_factors, colors):
T, center, v_avg, vol = run_parallel_piston_simulation(v_piston_test,
initial_temp_factor=temp_factor,
num_steps=3000)
displacement = center / BOX_SIZE_0
axes1[1, 0].plot(displacement, T, color=color, linewidth=2,
label=f'T₀ = {temp_factor:.1f}×T_ref')
axes1[1, 1].plot(displacement, v_avg, color=color, linewidth=2)
print(f"T₀ factor = {temp_factor:.1f}:")
print(f" T_start = {T[0]:.2f} K")
print(f" ΔT = {T[-1]-T[0]:.2f} K ({((T[-1]-T[0])/T[0]*100):.1f}%)")
axes1[1, 0].set_xlabel('Verplaatsing centrum (×BOX_SIZE)', fontsize=11)
axes1[1, 0].set_ylabel('Temperatuur [K]', fontsize=11)
axes1[1, 0].set_title('Effect starttemperatuur', fontsize=12, fontweight='bold')
axes1[1, 0].legend()
axes1[1, 0].grid(True, alpha=0.3)
axes1[1, 1].set_xlabel('Verplaatsing centrum (×BOX_SIZE)', fontsize=11)
axes1[1, 1].set_ylabel('Gemiddelde v_x [m/s]', fontsize=11)
axes1[1, 1].set_title('Drift snelheid (verschillende T₀)', fontsize=12, fontweight='bold')
axes1[1, 1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# Verklaring
print("VERKLARING VAN HET GEDRAG")
print("""
MICROSCOPISCH NIVEAU (balletjes met impuls en energie):
-------------------------------------------------------
1. Wanneer beide zuigers in dezelfde richting bewegen, verandert het
volume NIET, maar het centrum van het systeem verschuift wel.
2. Moleculen die botsen met de VOORSTE wand (bewegingsrichting zuiger):
- Krijgen EXTRA energie (elastische botsing met naderende wand)
- Hun snelheid in x-richting NEEMT TOE
3. Moleculen die botsen met de ACHTERSTE wand:
- VERLIEZEN energie (botsing met wegbewegende wand)
- Hun snelheid in x-richting NEEMT AF
4. Netto effect:
- Gemiddelde x-snelheid verschuift naar de bewegingsrichting (drift)
- Temperatuur STIJGT omdat de spreiding in snelheden toeneemt
- Dit is een NIET-EVENWICHTSPROCES
MACROSCOPISCH NIVEAU (P, V, T):
--------------------------------
1. Volume V = constant (beide zuigers bewegen gelijk)
2. Druk P: Er is geen netto arbeid op het systeem, maar er is wel
energietoevoer door de asymmetrische botsingen.
3. Temperatuur T: STIJGT omdat:
- Kinetische energie wordt toegevoegd via de voorste wand
- De achterwand kan niet alle energie afvoeren
- Er ontstaat een netto energietoevoer
4. Dit lijkt op "visceuze verwarming": wrijving tussen het gas en de
bewegende wanden zorgt voor temperatuurstijging.
CONCLUSIE:
----------
Bij parallelle zuigerbeweging:
- Snellere zuigers → meer temperatuurstijging
- Hogere T₀ → grotere absolute ΔT
- Het systeem komt NIET in thermisch evenwicht
- Dit is geen reversibel thermodynamisch proces
""")Experiment 1: Variatie zuigersnelheid
v_piston = 0.1×v₀:
ΔT = -0.00 K (-0.2%)
Eind <v_x> = 0.100 m/s
v_piston = 0.5×v₀:
ΔT = 0.21 K (42.2%)
Eind <v_x> = 0.475 m/s
v_piston = 1.0×v₀:
ΔT = 1.18 K (237.0%)
Eind <v_x> = 1.003 m/s
v_piston = 2.0×v₀:
ΔT = 2.75 K (549.7%)
Eind <v_x> = 1.390 m/s
Experiment 2: Variatie starttemperatuur
T₀ factor = 0.5:
T_start = 0.25 K
ΔT = 0.22 K (86.2%)
T₀ factor = 1.0:
T_start = 0.53 K
ΔT = 0.20 K (37.0%)
T₀ factor = 2.0:
T_start = 1.00 K
ΔT = 0.37 K (36.7%)
T₀ factor = 4.0:
T_start = 2.05 K
ΔT = 0.07 K (3.3%)
VERKLARING VAN HET GEDRAG
MICROSCOPISCH NIVEAU (balletjes met impuls en energie):
-------------------------------------------------------
1. Wanneer beide zuigers in dezelfde richting bewegen, verandert het
volume NIET, maar het centrum van het systeem verschuift wel.
2. Moleculen die botsen met de VOORSTE wand (bewegingsrichting zuiger):
- Krijgen EXTRA energie (elastische botsing met naderende wand)
- Hun snelheid in x-richting NEEMT TOE
3. Moleculen die botsen met de ACHTERSTE wand:
- VERLIEZEN energie (botsing met wegbewegende wand)
- Hun snelheid in x-richting NEEMT AF
4. Netto effect:
- Gemiddelde x-snelheid verschuift naar de bewegingsrichting (drift)
- Temperatuur STIJGT omdat de spreiding in snelheden toeneemt
- Dit is een NIET-EVENWICHTSPROCES
MACROSCOPISCH NIVEAU (P, V, T):
--------------------------------
1. Volume V = constant (beide zuigers bewegen gelijk)
2. Druk P: Er is geen netto arbeid op het systeem, maar er is wel
energietoevoer door de asymmetrische botsingen.
3. Temperatuur T: STIJGT omdat:
- Kinetische energie wordt toegevoegd via de voorste wand
- De achterwand kan niet alle energie afvoeren
- Er ontstaat een netto energietoevoer
4. Dit lijkt op "visceuze verwarming": wrijving tussen het gas en de
bewegende wanden zorgt voor temperatuurstijging.
CONCLUSIE:
----------
Bij parallelle zuigerbeweging:
- Snellere zuigers → meer temperatuurstijging
- Hogere T₀ → grotere absolute ΔT
- Het systeem komt NIET in thermisch evenwicht
- Dit is geen reversibel thermodynamisch proces