Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Arbeid

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:

Daarna voegen we code toe voor de dynamiek van de zuiger:

En vervolgens:

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 = 0               # Hoogte en lengte startvolume
N = 40                       # Aantal deeltjes
V_0 = 0                      # Startsnelheid van deeltjes
RADIUS = 0                   # Straal van moleculen
DT = 0                       # Tijdstap om geen botsing te missen
V_PISTON_0 = -0.1 * V_0      # Startsnelheid van zuiger 
# (negatief betekent zowel links als rechts naar binnen gericht)

#your code/answer
BOX_SIZE_0 = 10 #nm              # Hoogte en breedte startvolume
N = 40                            # Aantal deeltjes 
V_0 = 1                           # Startsnelheid van deeltjes (m/s)
RADIUS = 0.3                        # Straal van moleculen (nm)
DT = 0.1 * RADIUS / V_0 #ps                    # Tijdstap om geen botsing te missen
V_PISTON_0 = -0.1 * V_0      # Startsnelheid van zuiger (m/s)
# (negatief betekent zowel links als rechts naar binnen gericht)

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)
    # 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

Het 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.

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     # hoogte van beheersvolume
box_length = BOX_SIZE_0     # breedte van beheersvolume
impulse_outward = 0.0       # totale stoot van deeltjes naar buiten gericht
pressure = 0.0              # druk in beheersvolume
v_piston = V_PISTON_0       # huidige snelheid van zuiger 
work = 0.0                  # arbeid uitgevoerd door gas

De 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):
    """ botsingen met wanden links en rechts controleren en totale stoot bepalen """
    global impulse_outward, box_length
    if abs(particle.r[0]) + particle.R > box_length / 2:
        particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
        impulse_outward += abs(particle.momentum[0]) * 2
        particle.v[0] *= -1

