Die wesentliche Physik des Modells in
algorithmischer Form
Die folgenden, auf das wesentliche vereinfachten, Prozeduren geben kom-
mentierte Listings der verwendeten Algorithmen:
(Beachte: Die Koordinate x läuft hier vertikal: +x hinauf; y horizontal,
+y nach rechts)
procedure move_mol; {At every time step add the velocity }
var i : integer; { component to the x and y coordinates}
begin { here: x vertical, y horizontal }
for i:=1 to ib do begin { ib=256, number of molecules }
x[i]:=x[i]+vx[i]; { add vx[i] to x[i] }
y[i]:=y[i]+vy[i]; { add vy[i] to y[i] }
if x[i]<0 then begin { reflect at bottom wall, if collision}
x[i] :=-x[i]; { reflect x-component }
vx[i]:=-vx[i]; { reverse sign of velocity component }
pu:=pu+2*abs(vx[i]) { sum change of momentum for pressure }
end; { at bottom }
if x[i]>rb then begin { reflect at upper wall, if collision }
x[i]:=rb+rb+1-x[i]; { rb is upper boundary - 1 }
vx[i]:=-vx[i];
po:=po+abs(vx[i]); { first part of momentum change }
if pistflg=1 then { add momentum of piston}
vx[i]:=vx[i]+165*signad; { gives correct adiabate for Cp/Cv=5/3}
po:=po+abs(vx[i]) { second part of mom. change incl.pist}
end;
if y[i]<0 then begin { reflect at left wall, if collision }
y[i] :=-y[i];
vy[i]:=-vy[i]
end;
if y[i]>ub then begin { reflect at right wall, if collision }
y[i] :=ubb-y[i]; { ub=boundary - 1, ubb=2*ub+1 }
vy[i]:=-vy[i]
end;
if (((gal>0) and (hi(x[i])>5)) or { grav. acceleration postfix ! }
((gal<0) and (hi(x[i])<74))) { turn off for bottom/top 6 rows }
then vx[i]:=vx[i]-gal-gal; { add acceleration to component }
end
end; {move_mol}
procedure hit90(xc,yc: integer); {collision between i and j: }
var j : integer; { if i,j have the same grid- }
halfsum,temp : integer; { coordinates but not necessarily }
begin { the same pixel-coordinates. }
j:=net[xc,yc]; {j-particle is already at xc,yc }
halfsum:=round(0.5*(vx[i]+vy[i]+vx[j]+vy[j])); {floating operation! }
temp :=halfsum-vy[i]; {Determine the velocities after the}
vy[i]:=halfsum-vx[j]; { collision, which have to conserve}
vx[j]:=halfsum-vy[j]; { kinetic energy and the total mo- }
vy[j]:=halfsum-vx[i]; { mentum of each velocity component}
vx[i]:=temp { This algorithm is valid for a }
end; {hit90} { fixed scattering angle of 90ø. }
procedure plot_mol; {plot a molecule at its newly com- }
var m,n : integer; { puted place. If this is occupied,}
begin { we have a collision }
net:=net0; {mark the 50x80 cells as unoccupied}
for i:=1 to ib do begin {for all 256 molecules do... }
m:=198-hi(x[i] shl 1); {draw in frame of 100x160, thus }
n:=hi(y[i] shl 1)+1; { keeping pixel accuracy }
xc:=hi(x[i]); yc:=hi(y[i]); {get coordinates in 50x80 grid: }
if net[xc,yc]=0 then begin {if cell free, put i-particle in, }
net[xc,yc]:=i; {else we have a collision ¾¾¾× }
if (i=1) and (hitflg=0) { ½ }
then putpixel(n,m,15) {#1 molecule ½ }
else begin { ½ }
if (difflg=1) and (xcol[i]>=10240) {diffusion on ½ }
then color:=13 else color:=col; { ½ }
if col>0 then putpixel(n,m,color) {show molecule ½ }
end { ½ }
end { ½ }
else begin {compute collision velocities ¿ }
if wflg=1 then hit(xc,yc); {between i,j }
inc(cct); {update collision-counter }
if hitflg=1 then
putpixel(n,m,13) {show collision cell if switch on }
end
end
end; {plot_mol}
procedure hitrand(i,j:integer); {Scattering angle 0 < b < p}
var cx,cy,wx,wy,cs,si,beta: extended; {beta may be fixed <>90ø,too}
begin {Derivation see in Expl. (8)}
beta:=random*p; {random scattering angle }
cs:=cos(beta); si:=sin(beta); {trigonometric functions }
cx:=vx[i]+vx[j]; {Components of the unchanged}
cy:=vy[i]+vy[j]; { velocity of centr.of grav.}
wx:=(vx[j]-vx[i])*cs-(vy[j]-vy[i])*si; {Rotate relative velocity }
wy:=(vx[j]-vx[i])*si+(vy[j]-vy[i])*cs; { by scattering angle beta. }
vx[i]:=round(0.5*(cx-wx)); {Compute components of }
vx[j]:=vx[i]+round(wx); { scattered relative }
vy[i]:=round(0.5*(cy-wy)); { velocity. }
vy[j]:=vy[i]+round(wy); {This algorithm is valid }
end; {hitrand} { for any scattering angle. }
Flags: pistflg Kolbenstempel: 0 Ausgangslage, >0 weg davon
difflg Diffusionsflag. Wenn an, Moleküle oben magenta,sonst gelb
wflg Kollisionen ein/aus
hitflg Kollisionsorte werden angezeigt oder nicht
net[80,50] ist das Gitter, in dem mit i oder 0 die An/Abwesenheit
des Moleküls i markiert wird. Gezeichnet werden die Mole-
küle jedoch im Feld von 160x100 Pixel. Innerhalb des 2x2
Quadrats wird jedes Molekül exakt gezeichnet. Ein Zusam-
menstoss von i,j findet in allen folgenden Konfiguratio-
nen statt:
| Doppelbesetzung des gleichen Pixels |
|
Doppelbesetzung der gleichen Zelle |
Siehe auch Programm COLLIS90, welches die verwendeten Stoss- und Reflek-
tionsgesetze einzeln vorführt, sowie Erkl.(8) für die Stossphysik.
Zur Beachtung: Die algorithmische Behandlung der Stossereignisse enthält
eine Vereinfachung. Es wird immer nur am Ende eines Zeitinkrements (Bild-
wechsel Dt) abgefragt, ob eine Kollision vorliegt. Der Stoss zweier Par-
tikel, welche im Flug gleichzeitig eine Raumzelle passieren, wird nicht
erfasst. Um das zu leisten, müssten Kollisionzeiten auf Bruchteile eines
Ticks (=Zeitelement der Simulation) ermittelt werden, was erheblich mehr
Rechenaufwand verursachte. Im Grenzfall Dt -> dt geht unser Algorithmus
in den exakten über.
Diese Näherung bewirkt nur an einer Stelle eine Abweichung von der Beob-
achtung: Die Abhängigkeit der Stosszahl von der Geschwindigkeit wird bei
grösserer Geschwindigkeit der Teilchen immer weniger genau. Wir beschrän-
ken uns daher beim Tutorial DIFFUS darauf, die (korrekt simulierte)
Temperaturabhängigkeit der Diffusionskonstanten zu bestimmen und verzich-
ten darauf, ungenaue Diffusionszeiten anzugeben.
In move_mol oben wurde die recht komplizierte Behandlung der Reflektion
an der teilweise geöffneten Mittelwand weggelassen, da diese nichts zum
physikalischen Verständnis beiträgt. Auf dieser Diskette ist der Source-
Code als GASSIM.PAS beigelegt. Sie können dort die vollständigen Proze-
duren studieren.
Wenn die Mittelwand geschlossen oder nur teilweise offen ist und gleich-
zeitig mit 'T' oder 'A' der Kolben bis zur maximalen Kompression einge-
fahren wird, so ist im oberen Volumen der Abstand Stempel-Mittelwand
klein geworden. Bei raschen Molekülen kann y+vy*dt dann grösser als die-
ser Abstand werden. In diesem Fall entstehen Mehrfachreflektionen der Mo-
leküle an der Wand. Der Algorithmus behandelt nur Dreifachreflektionen.
Deshalb kann es in seltenen Fällen vorkommen, dass Moleküle bei obiger
Sachlage durch die geschlossene Wand hindurch'tunneln':
