Arbeid - Thermolab

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.

# Ex 1

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 (nm)
N = 40                          # Aantal deeltjes
V_0 = 1                         # Startsnelheid van deeltjes (nm per ps)
RADIUS = 0.3                    # Straal van moleculen (nm)
DT = (0.1 * RADIUS) / V_0       # Tijdstap om geen botsing te missen (s)
V_PISTON_0 = -0.1 * V_0         # Startsnelheid van zuiger 

print(V_PISTON_0)
# (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.

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.

# Ex 2

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

# Addition of the function for temperature, made in 5_druktemp
# Correction for the temperature
# ((Temperature * mass_factor * ((distance_factor/(time_factor)**2) * ConstantBoltzmann_factor) / length of array)
def temperature(particles) -> float:
    temp = 0.0
    for p in particles:
        temp += (p.kin_energy)
    return ((temp * 4.8107e-26 * ((1e-9/(2.5 * 1e-12))**2) * 1.380649e23) / len(particles))
        
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()

Nu implementeren we de zuiger door het toegestane gebied voor de gasdeeltjes bij elke tijdstap te verkleinen met een stap $2v_{\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())

# 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 $v_{\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:

# Ex 3

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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)


plt.plot(volumes, temperatures, '-r')

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.

# Ex 4

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         # zetten zuiger terug
v_piston = V_PISTON_0
create_particles(particles)     # resetten deeltjes 

# Addition of the factor for the pressure
pressure_factor = ((4.8107e-26 * (1e-9/(1e-12))) / (1e-9 * 1e-12))
# Pressure_factor = (mass_factor * velocity_factor) / (length_factor * time_factor) 

for i in range(num_steps):
    take_time_step(particles)
    volumes[i] = box_length * box_height
    pressures[i] = (pressure * pressure_factor)

# plot
plt.figure()
plt.xlabel(r'Volume (nm$^2$)')
plt.ylabel('Pressure (Bar)')

pressures = np.asarray(pressures)


plt.plot(volumes, pressures, '-r')

plt.show()

Als je simulatie klopt heeft deze grafiek de vorm van een machtsfunctie met een negatieve exponent.

Solution to Exercise 5

Uitgaande van de ideale gaswet, geldt het volgende: $pV = Nk_{B}T$ .

Hierin is

$p$ de druk,
$V$ het volume,
$N$ het aantal deeltjes,
$k_{B}$ de constante van Boltzmann,
$T$ de temperatuur

In de simulatie is er sprake van een bewegende zuiger, die het volume verandert. Hierdoor wordt er arbeid verricht en als gevolg verandert de kinetische energie van de deeltjes en dus vind er ook een verandering in de temperatuur plaats.
Daarbij is er geen sprake van warmteuitwisseling met de omgeving. Dat betekent dat elke volumeverkleining zorgt voor een volledige omzetting in toename van de kinetische energie en dus de temperatuur.
Met dit in gedachte, kan worden gezegd dat de temperatuur een omgekeerd rechtevenredig verband heeft met het volume, $T \propto V^{-1}$ .

Daarnaast stelt de ideale gaswet het volgende: $pV = Nk_{B}T$ .

Als vervolgens deze ideale gaswet wordt samengevoegd met het gevonden omgekeerd rechtevenredig verband, wordt het volgende verkregen.

pV = Nk_{B}\frac{1}{V}

(1)

Aangezien $N$ en $k_{B}$ een constante waarde hebben, worden deze versimpeld naar $Nk_{B} = a$ .
Wanneer vervolgens de $p$ wordt vrijgemaakt in de vergelijking, levert dit de volgende vergelijking.

p = \frac{a}{V^{2}}

(2)

Deze vergelijking toont dat de verwachte exponent -2 zal zijn.

# Ex 6

def power_law(vol, a, n):
    """ de fitfunctie voor het P,V-diagram """
    return a * (vol)**n   

#popt_pv, pcov_pv = curve_fit(power_law, volumes, pressures, p0=(2, -1), bounds=((-100, -10),(5000,0)))
popt_pv, pcov_pv = curve_fit(power_law, volumes, pressures, p0=(2, -1))

# Variables
a = popt_pv[0]
a_u = np.sqrt(pcov_pv[0,0])
n = popt_pv[1]
n_u = np.sqrt(pcov_pv[1,1])