En 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:
    # De oplossing van je vorige werkblad
    """" berekent temperatuur van deeltjes in lijst """
    if len(particles) == 0:
        return 0.0
    if len(particles) >0:
        total_kin_energy = sum(p.kin_energy for p in particles)
        avg_kin_energy = total_kin_energy / len(particles)
        # Gebruik de relatie tussen gemiddelde kinetische energie en temperatuur
        # E_kin = 1/2 m v^2 = 3/2 kB T
        # T = (2/3) * E_kin / kB
        T = (2/3) * avg_kin_energy  # in geschikte eenheden
        return T
         
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   # om pressure buiten de functie te kunnen gebruiken
    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 """
    for p in particles:
        p.update_position()
    handle_collisions(particles)
    handle_walls(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.

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=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.show()
<Figure size 640x480 with 1 Axes>

Nu implementeren we de zuiger door het toegestane gebied voor de gasdeeltjes bij elke tijdstap te verkleinen met een stap 2vpistondt2v_{\text{piston}}dt. 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(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     

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[i] = box_length * box_height
    temperatures[i] = temperature(particles)
    return dot,
    
ani = FuncAnimation(fig, update, frames=int(num_steps/2), init_func=init, blit=True, interval=50)
HTML(ani.to_jshtml())
Loading...
<Figure size 640x480 with 1 Axes>
# 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()
<Figure size 640x480 with 1 Axes>

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 vpistonv_{\text{piston}}, 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 zuiger

Nu 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         # 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)

plt.figure()
plt.xlabel('Volume')
plt.ylabel('Temperature')

temperatures = np.asarray(temperatures)

# Plot de data
plt.plot(volumes, temperatures, '-r', linewidth=2, label='Simulatie data')

# Titels en eenheden 
plt.xlabel('Volume ($nm^2$)')       # want de boxsize is in nm
# T = (2/3) * <E_kin> / kB dus T is in arbitrary units (A.U.) gebaseerd op gemiddelde kinetische energie
# T is dus uitgedrukt in een dimensieloos, relatief getal tov de gemiddelde omgevingstemperatuur
plt.ylabel(r'Temperatuur (A.U. gebaseerd op $\langle E_{kin} \rangle$)')
plt.title('Compressie: T-V diagram')
plt.grid(True, alpha=0.7)
plt.legend()
plt.show()
<Figure size 640x480 with 1 Axes>

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.

# breid deze code uit met jouw antwoord
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)

# reset variables
global pressure
pressure = 0.0
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)
    # RUIMTE VOOR UITWERKING
    # Volume in nm^2 (2D oppervlakte)
    volumes[i] = box_length * box_height
    pressures[i] = pressure
    
#your code/answer
pressures = np.asarray(pressures)

# PLOT GRAFIEK
plt.figure()
plt.xlabel('Volume ($nm^2$)')
plt.ylabel('Druk (A.U.)')
plt.plot(volumes, pressures, '-b')
plt.title('Compressie: P-V diagram')
plt.grid(True, alpha=0.7)
plt.show()
<Figure size 640x480 with 1 Axes>

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   

# RUIMTE VOOR VERDERE UITWERKING
# bounds = ((min_a, min_n), (max_a, max_n))
bounds = ((0, -5), (np.inf, 0))
# Gebruik het eerste datapunt voor de schatting
V_start = volumes[0]
P_start = pressures[0]
n_guess = -2.0

# De schatting voor a is dan P / V^n
a_guess = P_start / (V_start**n_guess)

p0 = [a_guess, n_guess] # beginwaarden voor a en n# initiële schatting voor a en n  
popt, pcov = curve_fit(power_law, volumes, pressures, p0=p0, bounds=bounds)
a_fit, n_fit = popt
u_a_fit, u_n_fit = np.sqrt(np.diag(pcov))
print(f"Fitresultaten: a = {a_fit:.3f} ± {u_a_fit:.3f}, n = {n_fit:.3f} ± {u_n_fit:.3f}")
# PLOT GRAFIEK MET FIT
plt.figure()
plt.xlabel('Volume ($nm^2$)')
plt.ylabel('Druk (A.U.)')
plt.plot(volumes, pressures, 'b.', label='Simulatie data')
vol_fit = np.linspace(min(volumes), max(volumes), 100)
plt.plot(vol_fit, power_law(vol_fit, *popt), 'r-', label='Power-law fit')
plt.title('Compressie: P-V diagram met fit')
plt.grid(True, alpha=0.7)
plt.legend()
plt.show()
Fitresultaten: a = 1180228.406 ± 110142.197, n = -3.403 ± 0.025
<Figure size 640x480 with 1 Axes>

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.

# Verlaag de zuigersnelheid aanzienlijk (bijv. 1% )
V_PISTON_VERIFY = -0.01 * V_0      

# Bereken het aantal stappen opnieuw voor deze snelheid
num_steps_v = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_VERIFY * DT))
print(f"Maximaal aantal stappen voor verificatie: {num_steps_v}")

# Initialiseer arrays en variabelen
particles = []
volumes_v = np.zeros(num_steps_v)
pressures_v = np.zeros(num_steps_v)
box_length = BOX_SIZE_0 
pressure = 0.0 # reset de globale druk
v_piston = V_PISTON_VERIFY

create_particles(particles)

# Run de simulatie
for i in range(num_steps_v):
    take_time_step(particles)
    volumes_v[i] = box_length * box_height
    pressures_v[i] = pressure

# Data opschonen (druk heeft tijd nodig om op te bouwen)
# We negeren de eerste 10% van de stappen omdat de 'druk-gemiddelde' functie moet stabiliseren
warmup = int(num_steps_v * 0.1)
v_clean = volumes_v[warmup:]
p_clean = pressures_v[warmup:]

# fit uitvoeren
try:
    n_guess = -2.0
    a_guess = p_clean[0] * (v_clean[0]**-n_guess) 
    popt_v, pcov_v = curve_fit(power_law, v_clean, p_clean, p0=[a_guess, n_guess], bounds=((0, -5), (np.inf, 0)))
    
    print(f"Verificatie Fitresultaten bij lage snelheid: n = {popt_v[1]:.3f} ± {np.sqrt(pcov_v[1,1]):.3f}")
    
    # Plot resultaat
    plt.figure()
    plt.plot(v_clean, p_clean, 'b.', alpha=0.3, label='Simulatie data (lage snelheid)')
    plt.plot(v_clean, power_law(v_clean, *popt_v), 'r-', label=f'Fit (n={popt_v[1]:.2f})')
    plt.xlabel('Volume ($nm^2$)')
    plt.ylabel('Druk (A.U.)')
    plt.title('Verificatie: Quasi-statische compressie')
    plt.legend()
    plt.show()
except Exception as e:
    print(f"Fit mislukt: {e}")
# resultaat geeft aan dat hoe lager de zuigersnelheid, hoe lager n wordt. En dus hoe dichter n bij het theoretische waarde van -2 ligt
Maximaal aantal stappen voor verificatie: 11111
Verificatie Fitresultaten bij lage snelheid: n = -3.584 ± 0.007
<Figure size 640x480 with 1 Axes>

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:

W=∫12PdVW = \int_1^2 P dV
(1)

Hier staan de 1 voor de begintoestand en de 2 voor de eindtoestand van het proces. In de differentiaalvorm wordt dit geschreven als:

δW=PdV\delta W = P dV
(2)

Zoals dat in het boek gebeurt, kiezen we voor de notatie met δW\delta W in plaats van dWdW, 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:

ΔW=PΔV=2DPΔA=PhvpistonΔt=ΔphΔthvpistonΔt=vpistonΔp\Delta W = P \Delta V \stackrel{2D}{=} P \Delta A = P h v_{\text{piston}} \Delta t = \frac{\Delta p}{h \Delta t} h v_{\text{piston}} \Delta t = v_{\text{piston}}\Delta p
(3)

waarin hh 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_velocity

Voordat 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.

# Instellingen
V_PISTON_0 = -0.1 * V_0
num_steps = 1000 

# Initialisatie
particles = []
create_particles(particles)
initial_energy = sum(p.kin_energy for p in particles)

# Reset de globale variabelen voor de simulatie
global work, v_piston, box_length
work = 0.0
box_length = BOX_SIZE_0
v_piston = V_PISTON_0

# arrays aan van gelijke lengte (num_steps)
cumulative_work = np.zeros(num_steps)
delta_energy_list = np.zeros(num_steps) 

# simulatie Loop
for i in range(num_steps):
    take_time_step(particles)
    cumulative_work[i] = work
    current_energy = sum(p.kin_energy for p in particles)
    delta_energy_list[i] = current_energy - initial_energy

# Plotten (nu hebben beide assen 1000 datapunten)
plt.figure(figsize=(7, 5))
plt.plot(cumulative_work, delta_energy_list, 'b-', label='Simulatie resultaat')

# voeg de ideale y=-x lijn toe voor controle
max_val = max(cumulative_work.max(), delta_energy_list.max())
plt.plot([0, max_val], [0, -max_val], 'r--', label='Theorie ($W = -\Delta U$)')

plt.xlabel('Verrichte Arbeid ($W$)')
plt.ylabel('Verandering in Energie ($\Delta E_{kin}$)')
plt.title('Validatie Eerste Hoofdwet')
plt.legend()
plt.grid(True)
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.show()
# dit is logisch want de eerste hoofdwet stelt dat de verrichte arbeid gelijk is aan de verandering in interne energie.
# de inwendige energie wordt gelijk gesteld aan de kinetische energie in dit model, met verwaarlozing van potentiële energie.
# Dit zorgt ervoor dat de verandering in kinetische energie altijd gelijk is aan de (-)verrichte arbeid.
<>:33: SyntaxWarning: invalid escape sequence '\D'
<>:36: SyntaxWarning: invalid escape sequence '\D'
<>:33: SyntaxWarning: invalid escape sequence '\D'
<>:36: SyntaxWarning: invalid escape sequence '\D'
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\2239643146.py:33: SyntaxWarning: invalid escape sequence '\D'
  plt.plot([0, max_val], [0, -max_val], 'r--', label='Theorie ($W = -\Delta U$)')
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\2239643146.py:36: SyntaxWarning: invalid escape sequence '\D'
  plt.ylabel('Verandering in Energie ($\Delta E_{kin}$)')
<Figure size 700x500 with 1 Axes>

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.

# instellingen
num_steps = 2000
half = num_steps // 2

# initialiseer lege lijst voor deeltjes
particles = []
create_particles(particles)     # resetten deeltjes 
current_initial_energy = sum(p.kin_energy for p in particles)

# reset variables
global pressure, work, box_length, v_piston
v_piston = V_PISTON_0
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0         # zetten zuiger terug

#your code/answer
# arrays van gelijke lengte (num_steps)
cumulative_work = np.zeros(num_steps)
delta_energy_list = np.zeros(num_steps)
print(num_steps//2)
#simulatie loop
for i in range(num_steps):
    if i == half:
        v_piston = -v_piston      # Richting omdraaien
        work = 0.0                # Arbeid resetten naar 0
        # Energie referentie resetten naar het warme punt
        current_initial_energy = sum(p.kin_energy for p in particles)
        
    take_time_step(particles)
    cumulative_work[i] = work
    current_energy = sum(p.kin_energy for p in particles)
    delta_energy_list[i] = current_energy - current_initial_energy

# Begint bij $(0,0)$ en gaat naar linksboven ($W$ negatief, $\Delta E$ positief).
# Begint bij $(0,0)$ en gaat naar rechtsonder ($W$ positief door de reset, $\Delta E$ negatief omdat het gas afkoelt).
# plotten van work-energy diagram voor compressie en decompressie
# compressie heeft kleur blauw, decompressie kleur oranje
plt.figure(figsize=(7, 5))
plt.plot(cumulative_work[:num_steps//2], delta_energy_list[:num_steps//2], 'r-', label='Compressie')
plt.plot(cumulative_work[num_steps//2:], delta_energy_list[num_steps//2:], 'b-', label='Decompressie')
# voeg de ideale y=-x lijn toe voor controle
max_val = max(cumulative_work.max(), delta_energy_list.max())
plt.plot([0, max_val], [0, -max_val], 'k--', label='Theorie ($W = -\Delta U$)')
plt.xlabel('Verrichte Arbeid ($W$)')
plt.ylabel('Verandering in Energie ($\Delta E_{kin}$)')
plt.title('Validatie Eerste Hoofdwet: Compressie en Decompressie')
plt.legend()
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.grid(True)
plt.show()
<>:44: SyntaxWarning: invalid escape sequence '\D'
<>:46: SyntaxWarning: invalid escape sequence '\D'
<>:44: SyntaxWarning: invalid escape sequence '\D'
<>:46: SyntaxWarning: invalid escape sequence '\D'
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\870441661.py:44: SyntaxWarning: invalid escape sequence '\D'
  plt.plot([0, max_val], [0, -max_val], 'k--', label='Theorie ($W = -\Delta U$)')
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\870441661.py:46: SyntaxWarning: invalid escape sequence '\D'
  plt.ylabel('Verandering in Energie ($\Delta E_{kin}$)')

1000
<Figure size 700x500 with 1 Axes>
# De eerste hoofdwet in de vorm $\Delta U = -W$ gaat ervan uit dat het proces quasistatisch verloopt (heel langzaam). 
# Als de zuiger extreem snel beweegt ten opzichte van de deeltjes, krijgt het gas geen tijd om een nieuw evenwicht te vinden
#your code/answer
num_steps = 400
half = num_steps // 2

# initialiseer lege lijst voor deeltjes
particles = []
create_particles(particles)     # resetten deeltjes 
current_initial_energy = sum(p.kin_energy for p in particles)

# reset variables
global pressure, work, box_length, v_piston
v_piston = 5 * V_PISTON_0
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0         # zetten zuiger terug

#your code/answer
# arrays van gelijke lengte (num_steps)
cumulative_work = np.zeros(num_steps)
delta_energy_list = np.zeros(num_steps)
print(num_steps//2)
#simulatie loop
for i in range(num_steps):
    if i == half:
        v_piston = -v_piston      # Richting omdraaien
        # current_initial_energy = current_energy
        
    take_time_step(particles)
    cumulative_work[i] = work
    current_energy = sum(p.kin_energy for p in particles)
    delta_energy_list[i] = current_energy - current_initial_energy
    
# plotten van work-energy diagram voor compressie en decompressie
# compressie heeft kleur blauw, decompressie kleur oranje
plt.figure(figsize=(7, 5))
plt.plot(cumulative_work[:num_steps//2], delta_energy_list[:num_steps//2], 'r--', label='Compressie')
plt.plot(cumulative_work[num_steps//2:], delta_energy_list[num_steps//2:], 'b-', label='Decompressie')
# voeg de ideale y=-x lijn toe voor controle
max_val = max(cumulative_work.max(), delta_energy_list.max())
plt.plot([0, max_val], [0, -max_val], 'k--', label='Theorie ($W = -\Delta U$)')
plt.xlabel('Verrichte Arbeid ($W$)')
plt.ylabel('Verandering in Energie ($\Delta E_{kin}$)')
plt.title('Validatie Eerste Hoofdwet: Compressie en Decompressie')
plt.legend()
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.grid(True)
plt.show()
# Door de hoge snelheid van de zuiger vindt het proces niet meer in quasi-evenwicht plaats, 
# waardoor de deeltjes de snelle wand tijdens decompressie niet kunnen "bijhouden" om hun energie weer volledig af te geven. 
# Hierdoor is het proces irreversibel geworden en keert het gas macroscopisch gezien niet terug naar zijn oorspronkelijke begintoestand.
<>:42: SyntaxWarning: invalid escape sequence '\D'
<>:44: SyntaxWarning: invalid escape sequence '\D'
<>:42: SyntaxWarning: invalid escape sequence '\D'
<>:44: SyntaxWarning: invalid escape sequence '\D'
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\41628895.py:42: SyntaxWarning: invalid escape sequence '\D'
  plt.plot([0, max_val], [0, -max_val], 'k--', label='Theorie ($W = -\Delta U$)')
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\41628895.py:44: SyntaxWarning: invalid escape sequence '\D'
  plt.ylabel('Verandering in Energie ($\Delta E_{kin}$)')

200
<Figure size 700x500 with 1 Axes>

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.

# Hier moet nog naar gekeken worden
# instellingen voor hoge snelheid met aangepaste DT
BOX_SIZE_0 = 10              # nm
N = 40                       # deeltjes
V_0 = 1                      # m/s
RADIUS = 0.3                 # nm
V_PISTON_FAST = 5 * (-0.1 * V_0) # -0.5 m/s
DT_FAST = (0.1 * RADIUS / V_0) / 5 # ps (DT ook gedeeld door 5 voor precisie)

num_steps = 400
half = num_steps // 2

# Initialiseer simulatie
particles = []
create_particles(particles)
initial_energy = sum(p.kin_energy for p in particles)

# Reset variabelen voor de loop
global work, box_length, v_piston, DT
work = 0.0
v_piston = V_PISTON_FAST
box_length = BOX_SIZE_0
DT = DT_FAST 
current_initial_energy = initial_energy

cumulative_work = np.zeros(num_steps)
delta_energy_list = np.zeros(num_steps)

# Simulatie loop
for i in range(num_steps):
    if i == half:
        v_piston = -v_piston  # Omdraaien naar decompressie
        work = 0.0            # Reset arbeid voor fase 2
        current_initial_energy = sum(p.kin_energy for p in particles) # Reset E-ref
        
    take_time_step(particles)
    
    cumulative_work[i] = work
    delta_energy_list[i] = sum(p.kin_energy for p in particles) - current_initial_energy

# Plotten
plt.figure(figsize=(8, 6))
plt.plot(cumulative_work[:half], delta_energy_list[:half], 'b-', label='Compressie (W<0, ΔE>0)')
plt.plot(cumulative_work[half:], delta_energy_list[half:], 'r-', label='Decompressie (W>0, ΔE<0)')

# De ideale y = -x lijn
limit = max(abs(cumulative_work).max(), abs(delta_energy_list).max())
plt.plot([-limit, limit], [limit, -limit], 'k--', alpha=0.5, label='Theorie: $\Delta U = -W$')

plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='black', lw=1)
plt.xlabel('Verrichte Arbeid ($W$) [J]')
plt.ylabel('Energieverandering ($\Delta E_{kin}$) [J]')
plt.title('Validatie 1e Hoofdwet: Hoge Snelheid met DT-correctie')
plt.legend()
plt.grid(True, linestyle=':')
plt.axis('equal') 
plt.show()
<>:48: SyntaxWarning: invalid escape sequence '\D'
<>:53: SyntaxWarning: invalid escape sequence '\D'
<>:48: SyntaxWarning: invalid escape sequence '\D'
<>:53: SyntaxWarning: invalid escape sequence '\D'
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\1685624025.py:48: SyntaxWarning: invalid escape sequence '\D'
  plt.plot([-limit, limit], [limit, -limit], 'k--', alpha=0.5, label='Theorie: $\Delta U = -W$')
C:\Users\Eline\AppData\Local\Temp\ipykernel_17088\1685624025.py:53: SyntaxWarning: invalid escape sequence '\D'
  plt.ylabel('Energieverandering ($\Delta E_{kin}$) [J]')
<Figure size 800x600 with 1 Axes>
# opdrachrt met constant volume
#Instellingen en SI-Conversie factoren
# In de simulatie is massa m=1 en V_0=1.
m = 6.64e-27          # kg
v_unit = 1000.0          # 1 unit snelheid = 1000 m/s
e_unit = m * v_unit**2 # Conversie factor voor Energie naar Joules

# Simulatie parameters
BOX_SIZE_0 = 10.0        # nm
V_PISTON_CONST = 0.2     # Snelheid van beide zuigers in dezelfde richting
DT = 0.03                # ps
num_steps = 1000

# functies
def take_time_step_same_direction(particles, left_wall, right_wall, v_piston):
    """ Verplaatst beide wanden in dezelfde richting. """
    new_left = left_wall + v_piston * DT
    new_right = right_wall + v_piston * DT
    
    for p in particles:
        p.update_position()
        
        # Botsing linkerzuiger (beweegt met v_piston)
        if p.r[0] - p.R < new_left:
            p.r[0] = new_left + p.R
            p.v[0] = 2 * v_piston - p.v[0]
            
        # Botsing rechterzuiger (beweegt ook met v_piston)
        if p.r[0] + p.R > new_right:
            p.r[0] = new_right - p.R
            p.v[0] = 2 * v_piston - p.v[0]
            
        # Stilstaande boven/onderwanden
        if abs(p.r[1]) + p.R > BOX_SIZE_0 / 2:
            p.r[1] = np.sign(p.r[1]) * (BOX_SIZE_0/2 - p.R)
            p.v[1] *= -1
            
    return new_left, new_right

# Run Simulatie
particles = []
create_particles(particles) # Gebruikt de functie van eerder

left_wall = -BOX_SIZE_0 / 2
right_wall = BOX_SIZE_0 / 2
temps = []
displacements = []

for i in range(num_steps):
    left_wall, right_wall = take_time_step_same_direction(particles, left_wall, right_wall, V_PISTON_CONST)
    
    # Bereken T in SI (Kelvin)
    # E_kin_avg = 1/2 * m * v^2 = 3/2 * kB * T 
    # We gebruiken hier de 2D relatie: T = E_kin_totaal / (N * kB)
    kB = 1.38e-23
    total_energy_si = sum(p.kin_energy for p in particles) * e_unit
    temp_si = total_energy_si / (N * kB)
    
    temps.append(temp_si)
    displacements.append(i * V_PISTON_CONST * DT)

# Plotten in SI-eenheden
plt.figure(figsize=(8, 5))
plt.plot(np.array(displacements), temps, 'g-')
plt.xlabel('Verplaatsing van het volume (nm)')
plt.ylabel('Temperatuur (K)')
plt.title('Temperatuur bij verplaatsing met constante inhoud ($dV=0$)')
plt.grid(True, linestyle=':')
plt.show()
<Figure size 800x500 with 1 Axes>
VSC
De Maxwell-Boltzmann snelheidsverdeling
VSC
Thermostaat