Entropie 🌶 - Thermolab

Doelen¶

In hoofdstukken 5 en 6 van het boek wordt er verder ingegaan op de reversibiliteit van processen. De klassieke mechanica eist alleen behoud van energie en impuls, maar er is geen richtingsvoorkeur in processen. Dat is in strijd met je dagelijkse ervaring: van een video-opname kun je vrijwel altijd zeggen of deze vooruit of achteruit wordt gespeeld.

Deze richtingsvoorkeur wordt het best gemodelleerd door de entropie. Dit is een kwantitatieve grootheid die vaak een beetje mysterieus wordt gevonden. In dit werkblad zullen we kijken of we wat meer over deze entropie te weten kunnen komen.

In dit werkblad komen de volgende onderdelen langs:

We herhalen de belangrijke stukken code met gekozen constanten en het gedrag van deeltjes.
We introduceren een nieuwe klasse voor het beheersvolume die de reeds ontwikkelde functies bevat.
We controleren of deze code gelijkwaardige resultaten geeft als de code in voorgaande werkbladen.
We bekijken een ingewikkelder modelsysteem van gekoppelde zuigers en kijken hier naar verschillende processen.
We nemen een reversibele en een niet-reversibele route en kijken naar het verschil in entropie.

Laden van eerdere code¶

Net als bij de eerdere simulaties, gaan we uit van dezelfde set van constanten.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

FAST_MODE = True  # Zet op False voor de volledige (langzamere) simulatie
FAST_MAX_STEPS = 6000  # cap voor lange loops (snelle modus)


# ----------------------------
# Constanten (gekozen zodat T ≈ 300 K en diameter ≈ 0.1 nm)
# ----------------------------
k_B = 1.38e-23                 # J/K
MASS = 4.65e-26                # kg (orde van grootte luchtmolecuul)
T0 = 300                       # K

BOX_SIZE_0 = 20e-9             # m  (hoogte en lengte startvolume)
N = 30 if FAST_MODE else 40
RADIUS = 0.05e-9               # m  (diameter 0.1 nm)

# kies V_0 zodat de kinetische energie overeenkomt met T0 (2D: <E_kin> = k_B T)
V_0 = np.sqrt(2 * k_B * T0 / MASS)

# tijdstap
DT = 1.0e-13                   # s

# startsnelheid van zuiger (negatief betekent beide zuigers naar binnen)
V_PISTON_0 = (-0.08 if FAST_MODE else -0.02) * V_0
# reproduceerbaarheid
np.random.seed(0)

En maken we daarnaast ook gebruik van een klasse die het gedrag van deeltjes omschrijft:

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 beheersvolume¶

In het boek wordt er vaak gesproken over een ‘control volume’. Dit is een reëel of denkbeeldig volume waaraan een aantal macroscopische eigenschappen worden toegekend. Het Nederlandse woord hiervoor is ‘beheersvolume’. Veel van de functies die we tot nu toe hebben ontwikkeld gaan over eigenschappen en gedrag van een dergelijk beheersvolume - deze hebben we slechts impliciet gebruikt. Om de code beter te structureren, maken we hier een klasse voor. Dit maakt het vooral makkelijk om bijvoorbeeld met twee beheersvolumes te rekenen (later in dit werkblad is dat ons doel).

Als je de code hieronder vergelijkt met die van de vorige werkbladen, zal je een aantal zaken opvallen:

Elke functie heeft een parameter self gekregen. Daarmee begrijpt Python dat het om een functie van de klasse gaat.
De variabelen die als global stonden vermeld zijn nu een variabele van de klasse en hoeven dus niet meer gedeclareerd te worden in elke functie. Ze moeten in de code wel telkens worden voorafgegaan door het woord self, om dit duidelijk te maken.
Om de eigenschappen heat, work en pressure alleen te laten lezen en niet te laten schrijven, maken we gebruik van een Python conventie waarbij we een ‘geheime’ variabele gebruiken die wordt voorafgegaan door een laag liggend streepje ‘_’ (Engels: underscore).