print('De variabel a heeft een waarde van (%.1f \u00B1 %.1f).' %(a, a_u))
print('De variabel n heeft een waarde van (%.1f \u00B1 %.1f).' %(n, n_u))

# Calculating the minamum and maximum value of the variable volumes
min_value_v = np.min(volumes)
max_value_v = np.max(volumes)

x_pv = np.linspace(0.95*min_value_v, 1.05*max_value_v, 1000)
y_pv = power_law(x_pv, *popt_pv)

#Plotting the curve fit and the data
plt.figure()
plt.plot(volumes, pressures, 'r-', label='Ex 4')
plt.plot(x_pv, y_pv, 'b.', label='Ex 6')
plt.xlabel(r'Volume (nm$^2$)')
plt.ylabel('Pressure (Bar)')
plt.legend()
plt.show()

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.

Ruimte voorverificatie¶

#your code/answer

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 = \int_1^2 P dV

(3)

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

\delta W = P dV

(4)

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

\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

(5)

waarin $h$ 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.

Solution to Exercise 8

In de simulatie gaat het om een gas in een piston. Hierbij moet voldaan worden aan de wet van behoud van energie.

\Delta E = Q - W

(6)

Echter, is bekend dat er geen warmte wordt uitgewissels, waardoor de term $Q$ wegvalt.
Daarbij kan de term $\Delta E$ worden omgeschreven naar de kinetische energie van de deeltjes, $\Delta KE$ , aangezien de arbeid volledig wordt omgezet in kinetische energie, zoals uitgelegd bij opgave 5.

Kortom, er geldt nu

\Delta KE = - W

(7)

Om te valideren of de simulatie voldoet aan de eerste hoofdwet van Thermodynamica, kan de grootheid kinetische energie, $\Delta KE$ worden uitgezet tegen de grootheid arbeid, $W$ . De arbeid kan berekent worden door de functie bepaald onder het kopje ‘Eerste Hoofdwet’.

Er wordt dan verwacht dat er een grafiek zal onstaan met een helling van -1.

# Ex 9

num_steps = 1000
particles = []
pressure = 0.0
work = 0.0

box_length = BOX_SIZE_0         # zetten zuiger terug
v_piston = V_PISTON_0
create_particles(particles)     # resetten deeltjes 

W = np.zeros(num_steps, dtype=float)
KE = np.zeros(num_steps, dtype=float)

KE0 = sum(p.kin_energy for p in particles)
for i in range(num_steps):
    take_time_step(particles)
    W[i] = work                      
    KE[i] = sum(p.kin_energy for p in particles) - KE0
    
#work = np.asarray(work)

# plot
plt.figure()
plt.xlabel(r'Work (J)')
plt.ylabel(r'$\Delta KE$ (J)')

plt.plot(W, KE, '-r', label=r'$\Delta KE = -W$')
plt.legend()
plt.show()

print('Uit het resultaat blijkt dat inderdaad de hypothese klopt.')

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.

# Ex 10

num_steps_1 = 1000
num_steps_2 = 1000

particles = []
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0         # zetten zuiger terug
v_piston = V_PISTON_0
create_particles(particles)     # resetten deeltjes

W_1 = np.zeros(num_steps_1, dtype=float)
KE_1 = np.zeros(num_steps_1, dtype=float)
W_2 = np.zeros(num_steps_2, dtype=float)
KE_2 = np.zeros(num_steps_2, dtype=float)

KE0 = sum(p.kin_energy for p in particles)

# Seperation for the two timesteps
for i in range(num_steps_1):
    take_time_step(particles)
    W_1[i] = work                      
    KE_1[i] = sum(p.kin_energy for p in particles) - KE0

for i in range(num_steps_2):
    take_time_step(particles)
    v_piston = -V_PISTON_0
    W_2[i] = work                      
    KE_2[i] = sum(p.kin_energy for p in particles) - KE0

# plot
plt.figure()
plt.xlabel(r'Work (J)')
plt.ylabel(r'$\Delta KE$ (J)')

plt.plot(W_1, KE_1, '-r', label=r'$\Delta KE = -W$')
plt.plot(W_2, KE_2, 'k:', label='Reversed')
plt.title('Exercise 10, Work against kinetic energy')
plt.legend()
plt.show()

timesteps1 = np.linspace(1, 1000, 1000)
timesteps2 = np.linspace(1001, 2000, 1000)

