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 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] *= -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(-box_length/2 + RADIUS, box_length/2 - RADIUS)
y = np.random.uniform(-box_height/2 + RADIUS, box_height/2 - RADIUS)
particles.append(ParticleClass(m=1.0, v=[vx, vy], r=[x, y], R=RADIUS))
def temperature(particles) -> float:
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 abs(particle.r[0]) + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = np.sign(particle.r[0]) * v_piston
relative_velocity = particle.v[0] - piston_velocity # stelsel zuiger
particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
impulse_outward += 2 * particle.m * abs(relative_velocity) # stoot gevoeld door zuigerNu kunnen we de simulatie opnieuw uitvoeren:
num_steps = round(2/3 * BOX_SIZE_0 / (2 * -V_PISTON_0 * DT))
particles = []
volumes = np.zeros(num_steps, dtype=float)
temperatures = np.zeros(num_steps, dtype=float)
box_length = BOX_SIZE_0 # 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])[1432622.79609873 1422568.56011635 1412604.73017532 1402730.31562656
1392944.3386161 1383245.83389494 1373633.84863198 1364107.44223014
1354665.6861456 1345307.66370995 1336032.46995544 1326839.21144302
1317727.00609333 1308694.98302042 1299742.28236827 1290868.05514997
1282071.46308963 1273351.67846674 1264707.88396329 1256139.27251318
1247645.04715429 1239224.42088282 1230876.61651006 1222600.8665215
1214396.41293824 1206262.50718055 1198198.40993378 1190203.39101631
1182276.72924971 1174417.71233093 1166625.63670658 1158899.80744921
1151239.53813556 1143644.15072677 1136112.97545051 1128645.35068491
1121240.62284444 1113898.14626752 1106617.28310592 1099397.40321595
1092237.8840513 1085138.11055761 1078097.47506867 1071115.37720426
1064191.2237696 1057324.42865634 1050514.41274508 1043760.60380949
1037062.4364218 1030419.35185986 1023830.79801556 1017296.22930473
1010815.10657835 1004386.89703522 998011.07413586 991687.11751788
985414.51291245 979192.75206224 973021.33264048 966899.75817129
960827.53795125 954804.18697213 948829.22584477 942902.18072415
937022.58323558 931189.970402 925403.88457234 919663.87335102
913969.48952842 908320.29101251 902715.84076137 897155.70671678
891639.46173884 886166.68354145 880736.95462884 875349.86223304
870004.99825217 864701.95918978 859440.34609496 854219.76450344
849039.82437945 843900.1400585 838800.330191 833740.01768664
828718.82965962 823736.39737473 818792.35619409 813886.34552477
809018.00876713 804186.9932639 799392.95024997 794635.53480299
789914.40579457 785229.22584226 780579.66126219 775965.38202236
771386.06169666 766841.37741951 762331.00984108 757854.64308332
753411.96469639 749002.66561592 744626.44012069 740282.98579108
735972.00346792 731693.19721211 727446.27426464 723230.9450073
719046.92292384 714893.92456175 710771.6694945 706679.88028438
702618.28244584 698586.60440924 694584.57748527 690611.93582974
686668.41640887 682753.75896512 678867.70598342 675010.00265792
671180.39685915 667378.63910167 663604.48251212 659857.68279775
656137.99821536 652445.18954064 648779.02003798 645139.25543063
641525.66387132 637938.01591322 634376.08448135 630839.64484433
627328.47458654 623842.35358063 620381.06396043 616944.39009419
613532.11855819 610144.03811071 606779.93966637 603439.61627072
600122.86307529 596829.4773129 593559.25827329 590312.00727912
587087.52766222 583885.62474025 580706.10579356 577548.78004244
574413.4586246 571299.95457303 568208.08279402 565137.66004564
562088.50491634 559060.4378039 556053.28089471 553066.85814316
550100.99525151 547155.5196498 544230.2604762 541325.04855751
538439.71638996 535574.09812022 532728.02952667 529901.34800094
527093.89252964 524305.50367636 521536.02356387 518785.29585654
516053.16574306 513339.4799193 510644.08657138 507966.83535904
505307.57739913 502666.16524937 500042.45289227 497436.2957193
494847.55051522 492276.07544261 489721.73002663 487184.37513993
484663.87298775 482160.08709327 479672.88228304 477202.12467267
474747.68165269 472309.42187452 469887.21523671 467480.9328713
465090.44713031 462715.6315725 460356.36095019 458012.51119628
455683.95941147 453370.58385157 451072.26391503 448788.88013054
446520.31414488 444266.44871085 442027.16767535 439802.35596765
437591.89958778 435395.685595 433213.60209653 431045.53823634
428891.38418403 426751.03112399 424624.37124453 422511.29772724
420411.70473646 418325.48740887 416252.54184315 414192.76508991
412146.05514155 410112.31092239 408091.43227887 406083.31996982
404087.87565695 402105.00189532 400134.60212404 398176.58065703
396230.84267385 394297.29421075 392375.8421517 390466.39421961
388568.8589676 386683.14577043 384809.16481597 382946.82709682
381096.04440199 379256.72930868 377428.79517421 375612.15612798
373806.72706352 372012.4236307 370229.16222799 368456.85999475
366695.43480373 364944.80525359 363204.89066145 361475.61105568
359756.88716858 358048.64042932 356350.79295685 354663.26755294
352985.98769524 351318.87753055 349661.861868 348014.86617242
346377.81655778 344750.63978063 343133.26323368 341525.61493947
339927.62354404 338339.21831071 336760.32911395 335190.88643328
333630.82134727 332080.06552759 330538.55123313 329006.21130418
327482.97915672 325968.78877666 324463.57471432 322967.27207879
321479.81653247 320001.14428564 318531.19209106 317069.89723871
315617.19755046 314173.03137495 312737.33758238 311310.05555947
309891.12520444 308480.48692201 307078.08161849 305683.85069694
304297.73605232 302919.68006678 301549.62560492 300187.51600915
298833.29509511 297486.90714704 296148.2969134 294817.4096023
293494.19087715 292178.58685232 290870.54408877 289570.00958986
288276.93079706 286991.25558582 285712.93226145 284441.909555
283178.13661924 281921.5630247 280672.13875563 279429.8142062
278194.54017653 276966.26786892 275744.94888406 274530.53521727
273322.97925479 272122.23377013 270928.25192042 269740.98724285
268560.39365109 267386.42543179 266219.0372411 265058.18410123
263903.82139706 262755.90487273 261614.39062836 260479.23511673
259350.39513998 258227.82784646 257111.49072747 256001.3416141
254897.33867416 253799.440409 252707.60565052 251621.79355809
250541.96361558 249468.07562837 248400.0897204 247337.9663313
246281.66621349 245231.15042932 244186.38034828 243147.31764416
242113.92429234 241086.16256705 240063.99503863 239047.38457089
238036.29431841 237030.687724 236030.52851602 235035.78070585
234046.40858535 233062.37672432 232083.64996803 231110.19343473
230141.97251325 229178.95286052 228221.10039923 227268.38131544
226320.76205624 225378.20932742 224440.69009119 223508.17156387
222580.6212137 221658.00675854 220740.29616372 219827.45763981
218919.45964049 218016.27086039 217117.860233 216224.19692851
215335.2503518 214450.99014034 213571.38616216 212696.40851385
211826.02751855 210960.21372396 210098.93790041 209242.1710389
208389.8843492 207542.04925793 206698.63740668 205859.62065017
205024.97105437 204194.66089471 203368.66265425 202546.94902187
201729.49289055 200916.26735556 200107.24571274 199302.40145679
198501.70827955 197705.14006831 196912.67090416 196124.27506029
195339.92700038 194559.60137699 193783.2730299 193010.91698458
192242.50845054 191478.02281985 190717.43566551 189960.72273998
189207.85997363 188458.82347323 187713.58952048 186972.13457055
186234.43525056 185500.46835819 184770.21086023 184043.63989114
183320.73275168 182601.46690747 181885.81998766 181173.76978354
180465.29424717 179760.37149006 179058.97978186 178361.09754898
177666.70337336 176975.77599114 176288.29429138 175604.23731481
174923.58425256 174246.31444491 173572.4073801 172901.84269306
172234.60016421 171570.6597183 170910.00142318 170252.60548865
169598.45226528 168947.52224326 168299.79605128 167655.25445536
167013.87835775 166375.64879582 165740.54694095 165108.55409743
164479.65170143 163853.82131984 163231.04464929 162611.30351507
161994.57987006 161380.85579376 160770.11349119 160162.33529195
159557.50364917 158955.60113852 158356.61045724 157760.51442316
157167.29597372 156576.938165 155989.42417083 155404.73728176
154822.8609042 154243.77855944 153667.47388279 153093.93062262
152523.13263948 151955.06390522 151389.70850208 150827.05062184
150267.07456495 149709.76473962 149155.10566106 148603.08195053
148053.6783346 147506.87964423 146962.67081404 146421.0368814
145881.96298569 145345.43436748 144811.43636772 144279.95442699
143750.97408466 143224.48097819 142700.46084231 142178.89950827
141659.78290313 141143.09704893 140628.82806205 140116.96215241
139607.48562274 139100.38486791 138595.64637417 138093.25671849
137593.20256779 137095.47067831 136600.04789489 136106.92115029
135616.07746451 135127.50394415 134641.18778168 134157.11625486
133675.27672601 133195.65664142 132718.24353068 132243.02500605
131769.98876182 131299.1225737 130830.41429818 130363.85187194
129899.42331121 129437.1167112 128976.92024546 128518.82216532
128062.81079929 127608.87455247 127157.00190597 126707.18141634
126259.40171502 125813.65150773 125369.91957396 124928.19476638
124488.46601031 124050.72230316 123614.95271391 123181.14638253
122749.29251952 122319.38040531 121891.39938979 121465.33889174
121041.18839839 120618.93746482 120198.57571355 119780.09283396
119363.47858182 118948.72277882 118535.81531206 118124.74613356
117715.50525979 117308.08277121 116902.46881175 116498.6535884
116096.6273707 115696.38049029 115297.90334049 114901.18637578
114506.22011142 114112.99512296 113721.50204579 113331.73157476
112943.67446368 112557.32152493 112172.66362902 111789.69170416
111408.39673585 111028.76976645 110650.8018948 110274.48427576
109899.80811984 109526.76469279 109155.34531519 108785.54136206
108417.34426246 108050.74549913 107685.73660804 107322.30917809
106960.45485064 106600.16531921 106241.43232906 105884.24767682
105528.60321015 105174.49082734 104821.90247699 104470.83015759
104121.26591722 103773.20185318 103426.6301116 103081.54288717
102737.93242271 102395.79100889 102055.11098386 101715.88473292
101378.10468819 101041.76332827 100706.85317792 100373.36680772
100041.29683376 99710.63591733 99381.37676455 99053.51212613
98727.03479699 98401.93761598 98078.21346557 97755.85527153
97434.85600266 97115.20867044 96796.90632878 96479.94207367
96164.30904294 95850.00041595 95537.00941326 95225.32929642
94914.95336761 94605.87496941 94298.08748448 93991.58433532
93686.35898397 93382.40493173 93079.71571891 92778.28492454
92478.10616613 92179.17309937 91881.47941788 91585.01885296
91289.78517333 90995.77218485 90702.9737303 90411.38368907
90120.99597698 89831.80454599 89543.80338394 89256.98651434
88971.34799611 88686.88192333 88403.58242501 88121.44366484
87840.45984099 87560.62518582 87281.9339657 87004.38048073
86727.95906457 86452.66408414 86178.48993947 85905.43106343
85633.4819215 85362.6370116 85092.8908638 84824.23804018
84556.67313456 84290.1907723 84024.7856101 83760.45233577
83497.18566805 83234.98035636 82973.83118064 82713.73295111
82454.68050809 82196.6687218 81939.69249214 81683.7467485
81428.82644957 81174.92658316 80922.04216596 80670.16824339
80419.2998894 80169.43220627 79920.56032443 79672.67940228
79425.78462598 79179.8712093 78934.93439341 78690.96944672
78447.97166469 78205.93636964 77964.85891059 77724.7346631
77485.55902906 77247.32743652 77010.03533956 76773.67821806
76538.25157759 76303.75094921 76070.17188928 75837.50997936
75605.76082597 75374.9200605 75144.98333899 74915.946342
74687.80477445 74460.55436544 74234.19086812 74008.71005951
73784.10774037 73560.37973502 73337.52189123 73115.53008
72894.40019547 72674.12815475 72454.70989778 72236.14138715
72018.41860801 71801.53756787 71585.4942965 71370.28484575
71155.90528943 70942.35172318 70729.6202643 70517.70705163
70306.60824541 70096.32002716 69886.83859949 69678.16018604
69470.28103129 69263.19740044 69056.9055793 68851.40187413
68646.68261153 68442.74413831 68239.58282135 68037.19504746
67835.57722332 67634.72577525 67434.6371492 67235.30781052
67036.73424392 66838.91295331 66641.84046166 66445.51331094
66249.92806194 66055.08129418 65860.96960581 65667.58961344
65474.93795208 65283.011275 65091.8062536 64901.31957734
64711.54795357 64522.48810749 64334.13678196 64146.49073745
63959.5467519 63773.30162064 63587.75215623 63402.8951884
63218.72756395 63035.24614658 62852.44781688 62670.32947212
62488.88802624 62308.1204097 62128.02356938 61948.59446848
61769.83008643 61591.72741879 61414.28347714 61237.49528898
61061.35989765 60885.87436221 60711.03575738 60536.84117338
60363.2877159 60190.37250598 60018.09267992 59846.44538915
59675.42780022 59505.03709461 59335.27046873 59166.12513375
58997.59831558 58829.68725472 58662.3892062 58495.70143952
58329.62123849 58164.14590122 57999.27273999 57834.99908117
57671.32226513 57508.23964618 57345.74859246 57183.84648587
57022.530722 56861.79870999 56701.64787253 56542.07564572
56383.07947901 56224.65683513 56066.80518997 55909.52203255
55752.80486493 55596.65120209 55441.05857193 55286.0245151
55131.54658501 54977.62234771 54824.2493818 54671.42527839
54519.14764102 54367.41408557 54216.22224019 54065.56974524
53915.45425321 53765.87342864 53616.82494805 53468.30649991
53320.3157845 53172.85051388 53025.90841183 52879.48721376
52733.58466664 52588.19852895 52443.3265706 52298.96657287
52155.11632833 52011.77364078 51868.93632519 51726.60220764
51584.76912523 51443.43492605 51302.59746908 51162.25462415
51022.40427187 50883.04430356 50744.17262121 50605.7871374
50467.88577523 50330.46646827 50193.52716052 50057.0658063
49921.08037024 49785.56882718 49650.52916216 49515.95937029
49381.85745678 49248.22143678 49115.04933543 48982.33918772
48850.08903845 48718.29694223 48586.96096333 48456.0791757
48325.6496629 48195.670518 48066.13984359 47937.05575167
47808.41636363 47680.21981019 47552.46423134 47425.14777629
47298.26860342 47171.82488021 47045.81478321 46920.236498
46795.0882191 46670.36814995 46546.07450283 46422.20549884
46298.75936784 46175.73434841 46053.12868775 45930.94064172
45809.1684747 45687.8104596 45566.86487779 45446.33001908
45326.2041816 45206.48567185 45087.17280458 44968.26390277
44849.7572976 44731.65132835 44613.94434244 44496.63469529
44379.72075034 44263.200879 44147.07346055 44031.33688217
43915.98953885 43801.02983337 43686.45617623 43572.26698563
43458.46068741 43345.03571502 43231.99050949 43119.32351936
43007.03320063 42895.11801678 42783.57643865 42672.40694446
42561.60801973 42451.17815727 42341.11585711 42231.41962647
42122.08797976 42013.11943846 41904.51253116 41796.26579346
41688.37776799 41580.84700431 41473.67205891 41366.85149518
41260.38388333 41154.2678004 41048.50183018 40943.0845632
40838.01459671 40733.29053459 40628.91098735 40524.87457209
40421.17991246 40317.82563864 40214.81038727 40112.13280145
40009.79153066 39907.7852308 39806.11256407 39704.772199
39603.76281038 39503.08307924 39402.73169281 39302.7073445
39203.00873384 39103.63456647 39004.58355412 38905.85441452
38807.44587143 38709.35665459 38611.58549964 38514.13114818
38416.99234764 38320.16785133 38223.65641835 38127.45681359
38031.56780769 37935.98817702 37840.71670362 37745.75217521
37651.09338512 37556.73913229 37462.68822122 37368.93946195
37275.49167004 37182.34366653 37089.49427788 36996.94233602
36904.68667823 36812.72614717 36721.05959085 36629.68586257
36538.6038209 36447.81232968 36357.31025797 36267.09648002
36177.16987526 36087.52932823 35998.17372861 35909.10197117
35820.31295572 35731.80558712 35643.57877522 35555.63143486
35467.96248583 35380.57085286 35293.45546557 35206.61525845
35120.04917085 35033.75614694 34947.7351357 34861.98509087
34776.50497096 34691.29373917 34606.35036343 34521.67381634
34437.26307515 34353.11712172 34269.23494254 34185.61552865
34102.25787568 34019.16098376 33936.32385754 33853.74550615
33771.42494319 33689.36118668 33607.55325907 33526.00018721
33444.70100229 33363.65473987 33282.86043983 33202.31714634
33122.02390786 33041.97977711 32962.18381104 32882.63507081
32803.33262179 32724.27553348 32645.46287958 32566.89373788
32488.56719029 32410.48232281 32332.63822549 32255.03399244
32177.66872179 32100.54151565 32023.65148014 31946.99772532
31870.57936522 31794.39551776 31718.44530477 31642.72785197
31567.24228893 31491.98774907 31416.96336964 31342.16829167
31267.60165998 31193.26262318 31119.15033359 31045.26394728
30971.602624 30898.16552723 30824.95182409 30751.96068535
30679.19128543 30606.64280234 30534.31441772 30462.20531676
30390.31468822 30318.64172441 30247.18562114 30175.94557776
30104.92079708 30034.11048539 29963.51385245 29893.13011144]
de machtconstante is -3.5236716254510267
Als je de simulatie een aantal keer uitvoert, zie je dat er een structureel verschil zit tussen de waarde die je verwacht en de waarde die je uit de fit krijgt.
Eerste hoofdwet¶
Een goede controle voor ons simulatie is te onderzoeken of deze voldoet aan de eerste hoofdwet van de thermodynamica.
Omdat er (nog) geen warmte wordt uitgewisseld, moeten we alleen de hoeveelheid arbeid bepalen. De arbeid wordt geleverd door de zuiger(s) op de deeltjes. Dus in de functie waar we de botsingen van de deeltjes met de wand behandelen, kunnen we ook de verrichte arbeid bepalen. De arbeid wordt gegeven\ door:
Hier staan de 1 voor de begintoestand en de 2 voor de eindtoestand van het proces. In de differentiaalvorm wordt dit geschreven als:
Zoals dat in het boek gebeurt, kiezen we voor de notatie met in plaats van , om aan te geven dat deze integraal afhankelijk is van het gekozen proces en niet alleen van het begin- en eindpunt. Voor ons experiment, waar het volume niet geleidelijk maar in stapjes verandert, kan je de hoeveelheid arbeid per tijdstip dus herschrijven tot:
waarin de hoogte is van de zuiger.
Dit betekent dat we de arbeid per tijdstap in onze code kunnen bepalen op hetzelfde moment dat we de druk bepalen met behulp van de gezamenlijke stoot van de moleculen. Daarvoor moeten we de code voor de functie left_right_collision nogmaals aanpassen:
def left_right_collision(particle: ParticleClass):
""" verzorgen van botsingen met wand links en rechts, die als zuiger kunnen bewegen """
global box_length, v_piston, impulse_outward, work
if abs(particle.r[0]) + particle.R > box_length / 2:
particle.r[0] = np.sign(particle.r[0]) * (box_length/2 - particle.R)
piston_velocity = np.sign(particle.r[0]) * v_piston
relative_velocity = particle.v[0] - piston_velocity # stelsel zuiger
particle.v[0] = -relative_velocity + piston_velocity # stelsel waarnemer
impulse_outward += 2 * particle.m * abs(relative_velocity)
work += 2 * particle.m * relative_velocity * piston_velocityVoordat we een simulatie draaien is het goed even stil te staan en na te denken. Net als bij een experiment in het lab doorlopen we bewust de onderzoekscyclus. We moeten daarvoor eerst een hypothese opstellen en daarna pas de simulatie uitvoeren.
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.