class ControlVolume:
    def __init__(self, length, height, N, v_piston, set_temp):
        """ maakt een beheersvolume (constructor) """
        self.length = length
        self.height = height 
        self.v_piston = v_piston
        self.set_temp = set_temp
        self.particles = []
        self._work = 0.0
        self._heat = 0.0
        self._impulse_out = 0.0
        self._pressure = 0.0
        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(-self.length/2 + RADIUS, self.length/2 - RADIUS)
            y = np.random.uniform(-self.height/2 + RADIUS, self.height/2 - RADIUS)
            self.particles.append(ParticleClass(m=MASS, v=[vx, vy], r=[x, y], R=RADIUS))

    def handle_collisions(self) -> None:
        """ alle onderlinge botsingen afhandelen voor deeltjes in lijst """
        num_particles = len(self.particles)
        for i in range(num_particles):
            for j in range(i+1, num_particles):
                if collide_detection(self.particles[i], self.particles[j]):
                    particle_collision(self.particles[i], self.particles[j])

    def piston_collision(self, particle: ParticleClass) -> None:
        """ verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
        if abs(particle.r[0]) + particle.R > self.length / 2:
            particle.r[0] = np.sign(particle.r[0]) * (self.length/2 - particle.R)
            piston_velocity = np.sign(particle.r[0]) * self.v_piston
            relative_velocity = particle.v[0] - piston_velocity  # stelsel zuiger
            particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
            self._impulse_out += 2 * particle.m * abs(relative_velocity)
            self._work += 2 * particle.m * relative_velocity * piston_velocity

    def thermostat_collision(self, particle: ParticleClass) -> None:
        """ verzorgen van botsingen met wand boven en onder, die als thermostaat kunnen werken """
        if abs(particle.r[1]) + particle.R > self.height / 2:
            temp_factor = (self.set_temp/self.temperature) if self.set_temp > 0 else 1.0
            particle.r[1] = np.sign(particle.r[1]) * (self.height/2 - particle.R)
            self._impulse_out += abs(particle.momentum[1]) * (1 + temp_factor**0.5) 
            self._heat += particle.kin_energy * (temp_factor - 1)
            particle.v *= temp_factor**0.5
            particle.v[1] *= -1

    def handle_walls(self) -> None:
        """ verzorgen van alle botsingen met wanden en aanpassen waarde voor druk """
        self._impulse_out = 0.0
        for p in self.particles:
            self.piston_collision(p)
            self.thermostat_collision(p)
        self._pressure = 0.95 * self._pressure + 0.05 * self._impulse_out / (self.circumference * DT)

    def take_time_step(self) -> None:
        self.length += 2 * self.v_piston * DT # zowel links als rechts zuiger
        for p in self.particles:
            p.update_position()
        self.handle_collisions()
        self.handle_walls()  

    @property
    def temperature(self) -> float:
        # In 2D geldt per deeltje: <E_kin> = k_B T
        if len(self.particles) == 0:
            return 0.0
        total_ke = 0.0
        for p in self.particles:
            total_ke += p.kin_energy
        temp = total_ke / (len(self.particles) * k_B)
        return temp

    @property
    def area(self) -> float:
        return self.length * self.height
    
    @property 
    def circumference(self) -> float:
        return 2 * (self.length + self.height)
    
    @property
    def heat(self) -> float:
        return self._heat
    
    @property
    def work(self) -> float:
        return self._work
    
    @property
    def pressure(self) -> float:
        return self._pressure

Herhaling oude test¶

Om zeker te zijn dat de code functioneert herhalen we eerst een berekening, waarvan we het antwoord al kennen. Dit is de simulatie van de druk en de temperatuur als functie van het volume voor een isotherm proces (waarbij de temperatuur dus constant wordt gehouden met behulp van een thermostaat).

num_steps = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))

volumes = np.zeros(num_steps, dtype=float)
pressures = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)

cv = ControlVolume(BOX_SIZE_0, BOX_SIZE_0, 40, V_PISTON_0, 300)

for i in range(num_steps):
    cv.take_time_step()
    volumes[i] = cv.area
    pressures[i] = cv.pressure
    temperatures[i] = cv.temperature

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
ax1.set_xlabel('Volume')
ax1.set_ylabel('Pressure')
ax2.set_xlabel('Volume')
ax2.set_ylabel('Temperature')
fig.tight_layout

ax1.plot(volumes, pressures, '-r')
ax2.plot(volumes, temperatures, '-b')

plt.show()

Systeem van twee gekoppelde zuigers¶

Om nu een studie te doen naar de entropie van een systeem maken we gebruik van twee gekoppelde zuigers, zie Figure 1. Dit zijn twee zuigers genummerd ‘1’ en ‘2’ die samen een constant volume hebben, maar waarbij de beweging van de zuiger(s) het ene volume verkleint en het andere volume evenveel vergroot. Dit systeem houden we voor de rest van dit weerkblad thermisch geïsoleerd ten opzichte van de omgeving, zodat we een goede boekhouding kunnen doen van de totale hoeveelheid energie.

Figure 1:De gekoppelde zuigers, waarbij het totaal volume constant is.

Isotherme verplaatsing¶

Als eerste kijken we naar het reversibele proces. We herhalen dat er geen thermisch contact is met de omgeving, maar er is wel thermisch contact tussen de twee volumes zodat ze constant dezelfde temperatuur hebben.

In de onderstaande code beginnen we met twee beheersvolumes cv1 en cv2, die aan het begin hetzelfde volume $V_0$ hebben. Dan verplaatsen we de zuiger(s) met een constante snelheid naar het punt waar $V_1=\frac{1}{5}V_0$ en $V_2= \frac{9}{5} V_0$ (wat verwacht je, conceptueel, wat er gebeurt met de verschillende grootheden?)

Hieronder vind je het verloop van de druk en de temperatuur van beide volumes tijdens dit proces.

# runnen van de simulatie
num_data = round(abs(BOX_SIZE_0 * 0.4 / (DT * V_PISTON_0)))
if FAST_MODE:
    num_data = min(num_data, FAST_MAX_STEPS)

volumes1 = np.zeros(num_data, dtype=float)
pressures1 = np.zeros(num_data, dtype=float)
temperatures1 = np.zeros(num_data, dtype=float)
volumes2 = np.zeros(num_data, dtype=float)
pressures2 = np.zeros(num_data, dtype=float)
temperatures2 = np.zeros(num_data, dtype=float)

cv1 = ControlVolume(BOX_SIZE_0, BOX_SIZE_0, (30 if FAST_MODE else 50), +V_PISTON_0, 0.5)
cv2 = ControlVolume(BOX_SIZE_0, BOX_SIZE_0, (30 if FAST_MODE else 50), -V_PISTON_0, 0.5)

for i in range(num_data):
    cv1.take_time_step()
    cv2.take_time_step()
    
    average_temp = (cv1.temperature + cv2.temperature) / 2
    cv1.set_temp = average_temp
    cv2.set_temp = average_temp
    
    volumes1[i] = cv1.area
    pressures1[i] = cv1.pressure
    temperatures1[i] = cv1.temperature
    
    volumes2[i] = cv2.area
    pressures2[i] = cv2.pressure
    temperatures2[i] = cv2.temperature

# plotten van de druk en temperatuur van beide volumes
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3))
ax1.set_xlabel('Volume')
ax1.set_ylabel('Pressure')

ax2.set_xlabel('Volume')
ax2.set_ylabel('Temperature')
fig.tight_layout

#druk als functie van volume
ax1.plot(volumes1, pressures1, '-r', label='cv1')
ax1.plot(volumes2, pressures2, '-m',label='cv2')
ax1.legend()

#temperatuur als functie van volume
ax2.plot(volumes1, temperatures1, '-b',label='cv1')
ax2.plot(volumes2, temperatures2, '-g',label='cv2')
ax2.legend()

plt.show()

{solution} ex-entropy-04 Als de zuiger uitwijkt, wordt het ene volume kleiner en het andere groter. Omdat $P \propto T/V$ stijgt de druk aan de kleine kant; om de zuiger verder te verplaatsen moet er netto arbeid op het totale systeem worden verricht.
Die arbeid gaat naar de interne energie, dus $U$ stijgt en daarmee ook $T$ .

Voor grotere uitwijking wordt het volume van de kleine kant heel klein, waardoor $P$ sterk toeneemt en de benodigde arbeid per extra stap groter wordt. Dat zie je ook in de afleiding:

\frac{dT}{dx}=\frac{x}{1-x^2}T,

(1)

waarbij de factor $\frac{1}{1-x^2}$ snel groter wordt als $|x|\to 1$ .

Dit systeem is zo gekozen dat het mogelijk is om de temperatuur analytisch te bepalen met behulp van een wiskundige afleiding. We gaan ervan uit dat het systeem bij de start helemaal symmetrisch is (zoals in de simulatie) en definiëren hierbij $x$ als de positie van de zuiger. Bij $x=-1$ staat de zuiger helemaal aan de kant van cv1, zodat $V_1=0$ en $V_2=2V_0$ . Bij $x=+1$ staat de zuiger helemaal aan de andere kant, zodat $V_1=2V_0$ en $V_2=0$ . Dan moet dus gelden dat:

V_1 = V_0 (1+x) \text{ en } V_2 = V_0 (1-x).

(2)

Omdat ons gas twee vrijheidsgraden heeft, wordt de interne energie van elk van de volumes ( $i \in {1,2}$ ) gegeven door:

U_i = \frac{2}{2} N_i k_B T.

(3)

Met het aantal deeltjes aan beide kanten gelijk betekent dit voor de interne energie van het totale systeem

dU = dU_1 + dU_2 = 2N k_B dT.

(4)

Voor het verplaatsen van de zuiger geldt voor beide kanten:

\delta W_i = P_i dV_i = \frac{N k_B T}{V_i} dV_i.

(5)

Als we deze twee formules samennemen dan levert dit:

-N k_B T \left( \frac{dV_1}{V_1} + \frac{dV_2}{V_2} \right) = 2 N k_B dT.

(6)

Vullen we de startformules voor $V_1$ en $V_2$ in, dan wordt dit:

\frac{dT}{dx} = \frac{x}{1-x^2}T.

(7)

{solution} ex-entropy-06 We kiezen als systeem het totale gekoppelde systeem (cv1 + cv2) en nemen aan dat er geen warmte-uitwisseling is met de omgeving. Dus $\delta Q_{\text{omg}}=0$ .
Met de eerste hoofdwet $dU = \delta Q - \delta W$ volgt dan direct:

dU = -\delta W.

(8)

De warmte die in de code tussen cv1 en cv2 “verplaatst” wordt is intern en telt voor het totale systeem niet als warmte naar de omgeving.

# Theorie: dT/dx = x/(1-x^2) T  =>  T(x) = T0 / sqrt(1 - x^2)

# x uit de volumes (totaal volume constant)
x = (volumes1 - volumes2) / (volumes1 + volumes2)

T_avg = 0.5 * (temperatures1 + temperatures2)
T0_sim = T_avg[0]
T_theory = T0_sim / np.sqrt(1 - x**2)

plt.figure(figsize=(7,4))
plt.plot(x, temperatures1, '.', ms=3, alpha=0.4, label='cv1 (sim)')
plt.plot(x, temperatures2, '.', ms=3, alpha=0.4, label='cv2 (sim)')
plt.plot(x, T_theory, '-', lw=2, label='theorie $T_0/\sqrt{1-x^2}$')
plt.xlabel('$x$')
plt.ylabel('$T$')
plt.legend()
plt.show()

print(f"T0 (sim) ≈ {T0_sim:.2f}")
print(f"T_theorie(x=-0.8) ≈ {T0_sim/np.sqrt(1-0.8**2):.2f}")

T0 (sim) ≈ 300.00
T_theorie(x=-0.8) ≈ 500.00

Om te verifiëren dat dit proces reversibel is moeten we het proces ook terug kunnen uitvoeren.

# Beweeg de zuiger terug en vergelijk heen- en terugweg

# begin vanuit de eindtoestand van de vorige simulatie (cv1 en cv2 bestaan al)
cv1.v_piston *= -1
cv2.v_piston *= -1

num_back = num_data
volumes1_back = np.zeros(num_back, dtype=float)
temperatures1_back = np.zeros(num_back, dtype=float)
volumes2_back = np.zeros(num_back, dtype=float)
temperatures2_back = np.zeros(num_back, dtype=float)

for i in range(num_back):
    cv1.take_time_step()
    cv2.take_time_step()

    average_temp = (cv1.temperature + cv2.temperature) / 2
    cv1.set_temp = average_temp
    cv2.set_temp = average_temp

    volumes1_back[i] = cv1.area
    temperatures1_back[i] = cv1.temperature
    volumes2_back[i] = cv2.area
    temperatures2_back[i] = cv2.temperature

x_fwd = (volumes1 - volumes2) / (volumes1 + volumes2)
x_back = (volumes1_back - volumes2_back) / (volumes1_back + volumes2_back)

T_fwd = 0.5 * (temperatures1 + temperatures2)
T_back = 0.5 * (temperatures1_back + temperatures2_back)

plt.figure(figsize=(7,4))
plt.plot(x_fwd, T_fwd, '-', label='heenweg')
plt.plot(x_back, T_back, '--', label='terugweg')
plt.xlabel('$x$')
plt.ylabel('$T$')
plt.legend()
plt.show()

print(f"Start (heen): x≈{x_fwd[0]:.3f}, T≈{T_fwd[0]:.2f}")
print(f"Einde (terug): x≈{x_back[-1]:.3f}, T≈{T_back[-1]:.2f}")

Start (heen): x≈-0.000, T≈300.00
Einde (terug): x≈0.000, T≈732.59

Solution to Exercise 11

De terugweg geeft meestal een kleine hysterese: de curve $T(x)$ op de terugweg ligt net iets anders dan op de heenweg en je eindigt niet exact op dezelfde $T$ als waarmee je begon.

Oorzaken:

het proces is niet perfect quasi-statisch (zuiger beweegt met eindige snelheid);
tijdens de beweging ontstaan kleine temperatuurverschillen en ruis, waardoor er entropie wordt geproduceerd;
numerieke effecten door discrete DT en overlap-detectie (soms “net te laat” botsen).

Reversibiliteit verbeteren:

kleinere |V_PISTON_0| (langzamer verplaatsen),
kleinere DT,
meer deeltjes (groter $N$ ) en/of groter volume (minder ruis),
eventueel extra “wachttijd” inbouwen zodat de volumes tussendoor kunnen uitmiddelen.

Helemaal perfect krijg je het praktisch niet, omdat er altijd fluctuaties zijn bij eindig $N$ en omdat de simulatie discreet is.

Adiabatische verplaatsing zuiger¶

We kunnen in een vergelijkbare eindsituatie komen door de zuiger eerst adiabatisch te verplaatsen (dat wil zeggen: zonder onderling thermisch contact tussen de volumes) en pas daarna het thermische contact toe te staan. Om dat in onze simulatie te modelleren concentreren we ons eerst op de adiabatische verplaatsing van de zuigerwand en nemen we het thermische contact nog niet mee.

ex-entropy-11 Voor een 2D ideaal (monoatomisch) gas geldt $\gamma = 2$ (want $C_V = k_B$ en $C_P = 2k_B$ per deeltje).
Bij een adiabatisch proces geldt:

T V^{\gamma-1} = \text{const} \quad \Rightarrow \quad T V = \text{const}.

(9)

We verplaatsen naar $x=-0.8$ , dus

V_1 = 0.2 V_0, \qquad V_2 = 1.8 V_0.

(10)

Daarmee:

T_1 = T_0 \frac{V_0}{0.2V_0} = 5T_0, \qquad T_2 = T_0 \frac{V_0}{1.8V_0} = \frac{5}{9}T_0.

(11)


```{exercise}
:label: ex-entropy-12
Maak een code voor een adiabatische verplaatsing voor de zuigerwand van de positie $x=0$ naar de positie $x=-0.8$.
Je kan hiervoor het best starten vanuit de code voor het isotherme proces en deze aanpassen.

# Adiabatische verplaatsing: geen thermisch contact tijdens het bewegen
# => zet set_temp = 0 zodat de thermostaat niets doet (temp_factor = 1)

num_data = round(abs(BOX_SIZE_0 * 0.4 / (DT * V_PISTON_0)))
if FAST_MODE:
    num_data = min(num_data, FAST_MAX_STEPS)

volumes1_ad = np.zeros(num_data, dtype=float)
temperatures1_ad = np.zeros(num_data, dtype=float)
volumes2_ad = np.zeros(num_data, dtype=float)
temperatures2_ad = np.zeros(num_data, dtype=float)

cv1 = ControlVolume(BOX_SIZE_0, BOX_SIZE_0, (30 if FAST_MODE else 50), +V_PISTON_0, 0.0)   # cv1 wordt kleiner
cv2 = ControlVolume(BOX_SIZE_0, BOX_SIZE_0, (30 if FAST_MODE else 50), -V_PISTON_0, 0.0)   # cv2 wordt groter

T0_ad = 0.5 * (cv1.temperature + cv2.temperature)

for i in range(num_data):
    cv1.take_time_step()
    cv2.take_time_step()

    volumes1_ad[i] = cv1.area
    temperatures1_ad[i] = cv1.temperature
    volumes2_ad[i] = cv2.area
    temperatures2_ad[i] = cv2.temperature

x_ad = (volumes1_ad - volumes2_ad) / (volumes1_ad + volumes2_ad)

plt.figure(figsize=(7,4))
plt.plot(x_ad, temperatures1_ad, '-', label='cv1 (adiabatisch)')
plt.plot(x_ad, temperatures2_ad, '-', label='cv2 (adiabatisch)')
plt.xlabel('$x$')
plt.ylabel('$T$')
plt.legend()
plt.show()

# Theoretische verwachting voor 2D ideaal gas: γ = 2 => T V^(γ-1) = T V = const
V0 = BOX_SIZE_0 * BOX_SIZE_0
T1_pred = T0_ad * (V0 / (0.2*V0))
T2_pred = T0_ad * (V0 / (1.8*V0))

print(f"T0 ≈ {T0_ad:.2f}")
print(f"Sim eind: T1≈{temperatures1_ad[-1]:.2f}, T2≈{temperatures2_ad[-1]:.2f}")
print(f"Theorie:  T1≈{T1_pred:.2f} (=5 T0), T2≈{T2_pred:.2f} (=5/9 T0)")

T0 ≈ 300.00
Sim eind: T1≈2465.38, T2≈208.15
Theorie:  T1≈1500.00 (=5 T0), T2≈166.67 (=5/9 T0)

Na de adiabatische verplaatsing moeten de twee volumes alsnog in thermisch evenwicht gebracht worden. Dat doen we met de code hieronder.

times = []
temp1 = []
temp2 = []

average_temp = (cv1.temperature + cv2.temperature) / 2
cv1.set_temp = average_temp
cv2.set_temp = average_temp
cv1.v_piston = 0.0
cv2.v_piston = 0.0
time = 0.0

max_iter = 5000 if FAST_MODE else 200000
iter_count = 0

while (cv2.temperature < 0.99 * cv2.set_temp) and (iter_count < max_iter):
    time += DT
    iter_count += 1
    cv1.take_time_step()
    cv2.take_time_step()
    temp1.append(cv1.temperature)
    temp2.append(cv2.temperature)
    times.append(time)

if iter_count >= max_iter:
    print('Let op: maximaal aantal stappen bereikt; evenwicht mogelijk nog niet volledig (FAST_MODE).')

plt.figure()
plt.xlabel('$t$')
plt.ylabel('$T$')

plt.plot(times, temp1, '-r', label='cv1')
plt.plot(times, temp2, 'b-', label='cv2')
plt.legend()

plt.show()

Vergelijk van de twee experimenten¶

We hebben de zuiger in twee experimenten verschoven:

We begonnen in toestand 0, waarbij de zuiger in het midden stond en naar toestand 1 bewoog, met thermisch contact tussen de twee volumes.
We begonnen in toestand 0 en verplaatsten de zuiger zonder thermisch contact (adiabatisch) naar toestand 2. Daarna brachten we de twee volumes in thermisch contact en zonder volumeverandering naar evenwicht kwamen in toestand 3.

Deze processen zijn schematisch weergegeven in Figure 2.

Schematische weergave van de twee processen in een (p,V)-diagram met de verschillende toestanden genummerd weergegeven. — Figure 2:Schematische weergave van de twee processen in een $(p,V)$ -diagram met de verschillende toestanden genummerd weergegeven.

Uit de thermodynamische formules en uit de simulatie kun je concluderen dat de eindtemperatuur na deze twee experimenten verschillend is. Omdat we het systeem nooit in thermisch contact met de omgeving hebben gebracht, moet de energie hiervoor uit de arbeid komen die de zuigerwand heeft verricht. Dit laat zien dat de arbeid die nodig is van de ene in de andere toestand te komen, afhankelijk is van het gekozen proces.

ex-entropy-16 We kiezen als systeem steeds het totale systeem (cv1 + cv2), zonder warmte-uitwisseling met de omgeving. Dan geldt:

W_{\text{op}} = \Delta U.

(15)

Proces 0 → 1 (met thermisch contact tijdens bewegen)
De theorie geeft bij $x=-0.8$ :

T_1 = \frac{T_0}{\sqrt{1-x^2}} = \frac{T_0}{\sqrt{1-0.8^2}} = \frac{5}{3}T_0.

(16)

Totale interne energie: $U = 2N k_B T$ (met $N$ het aantal deeltjes per volume).
Dus:

\Delta U_{01} = 2N k_B (T_1 - T_0) = 2N k_B\left(\frac{5}{3}-1\right)T_0 = \frac{4}{3}N k_B T_0.

(17)

Proces 0 → 2 → 3 (adiabatisch bewegen, daarna alleen warmte-uitwisseling bij stilstaande zuiger)
Na de adiabatische stap:

T_1 = 5T_0, \qquad T_2 = \frac{5}{9}T_0

(18)

en dus

U_2 = N k_B(5T_0) + N k_B\left(\frac{5}{9}T_0\right) = \frac{50}{9}N k_B T_0.

(19)

Start: $U_0 = 2N k_B T_0 = \frac{18}{9}N k_B T_0$ .
Dus:

\Delta U_{02} = \left(\frac{50}{9}-\frac{18}{9}\right)N k_B T_0 = \frac{32}{9}N k_B T_0.

(20)

De stap 2→3 heeft geen arbeid (zuiger staat stil), dus $W_{023}=\Delta U_{02}$ .

Verschil in verrichte arbeid

\Delta W = W_{023} - W_{01} = \left(\frac{32}{9} - \frac{4}{3}\right)N k_B T_0 = \frac{20}{9}N k_B T_0.

(21)

Met in de simulatie $N=50$ per volume en $T_0\approx 300\,\text{K}$ :

\Delta W \approx \frac{20}{9}\cdot 50 \cdot k_B \cdot 300 \approx 2.3\cdot 10^{-19}\,\text{J}.

(22)


Een reversibele route is dan het 'goedkoopste' proces om in een toestand te komen. Voor het experiment zou het perfect uitgevoerde isotherme proces reversibel zijn en daarom de laagste temperatuur opleveren bij $x=-0.8$.

```{exercise} entropy production
:label: ex-entropy-17
Bereken uit de theoretische, thermodynamische formules hoeveel entropie er vrijkomt bij het tweede proces.

ex-entropy-17 De entropieproductie is gelijk aan de totale entropieverandering van het geïsoleerde systeem (er is geen warmte-uitwisseling met de omgeving).
Voor een 2D ideaal gas geldt per volume:

\Delta S = N k_B \ln\frac{T_f}{T_i} + N k_B \ln\frac{V_f}{V_i}.

(23)

Voor toestand 0: $T_i=T_0$ , $V_i=V_0$ .
Voor toestand 3: $T_f=T_{eq}=\frac{25}{9}T_0$ , $V_1=0.2V_0$ , $V_2=1.8V_0$ .

Sommeer voor beide volumes:

\Delta S_{sys} = 2N k_B \ln\frac{25}{9} + N k_B \ln(0.2) + N k_B \ln(1.8) = 2N k_B \ln\frac{25}{9} + N k_B \ln(0.36).

(24)

Omdat $0.36 = \frac{9}{25}$ , valt dit mooi samen:

\Delta S_{sys} = 2N k_B \ln\frac{25}{9} + N k_B \ln\frac{9}{25} = N k_B \ln\frac{25}{9}.

(25)

Dus de geproduceerde entropie is:

S_{gen} = N k_B \ln\frac{25}{9}.

(26)

Willen we maximale arbeid halen in de stap van toestand 2 naar toestand 3, dan kunnen we dat doen met twee reversibele warmtepompen die de twee beheersvolumes koppelen aan een denkbeeldig bad. Entropie is een toestandsgrootheid, zodat de entropieproductie in dit proces net zo groot moet zijn als in onze zuiger van daarnet. Voor een reversibel proces zoals we met deze warmtepompen uitvoeren is er geen entropieproductie. Dan moet er dus gelden dat:

\Delta S_{\text{sys}} + \Delta S_{\text{omg}} = 0

(27)

Als we de omgevingstemperatuur op een waarde stellen van $T_{\text{omg}}$ dan moet er dus gelden:

\Delta S_{\text{omg}} = -\frac{Q_{\text{rev}}}{T_{\text{omg}}} \Rightarrow Q_{\text{rev}} = T_{\text{omg}} \Delta S_{\text{sys}}

(28)

Dit betekent dat het systeem tijdens dit reversibele proces deze hoeveelheid $Q_{\text{rev}}$ moet accepteren van de omgeving om de eindtoestand te bereiken.

Omdat we daarnaast van de eerste hoofdwet weten dat:

\Delta U_{\text{sys}} = Q_{\text{rev}} - W_{\text{rev}}

(29)

en er moet gelden dat de interne energie van het systeem niet verandert tijdens dit temperatuurverloop kunnen we concluderen dat:

W_{\text{max}} = W_{\text{rev}} = T_{\text{omg}} \Delta S_{\text{sys}} = N k_B T_{\text{omg}} \ln{\frac{25}{9}}

(30)

In het proces dat bij onze simulatie heeft plaatsgevonden is er geen enkele arbeid door het systeem uitgevoerd. We zien dus dat de entropie die we in dit geval hebben berekend een maat is voor de arbeid die we hadden kunnen winnen door op dezelfde eindtoestand te komen via een reversibel proces.

We kunnen dit ook andersom zeggen: de hoeveelheid entropie die tijdens het proces gecreëerd wordt is een maat voor de hoeveelheid beschikbare arbeid die we tijdens het proces hebben verloren door het niet perfect reversibel uit te voeren.