# Another plot
plt.figure()
plt.xlabel('Timesteps')
plt.ylabel('Work (J)')
plt.plot(timesteps1, W_1, '-r', label='Normal')
plt.plot(timesteps2, W_2, 'k:', label='Reversed')
plt.title('Exercise 10, Work against timesteps')
plt.legend()
plt.show()

# Ex 11

num_steps_1_5 = 200
num_steps_2_5 = 200

particles = []
pressure = 0.0
work = 0.0
box_length = BOX_SIZE_0         # zetten zuiger terug
v_piston = 5 * V_PISTON_0
create_particles(particles)     # resetten deeltjes

W_1_5 = np.zeros(num_steps_1_5, dtype=float)
KE_1_5 = np.zeros(num_steps_1_5, dtype=float)
W_2_5 = np.zeros(num_steps_2_5, dtype=float)
KE_2_5 = np.zeros(num_steps_2_5, dtype=float)

KE0_5 = sum(p.kin_energy for p in particles)

# Seperation for the two timesteps
for i in range(num_steps_1_5):
    take_time_step(particles)
    W_1_5[i] = work                      
    KE_1_5[i] = sum(p.kin_energy for p in particles) - KE0_5

for i in range(num_steps_2_5):
    take_time_step(particles)
    v_piston = -5 * V_PISTON_0
    W_2_5[i] = work                      
    KE_2_5[i] = sum(p.kin_energy for p in particles) - KE0_5

# plot
plt.figure()
plt.xlabel(r'Work (J)')
plt.ylabel(r'$\Delta KE$ (J)')

plt.plot(W_1_5, KE_1_5, '-r', label=r'$\Delta KE = -W$')
plt.plot(W_2_5, KE_2_5, 'k:', label='Reversed')
plt.title('Exercise 11, Work against kinetic energy')
plt.legend()
plt.show()

# Another Plot
timesteps1_5 = np.linspace(1, 200, 200)
timesteps2_5 = np.linspace(201, 400, 200)

# Another plot
plt.figure()
plt.xlabel('Timesteps')
plt.ylabel('Work (J)')
plt.plot(timesteps1_5, W_1_5, '-r', label='Normal')
plt.plot(timesteps2_5, W_2_5, 'k:', label='Reversed')
plt.title('Exercise 11, Work against timesteps')
plt.legend()
plt.show()

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.

# Ex 12

# Normal velocity
time1 = np.abs(timesteps1 * V_PISTON_0)
time2 = np.abs(timesteps2 * -V_PISTON_0)

plt.figure()
plt.xlabel('Time (ps)')
plt.ylabel('Work (J)')
plt.plot(time1, W_1, '-r', label='Normal')
plt.plot(time2, W_2, 'k:', label='Reversed')
plt.title('Exercise 10, Work against timesteps')
plt.legend()
plt.show()

# Velocity times 5
time1_5 = np.abs(timesteps1_5 * (5* V_PISTON_0))
time2_5 = np.abs(timesteps2_5 * (5* -V_PISTON_0))

plt.figure()
plt.xlabel('Time (ps)')
plt.ylabel('Work (J)')
plt.plot(time1_5, W_1_5, '-r', label='Normal')
plt.plot(time2_5, W_2_5, 'k:', label='Reversed')
plt.title('Exercise 11, Work against timesteps')
plt.legend()
plt.show()

# Ex 14

# Take_time_step functie veranderen, zodat de twee zuigers in dezelfde richting bewegen 
def take_time_step(particles):
    """ zet tijdstap voor een lijst deeltjes en verwerk alle botsingen onderling en met wanden """
    global box_length, v_piston, zuiger, speed_factor
    #box_length += 2 * v_piston * DT # zowel links als rechts zuiger
    zuiger +=  v_piston * DT * speed_factor         #Nieuwe code voor de beweging van de zuiger
    for p in particles:
        p.update_position()
    handle_walls(particles)  
    handle_collisions(particles)

# Create_particles functie veranderen, zodat de x-coordinaat rekening houdt met het feit dat de zuigers bewegen
def create_particles(particles):
    """ Leegmaken en opnieuw aanmaken van deeltjes  in lijst """
    global box_length, box_height, zuiger
    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+zuiger, box_length/2 - RADIUS+zuiger)
        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))

