Doelen¶
Nu we de microscopische grootheden van de moleculen hebben verbonden aan de macroscopische grootheden van het gas kunnen we de thermodynamica van het gas echt bestuderen met onze simulatie. In dit werkblad kijken we hoe de temperatuur en de druk veranderen onder invloed van een zuiger die het volume verandert.
Eerst herhalen we de delen van de code die we nodig hebben:
klasse voor het deeltje met bijbehorende functies
variabelen en randcondities van de controle volume
functies voor (een lijst) deeltjes
Daarna voegen we code toe voor de dynamiek van de zuiger:
zuiger implementeren in volume en dynamische formules
En vervolgens:
bestuderen van temperatuur en druk als functie van volume
onderzoeken of we terug kunnen keren naar startcondities
In onderstaande animatie laten we het proces zien dat je gaat programmeren.
Laden van eerdere code¶
We beginnen weer met de noodzakelijke pakketten en de constanten. Daar voegen we nu een constante aan toe: de startsnelheid van de zuiger.
# ruimte voor uitwerking
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.optimize import curve_fit
BOX_SIZE_0 = 10.0 # Hoogte/breedte startvolume (in nm)
N = 40 # Aantal deeltjes
V_0 = 410.0 # Startsnelheid van deeltjes (in nm per ps)
RADIUS = 0.3 # Straal van moleculen (in nm)
DT = 5e-5 # Tijdstap (in ps)
V_PISTON_0 = -0.1 * V_0 # Startsnelheid van zuiger
# (negatief betekent zowel links als rechts naar binnen gericht)
#your code/answer
Zoals altijd laden we de klasse voor de gasmoleculen en de functies voor hun onderlinge interactie:
class ParticleClass:
def __init__(self, m, v, r, R):
""" maakt een deeltje (constructor) """
self.m = m
self.v = np.array(v, dtype=float)
self.r = np.array(r, dtype=float)
self.R = R
def update_position(self):
""" verandert positie voor één tijdstap """
self.r += self.v * DT
@property
def momentum(self):
return self.m * self.v
@property
def kin_energy(self):
return 1/2 * self.m * np.dot(self.v, self.v)
def collide_detection(p1: ParticleClass, p2: ParticleClass) -> bool:
""" Geeft TRUE als de deeltjes overlappen """
dx = p1.r[0]-p2.r[0]
dy = p1.r[1]-p2.r[1]
rr = p1.R + p2.R
return dx**2+dy**2 < rr**2
def particle_collision(p1: ParticleClass, p2: ParticleClass):
""" past snelheden aan uitgaande van overlap """
m1, m2 = p1.m, p2.m
delta_r = p1.r - p2.r
delta_v = p1.v - p2.v
dot_product = np.dot(delta_r, delta_v)
# 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_rHet volume en de randvoorwaarden zullen we moeten aanpassen aan onze simulatie met bewegende zuiger: Het volume zal nu niet meer altijd een vierkant zijn.
De simulatie bestaat uit een volume met links en rechts een bewegende wand: de zuiger.
Laten we aannemen dat de zuiger altijd in de horizontale richting verplaatst en het volume symmetrisch houdt ten opzichte van de oorsprong, d.w.z. er is een zuiger aan de linker wand die een tegengestelde verplaatsing heeft aan die in de rechter wand.
We maken eerst een aantal variabelen aan die bij het volume horen:
box_height = BOX_SIZE_0 # 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 gasDe 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 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
elif particle.r[0] - particle.R < -7+box_length/2:
impulse_outward += abs(particle.momentum[0]) * 2
particle.v[0] *= -1En 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(-7+box_length/2, 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:
energy_array = []
for p in particles:
energy = 4.81*10**(-26)*p.v**2
energy_array = np.append(energy_array,energy)
energy_average = np.mean(energy_array)
temp = 0.0
k = 1.38*10**-23
f = 2
temp = (2*energy_average)/(k*f)
return temp
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 . 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...
# we kiezen het aantal datapunten zodat het volume tot 1/3 van het begin volume reduceert
num_steps = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
box_length = BOX_SIZE_0 # zetten zuiger terug
v_piston = V_PISTON_0
create_particles(particles) # resetten deeltjes
for i in range(num_steps):
take_time_step(particles)
volumes[i] = box_length * box_height
temperatures[i] = temperature(particles)
temperatures = np.asarray(temperatures)
plt.figure()
plt.xlabel('Volume')
plt.ylabel('Temperature')
plt.plot(volumes, temperatures, '-r')
plt.show()
Dit kan niet kloppen. Hier zien we dat de temperatuur vrijwel constant is (let op de vermenigvuldigingsfactor vermeld aan de bovenkant van de verticale as). Maar de zuiger voert arbeid uit, op baiss van de wet van behoud van energie betekent zou de temperatuur moet veranderen!
Om het model kloppend te maken moeten we kijken naar de botsing van de deeltjes met de wand. In de vorige werkbladen stonden de wanden stil en veranderde de snelheid van de deeltjes alleen van teken in de component loodrecht op de wand. Nu dat de wanden een zuiger zijn en snelheid hebben, moeten we daarvoor corrigeren.
De snelheid van de deeltjes klapt nog steeds om van teken in het referentiestelsel van de wand, maar omdat de wand beweegt ten opzichte van het volume met snelheid , wordt de juiste functie:
def left_right_collision(particle: ParticleClass):
""" verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
global box_length, v_piston, impulse_outward
if particle.r[0] + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = 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)
elif particle.r[0] - particle.R < -7+box_length/2:
piston_velocity = 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) 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)
#your code/answer
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.
pressure = 0
pressure = 0.95 * pressure + 0.05 * impulse_outward / ((2 * box_length + 2 * box_height) * DT)
# 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)
druk = 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
druk[i] = pressure
plt.figure()
plt.xlabel('Volume')
plt.ylabel('Pressure')
druk = np.asarray(druk)
#your code/answer
plt.plot(volumes, druk, '-r')
plt.show()
Als je simulatie klopt heeft deze grafiek de vorm van een machtsfunctie met een negatieve exponent.
# pv^2 = constantdef power_law(vol, a, n):
""" de fitfunctie voor het P,V-diagram """
return a * (vol)**n
popt, cov = curve_fit(power_law, volumes, druk, p0=(1, -1))
x1 = np.linspace(min(volumes), max(volumes), 1000)
y1 = power_law(x1, *popt)
print(y1)
plt.figure()
plt.plot(volumes, druk, '-r')
plt.plot(x1, y1, '--k')
plt.xlabel('Volume (m³)')
plt.ylabel('Pressure (Pa)')
plt.show()
print("de machtconstante is", popt[1])[411258.43972676 410053.07863561 408853.64415648 407660.09542567
406472.39194177 405290.4935617 404114.36049684 402943.95330924
401779.23290779 400620.16054455 399466.69781105 398318.80663463
397176.44927489 396039.58832014 394908.1866839 393782.20760143
392661.61462636 391546.37162728 390436.44278445 389331.79258649
388232.38582716 387138.18760213 386049.16330585 384965.27862841
383886.49955244 382812.79235009 381744.12358001 380680.46008437
379621.76898594 378568.01768516 377519.17385731 376475.20544966
375436.08067869 374401.7680273 373372.23624211 372347.45433075
371327.3915592 370312.01744917 369301.30177548 368295.21456351
367293.72608665 366296.80686381 365304.42765693 364316.55946854
363333.17353931 362354.24134573 361379.73459768 360409.62523612
359443.88543081 358482.487578 357525.40429817 356572.60843384
355624.07304736 354679.7714187 353739.67704336 352803.7636302
351872.00509937 350944.37558022 350020.84940925 349101.4011281
348186.00548153 347274.63741544 346367.27207493 345463.88480235
344564.45113539 343668.94680519 342777.34773448 341889.63003571
341005.77000926 340125.7441416 339249.52910351 338377.10174834
337508.43911023 336643.51840242 335782.31701551 334924.8125158
334070.98264359 333220.80531159 332374.25860322 331531.32077104
330691.97023513 329856.18558157 329023.94556078 328195.22908609
327370.01523213 326548.28323335 325730.01248256 324915.18252941
324103.77307895 323295.76399018 322491.13527465 321689.867095
320891.93976362 320097.33374123 319306.02963551 318518.0081998
317733.25033171 316951.7370718 316173.44960233 315398.36924589
314626.47746419 313857.75585674 313092.18615963 312329.75024428
311570.43011622 310814.20791389 310061.06590738 309310.98649734
308563.95221374 307819.94571471 307078.94978542 306340.94733693
305605.92140505 304873.85514926 304144.73185158 303418.5349155
302695.24786487 301974.85434288 301257.33811096 300542.68304778
299830.87314817 299121.89252212 298415.72539377 297712.35610042
297011.76909149 296313.94892759 295618.88027954 294926.54792739
294236.93675947 293550.03177146 292865.81806546 292184.28084907
291505.40543446 290829.17723747 290155.58177674 289484.6046728
288816.23164719 288150.44852163 287487.2412171 286826.59575307
286168.49824658 285512.93491145 284859.89205746 284209.35608951
283561.31350683 282915.75090218 282272.65496105 281632.01246087
280993.81027025 280358.0353482 279724.6747434 279093.71559339
278465.14512388 277838.95064796 277215.11956543 276593.63936203
275974.49760873 275357.68196104 274743.18015828 274130.98002293
273521.06945988 272913.43645578 272308.0690784 271704.95547588
271104.08387615 270505.4425862 269909.01999151 269314.80455533
268722.7848181 268132.94939678 267545.28698425 266959.78634869
266376.43633296 265795.22585398 265216.14390216 264639.17954078
264064.32190541 263491.56020333 262920.88371294 262352.28178318
261785.74383302 261221.2593508 260658.81789377 260098.40908748
259540.02262525 258983.64826763 258429.27584185 257876.89524134
257326.49642513 256778.06941738 256231.60430686 255687.0912464
255144.52045244 254603.88220449 254065.16684464 253528.36477708
252993.46646759 252460.46244309 251929.34329113 251400.09965941
250872.72225535 250347.20184561 249823.52925558 249301.69536901
248781.69112746 248263.50752995 247747.13563244 247232.56654741
246719.79144344 246208.80154477 245699.58813086 245192.14253598
244686.45614878 244182.52041186 243680.3268214 243179.86692671
242681.13232982 242184.11468511 241688.80569891 241195.19712906
240703.28078458 240213.04852522 239724.49226113 239237.60395245
238752.37560894 238268.79928959 237786.86710229 237306.57120342
236827.90379748 236350.8571368 235875.42352109 235401.59529714
234929.36485847 234458.72464494 233989.66714246 233522.1848826
233056.27044229 232591.91644343 232129.11555263 231667.86048081
231208.14398292 230749.95885759 230293.29794682 229838.15413565
229384.52035185 228932.38956562 228481.75478925 228032.60907683
227584.94552394 227138.75726737 226694.03748476 226250.77939438
225808.97625477 225368.62136447 224929.70806177 224492.22972433
224056.179769 223621.55165146 223188.33886597 222756.53494509
222326.13345942 221897.12801728 221469.5122645 221043.2798841
220618.42459605 220194.940157 219772.82036003 219352.05903436
218932.65004512 218514.5872931 218097.86471445 217682.47628051
217268.41599746 216855.67790616 216444.25608187 216034.14463399
215625.33770586 215217.82947449 214811.61415032 214406.68597703
214003.03923124 213600.66822234 213199.56729221 212799.73081505
212401.1531971 212003.82887645 211607.75232279 211212.91803724
210819.32055208 210426.95443055 210035.81426666 209645.89468493
209257.19034023 208869.69591755 208483.40613177 208098.31572749
207714.4194788 207331.7121891 206950.18869087 206569.84384552
206190.67254312 205812.66970227 205435.83026987 205060.14922096
204685.62155847 204312.24231309 203940.00654306 203568.90933397
203198.9457986 202830.1110767 202462.40033486 202095.80876629
201730.33159063 201365.96405381 201002.70142788 200640.53901076
200279.47212615 199919.49612333 199560.60637696 199202.79828696
198846.06727829 198490.40880084 198135.81832921 197782.29136257
197429.8234245 197078.41006285 196728.04684951 196378.72938032
196030.45327488 195683.21417642 195337.00775158 194991.82969035
194647.67570584 194304.54153415 193962.42293424 193621.31568776
193281.2155989 192942.11849425 192604.02022267 192266.9166551
191930.80368447 191595.6772255 191261.53321463 190928.3676098
190596.17639038 190264.95555697 189934.70113132 189605.40915616
189277.07569506 188949.69683232 188623.26867282 188297.78734188
187973.24898517 187649.64976852 187326.98587783 187005.25351894
186684.44891749 186364.56831879 186045.60798773 185727.5642086
185410.43328501 185094.21153976 184778.89531468 184464.48097057
184150.96488705 183838.34346242 183526.61311359 183215.77027592
182905.81140312 182596.73296714 182288.53145804 181981.2033839
181674.74527069 181369.15366215 181064.42511971 180760.55622233
180457.54356644 180155.38376582 179854.07345146 179553.6092715
179253.98789106 178955.20599223 178657.26027385 178360.14745152
178063.8642574 177768.40744016 177473.77376488 177179.96001292
176886.96298185 176594.77948532 176303.40635299 176012.84043042
175723.07857895 175434.11767566 175145.9546132 174858.58629977
174572.00965896 174286.22162969 174001.21916613 173716.99923758
173433.55882836 173150.89493779 172869.00458002 172587.88478399
172307.53259332 172027.94506624 171749.11927546 171471.05230814
171193.74126576 170917.18326404 170641.37543289 170366.31491628
170091.99887216 169818.4244724 169545.58890272 169273.48936254
169002.12306496 168731.48723668 168461.57911787 168192.39596214
167923.93503641 167656.19362089 167389.16900896 167122.85850709
166857.2594348 166592.36912454 166328.18492164 166064.70418422
165801.92428312 165539.84260183 165278.45653642 165017.76349544
164757.76089988 164498.44618306 164239.8167906 163981.87018032
163724.60382218 163468.01519818 163212.10180236 162956.86114065
162702.29073085 162448.38810253 162195.150797 161942.57636721
161690.66237769 161439.40640449 161188.80603512 160938.85886844
160689.56251466 160440.91459523 160192.91274279 159945.55460109
159698.83782497 159452.76008023 159207.31904363 158962.51240278
158718.33785611 158474.79311279 158231.87589269 157989.58392628
157747.9149546 157506.86672922 157266.43701211 157026.62357566
156787.42420257 156548.8366858 156310.85882853 156073.4884441
155836.72335592 155600.56139745 155365.00041213 155130.03825333
154895.67278428 154661.90187801 154428.72341733 154196.13529474
153964.13541241 153732.72168206 153501.892025 153271.64437199
153041.97666325 152812.88684837 152584.37288627 152356.43274514
152129.06440242 151902.26584471 151676.03506773 151450.37007628
151225.26888418 151000.72951423 150776.74999815 150553.32837653
150330.46269879 150108.15102313 149886.39141648 149665.18195444
149444.52072125 149224.40580975 149004.83532129 148785.80736573
148567.32006139 148349.37153495 148131.95992147 147915.08336433
147698.74001515 147482.92803377 147267.64558822 147052.89085464
146838.66201726 146624.95726836 146411.77480821 146199.11284504
145986.96959498 145775.34328204 145564.23213805 145353.63440262
145143.54832312 144933.97215459 144724.90415977 144516.34260897
144308.28578012 144100.73195867 143893.67943755 143687.12651717
143481.07150535 143275.51271728 143070.44847549 142865.87710982
142661.79695734 142458.20636238 142255.10367643 142052.48725811
141850.35547319 141648.70669447 141447.53930179 141246.851682
141046.64222889 140846.90934318 140647.65143248 140448.86691124
140250.55420071 140052.71172894 139855.33793071 139658.4312475
139461.99012748 139266.01302543 139070.49840274 138875.44472738
138680.85047384 138486.71412311 138293.03416263 138099.8090863
137907.03739439 137714.71759355 137522.84819676 137331.42772329
137140.45469868 136949.9276547 136759.84512933 136570.20566671
136381.00781712 136192.25013694 136003.93118863 135816.04954069
135628.60376764 135441.59244997 135255.01417413 135068.86753247
134883.15112325 134697.86355058 134513.00342441 134328.56936047
134144.55998028 133960.97391108 133777.80978584 133595.0662432
133412.74192746 133230.83548855 133049.34558197 132868.27086883
132687.61001573 132507.36169483 132327.52458373 132148.09736552
131969.07872871 131790.46736719 131612.26198026 131434.46127253
131257.06395395 131080.06873977 130903.47435049 130727.27951185
130551.48295483 130376.08341556 130201.07963537 130026.47036069
129852.25434309 129678.43033921 129504.99711076 129331.95342446
129159.29805208 128987.02977033 128815.14736093 128643.64961048
128472.53531054 128301.80325752 128131.45225272 127961.48110226
127791.88861708 127622.67361292 127453.83491027 127285.37133438
127117.28171522 126949.56488744 126782.21969037 126615.24496801
126448.63956898 126282.40234648 126116.53215833 125951.02786689
125785.88833907 125621.11244628 125456.69906445 125292.64707396
125128.95535964 124965.62281077 124802.64832102 124640.03078845
124477.7691155 124315.86220892 124154.30897983 123993.1083436
123832.25921992 123671.76053273 123511.61121021 123351.81018475
123192.35639294 123033.24877557 122874.48627757 122716.06784802
122557.9924401 122400.25901111 122242.86652242 122085.81393947
121929.10023172 121772.72437268 121616.68533984 121460.98211467
121305.61368264 121150.57903313 120995.87715945 120841.50705883
120687.46773239 120533.7581851 120380.37742583 120227.32446722
120074.59832578 119922.19802179 119770.12257933 119618.37102621
119466.94239402 119315.83571807 119165.05003735 119014.58439459
118864.43783615 118714.60941208 118565.09817605 118415.90318536
118267.02350093 118118.45818724 117970.20631237 117822.26694794
117674.63916913 117527.32205461 117380.31468658 117233.61615073
117087.22553622 116941.14193566 116795.36444511 116649.89216406
116504.72419541 116359.85964543 116215.29762381 116071.03724356
115927.07762107 115783.41787605 115640.05713153 115496.99451382
115354.22915255 115211.7601806 115069.58673411 114927.70795246
114786.12297826 114644.83095732 114503.83103867 114363.12237451
114222.70412019 114082.57543423 113942.7354783 113803.18341717
113663.91841873 113524.93965397 113386.24629696 113247.83752483
113109.71251779 112971.87045904 112834.31053486 112697.03193451
112560.03385026 112423.31547736 112286.87601404 112150.71466148
112014.8306238 111879.22310808 111743.89132428 111608.83448529
111474.0518069 111339.54250774 111205.30580935 111071.34093611
110937.64711522 110804.22357673 110671.06955351 110538.18428121
110405.5669983 110273.21694599 110141.13336829 110009.31551194
109877.76262645 109746.47396401 109615.44877958 109484.68633079
109354.18587798 109223.94668415 109093.96801499 108964.24913884
108834.78932668 108705.58785213 108576.64399142 108447.95702341
108319.52622955 108191.35089386 108063.43030296 107935.76374603
107808.35051479 107681.18990352 107554.281209 107427.62373058
107301.21677007 107175.05963181 107049.1516226 106923.49205174
106798.08023098 106672.91547454 106547.99709906 106423.32442364
106298.89676977 106174.71346139 106050.77382482 105927.07718876
105803.62288431 105680.41024494 105557.43860647 105434.70730707
105312.21568726 105189.96308988 105067.94886008 104946.17234535
104824.63289546 104703.32986246 104582.2626007 104461.43046678
104340.83281958 104220.46902021 104100.33843205 103980.44042067
103860.77435391 103741.33960177 103622.1355365 103503.16153251
103384.41696641 103265.90121698 103147.61366517 103029.55369408
102911.72068895 102794.11403718 102676.73312828 102559.57735388
102442.64610773 102325.93878567 102209.45478565 102093.19350769
101977.15435388 101861.33672839 101745.74003745 101630.36368931
101515.2070943 101400.26966476 101285.55081505 101171.04996155
101056.76652265 100942.69991874 100828.84957217 100715.21490732
100601.79535049 100488.59033 100375.59927607 100262.82162091
100150.25679865 100037.90424536 99925.76339902 99813.83369953
99702.11458873 99590.6055103 99479.30590986 99368.2152349
99257.33293478 99146.65846072 99036.19126583 98925.93080505
98815.87653516 98706.0279148 98596.38440442 98486.9454663
98377.71056454 98268.67916504 98159.85073549 98051.2247454
97942.80066604 97834.57797047 97726.55613351 97618.73463176
97511.11294356 97403.690549 97296.46692993 97189.4415699
97082.61395423 96975.98356992 96869.5499057 96763.31245202
96657.27070101 96551.4241465 96445.77228401 96340.31461071
96235.05062548 96129.97982886 96025.10172301 95920.41581179
95815.92160067 95711.61859678 95607.50630886 95503.58424731
95399.8519241 95296.30885285 95192.95454877 95089.78852867
94986.81031096 94884.01941562 94781.41536422 94678.9976799
94576.76588737 94474.7195129 94372.85808432 94271.18113099
94169.68818382 94068.37877526 93967.2524393 93866.30871143
93765.54712866 93664.96722954 93564.56855409 93464.35064385
93364.31304183 93264.45529257 93164.77694205 93065.27753773
92965.95662858 92866.81376497 92767.84849879 92669.06038333
92570.44897336 92472.01382509 92373.75449614 92275.67054559
92177.76153391 92080.02702303 91982.46657624 91885.07975829
91787.86613529 91690.82527478 91593.95674565 91497.26011822
91400.73496415 91304.38085651 91208.1973697 91112.18407951
91016.34056308 90920.6663989 90825.16116682 90729.82444802
90634.65582501 90539.65488165 90444.82120311 90350.15437589
90255.65398781 90161.31962799 90067.15088685 89973.14735614
89879.30862888 89785.63429939 89692.12396326 89598.77721739
89505.59365994 89412.57289033 89319.71450926 89227.01811869
89134.48332183 89042.10972314 88949.89692834 88857.84454438
88765.95217945 88674.21944297 88582.64594559 88491.23129918
88399.97511683 88308.87701284 88217.93660273 88127.1535032
88036.52733218 87946.05770876 87855.74425326 87765.58658715
87675.58433311 87585.73711497 87496.04455775 87406.50628763
87317.12193196 87227.89111924 87138.81347913 87049.88864244
86961.11624113 86872.49590828 86784.02727814 86695.70998606
86607.54366854 86519.52796321 86431.66250878 86343.94694513
86256.38091321 86168.96405509 86081.69601396 85994.57643407
85907.60496081 85820.78124063 85734.10492107 85647.57565076
85561.19307942 85474.95685781 85388.86663779 85302.92207227
85217.12281525 85131.46852174 85045.95884785 84960.59345071
84875.37198852 84790.29412051 84705.35950694 84620.56780912
84535.91868938 84451.4118111 84367.04683864 84282.82343743
84198.74127388 84114.80001542 84030.99933049 83947.33888854
83863.81836002 83780.43741637 83697.19573002 83614.09297441
83531.12882394 83448.30295401 83365.61504099 83283.06476223
83200.65179605 83118.37582174 83036.23651954 82954.23357069
82872.36665733 82790.6354626 82709.03967057 82627.57896627
82546.25303565 82465.06156563 82384.00424403 82303.08075965
82222.29080217 82141.63406224 82061.11023139 81980.71900211]
de machtconstante is -1.4685520594522676
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
# Eerste hoofdwet¶
Een goede controle voor ons simulatie is te onderzoeken of deze voldoet aan de eerste hoofdwet van de thermodynamica.
Omdat er (nog) geen warmte wordt uitgewisseld, moeten we alleen de hoeveelheid arbeid bepalen. De arbeid wordt geleverd door de zuiger(s) op de deeltjes. Dus in de functie waar we de botsingen van de deeltjes met de wand behandelen, kunnen we ook de verrichte arbeid bepalen. De arbeid wordt gegeven\ door:
Hier staan de 1 voor de begintoestand en de 2 voor de eindtoestand van het proces. In de differentiaalvorm wordt dit geschreven als:
Zoals dat in het boek gebeurt, kiezen we voor de notatie met in plaats van , om aan te geven dat deze integraal afhankelijk is van het gekozen proces en niet alleen van het begin- en eindpunt. Voor ons experiment, waar het volume niet geleidelijk maar in stapjes verandert, kan je de hoeveelheid arbeid per tijdstip dus herschrijven tot:
waarin de hoogte is van de zuiger.
Dit betekent dat we de arbeid per tijdstap in onze code kunnen bepalen op hetzelfde moment dat we de druk bepalen met behulp van de gezamenlijke stoot van de moleculen. Daarvoor moeten we de code voor de functie left_right_collision nogmaals aanpassen:
def left_right_collision(particle: ParticleClass):
""" verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
global box_length, v_piston, impulse_outward, work
if abs(particle.r[0]) + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = np.sign(particle.r[0]) * v_piston
relative_velocity = particle.v[0] - piston_velocity # stelsel zuiger
particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
impulse_outward += 2 * particle.m * abs(relative_velocity)
work += 2 * particle.m * relative_velocity * piston_velocityVoordat we een simulatie draaien is het goed even stil te staan en na te denken. Net als bij een experiment in het lab doorlopen we bewust de onderzoekscyclus. We moeten daarvoor eerst een hypothese opstellen en daarna pas de simulatie uitvoeren.
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
num_steps = 1000
tijd = np.zeros(num_steps, dtype=float)
work1 = np.zeros(num_steps, dtype=float)
create_particles(particles) # resetten deeltjes
for i in range(num_steps):
take_time_step(particles)
tijd[i] = DT + tijd[i-1]
work1[i] = work
plt.figure()
plt.xlabel('Time')
plt.ylabel('work')
work1 = np.asarray(work1)
#your code/answer
plt.plot(tijd, work1, '-r')
plt.show()
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.
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
num_steps = 1000
tijd = np.zeros(num_steps, dtype=float)
work1 = 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)
tijd[i] = DT + tijd[i-1]
work1[i] = work
plt.figure()
plt.xlabel('Time')
plt.ylabel('work')
plt.plot(tijd, work1, '-r')
plt.show()
work1 = np.asarray(work1)
v_piston = -V_PISTON_0 # zuiger richting wordt omgedraaid
num_steps = 1000
tijd2 = np.zeros(num_steps, dtype=float)
tijd2[0] = tijd[-1]
work1 = np.zeros(num_steps, dtype=float)
create_particles(particles) # resetten deeltjes
for i in range(num_steps):
take_time_step(particles)
tijd2[i] = DT + tijd2[i-1]
work1[i] = work
plt.figure()
plt.xlabel('Time')
plt.ylabel('work')
plt.plot(tijd2, work1, '-r')
plt.show()
work1 = np.asarray(work1)
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
num_steps = 200
tijd = np.zeros(num_steps, dtype=float)
work1 = 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)
tijd[i] = DT + tijd[i-1]
work1[i] = work
plt.figure()
plt.xlabel('Time')
plt.ylabel('work')
plt.plot(tijd, work1, '-r')
plt.show()
work1 = np.asarray(work1)
v_piston = -V_PISTON_0 # zuiger richting wordt omgedraaid
num_steps = 200
tijd2 = np.zeros(num_steps, dtype=float)
tijd2[0] = tijd[-1]
work1 = np.zeros(num_steps, dtype=float)
create_particles(particles) # resetten deeltjes
for i in range(num_steps):
take_time_step(particles)
tijd2[i] = DT + tijd2[i-1]
work1[i] = work
plt.figure()
plt.xlabel('Time')
plt.ylabel('work')
plt.plot(tijd2, work1, '-r')
plt.show()
work1 = np.asarray(work1)
#your code/answer
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.