# Left_right_collision aanpassen, zodat er rekening wordt gehouden met de zuigers die dezelfde kant op bewegen
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, zuiger

    if particle.r[0] + particle.R > box_length / 2 + zuiger:
        particle.r[0] = np.sign(particle.r[0]) * (box_length/2 + zuiger - 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 += abs(particle.momentum[0]) * 2
        impulse_outward += 2 * particle.m * abs(relative_velocity)
        work += 2 * particle.m * relative_velocity * piston_velocity
    if particle.r[0] + particle.R < -(box_length / 2) + zuiger:
        particle.r[0] = (-box_length/2 + zuiger - 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 += abs(particle.momentum[0]) * 2
        impulse_outward += 2 * particle.m * abs(relative_velocity)
        work += 2 * particle.m * relative_velocity * piston_velocity

speed_factor = 2
zuiger = 2

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

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


# Num_steps worden bepaald, zodat de zuigers minimaal over de helft van de lengte van het volume bewegen
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 - 4         # zetten zuiger terug
v_piston = V_PISTON_0
zuiger = 2
speed_factor = 2

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

# Plots of the temperature with different speed_factors

# First a speed_factor equal to 0.5
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 0.5
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = \frac{1}{2}v_{\text{Piston}}$")
plt.title('Speed factor equal to 0.5')
plt.legend()


# Secondly with a speed_factor equal to 1
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 1
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = v_{\text{Piston}}$")
plt.title('Speed factor equal to 1')
plt.legend()


#Thirdly with a speed_factor equal to 1.5
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 1.5
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = \frac{3}{2}v_{\text{Piston}}$")
plt.title('Speed factor equal to 1.5')
plt.legend()
plt.show()


# Fourtly a speed_factor equal to 2
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 2
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = 2v_{\text{Piston}}$")
plt.title('Speed factor equal to 2')
plt.legend()
plt.show()


# Last of all a speed_factor equal to 4
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 4
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label=r"$v = 4v_{\text{Piston}}$")
plt.legend()
plt.title('Speed factor equal to 4')
plt.show()

Wanneer de speed factor wordt aangepast, oftewel de snelheid van de zuiger wordt verandert, dan heeft dit invloed op de temperatuur.
Zoals te zien is in bovenstaande plots, wordt de temperatuur hoger naarmate de snelheid van de zuigers groter wordt.

# Plots of the temperature with different start temperatures

# First with a start temperature of 100 K
def temperature(particles) -> float:
    temp = 100.0
    for p in particles:
        temp += (p.kin_energy)
    return ((temp * 4.8107e-26 * ((1e-9/(2.5 * 1e-12))**2) * 1.380649e23) / len(particles))

# Now with a start temperature of 100 K, firstly a speed_factor equal to 0.5 is chosen
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 0.5
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = \frac{1}{2}v_{\text{Piston}}$")
plt.title('Start temperature of 100K with a speed factor equal to 0.5')
plt.legend()
plt.show()

# Still with a start temperature of 100 K, secondly a speed_factor equal to 2 is chosen
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 2
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = 2v_{\text{Piston}}$")
plt.title('Start temperature of 100K with a speed factor equal to 2')
plt.legend()
plt.show()


# Now still with a start temperature of 100K, a speed_factor of 4 is chosen
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 4
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label=r"$v = 4v_{\text{Piston}}$")
plt.legend()
plt.title('Start temperature of 100K with a speed factor equal to 4')
plt.show()



# Secondly, the start temperature is changed to 1000 K
def temperature(particles) -> float:
    temp = 10000.0
    for p in particles:
        temp += (p.kin_energy)
    return ((temp * 4.8107e-26 * ((1e-9/(2.5 * 1e-12))**2) * 1.380649e23) / len(particles))

# Now with a start temperature of 1000 K, firstly a speed_factor equal to 0.5 is chosen
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 0.5
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = \frac{1}{2}v_{\text{Piston}}$")
plt.title('Start temperature of 10000K with a speed factor equal to 0.5')
plt.legend()
plt.show()

# Still with a start temperature of 1000 K, secondly a speed_factor equal to 2 is chosen
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 2
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label = r"$v = 2v_{\text{Piston}}$")
plt.title('Start temperature of 10000K with a speed factor equal to 2')
plt.legend()
plt.show()

# Now still with a start temperature of 1000K, a speed_factor of 4 is chosen
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
speed_factor = 4
zuiger = 2
box_length = BOX_SIZE_0 - 4         # 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(r'Volume (nm$^2$)')
plt.ylabel('Temperature (K)')

temperatures = np.asarray(temperatures)

plt.plot(volumes, temperatures, '-r', label=r"$v = 4v_{\text{Piston}}$")
plt.legend()
plt.title('Start temperature of 10000K with a speed factor equal to 4')
plt.show()

Wanneer de start temperatuur wordt aangepast, dan heeft dit invloed op de temperatuur.
Zoals te zien is in bovenstaande plots, wordt de temperatuur hoger naarmate de start temperatuur ook hoger is.
Echter, moet de start temperatuur redelijk erg verhoogt worden om een duidelijk verschil te zien in de plots. Vandaar dat er gekozen is voor een start temperatuur van 100K en een start temperatuur van 10000K, oftewel 100x zo hoog. Met deze verhoging is er bij een snelheid van 4 maal de piston snelheid een verhoging van ongeveer 1,2x te zien in de plots.
Daarbij lijkt het erop dat een hogere start temperatuur meer invloed heeft bij een lagere snelheid van de zuigers, aangezien bij de laagste snelheid van de zuigers de temperatuur het meest verschilt tussen de twee plots.

Verklaring Gedrag
Er is gekeken naar de invloed van de snelheid van de zuigers op de temperatuur en naar de invloed van de start temperatuur op de temperatuur. Een verklaring van het gedrag wordt zowel op macroscopisch als op microscopisch niveau genoemd.

Microscopisch niveau
Ten eerste zorgt een hogere snelheid van de zuigers op microscopisch niveau ervoor dat de deeltjes harder tegen de wand kunnen botsen en vervolgens ook harder tegen andere deeltjes. Deze botsingen zullen dan zorgen voor een hogere kinetische energie. Aangezien de temperatuur wordt bepaald door middel van de kinetische energie met behulp van onderstaande vergelijking.

T = \frac{1}{2} \frac{m v^{2}}{k_{B}} = \frac{E_{Kin}}{k_{B}}

(8)

Kortom, als de kinetische energie, $E_{Kin}$ , toeneemt, zal ook de temperatuur, $T$ , toenemen, aangezien de constante van boltzmann, $k_{B}$ , niet zal veranderen. Dus een hogere snelheid van de zuigers zorgt voor een hogere temperatuur.

Ten tweede zorgt een hogere start temperatuur van de simulatie op microscopisch niveau ervoor dat de deeltjes sneller gaan bewegen en daardoor harder op elkaar of de wanden kunnen botsen. Ook hiervoor geldt dan vervolgens dat dit resulteert in een hogere kinetische energie, wat als gevolg heeft dat de eind temperatuur ook hoger is, wegens dezelfde reden als eerder uitgelegd.
Kortom, een hogere start temperatuur van de simulatie zorgt voor een hogere eind temperatuur.

Macroscopisch niveau

Om een verklaring te geven op macroscopisch niveau kunnen we gebruik maken van klasieke mechanica. Wanneer de zuigers allebij met dezelfde constante snelheid en zelfde richting bewegen veranderd het volume niet. Dit volume kan je dus zien als een doosje dat je voorduwt in de ruimte met een bepaalde snelheid. Dat is hetzelfde als aangeven dat het doosje stilstaat en de ruimte eromheen beweegt in tegengestelde snelheid. En omdat het volume geisoleerd is van de buitenomgeving zouden we moeten verwachten dat de temperatuur niet stijgt. Immers veranderd het volume niet en wordt er geen arbeid aan het systeem toegevoegd. Waarom zien we dan in onze simulatie een verandering van temperatuur? Als startvoorwaarde hebben we de deeltjes in het gas allemaal een start snelheid gegeven. Dit is alleen relatief tenopzichte van een stilstaand volume. Wanneer het doosje beweegt maar je nogsteeds kijkt vanaf het referentie systeem van de ruimte krijg je foutieve resultaten zoals we nu zien in vorrige grafieken en de verklaring op microscopisch niveau. Wanneer we opnieuw een N-aantal deeltjes zouden aanmaken allemaal met een snelheid Vx+zuiger_snelheid, en we daarna de simulatie zouden draaien blijkt dat de temperatuur wel constant zou blijven. Zoals we ook zouden verwachten.