F?r vi g?r l?s p? simuleringen er det noen begreper, analytisk og numerisk l?sning, som vi kommer til ? bruke en del. En analytisk l?sning er n?r vi kommer fram til et svar ved ? gj?re utregninger for h?nd uten tanke p? presisjon. En numerisk l?sning er n?r vi lager et program som kan g? av seg selv, men finner et tiln?rmet svar som vi bruker videre. Ved ? bruke smarte algoritmer, sekvenser av forh?ndsbestemte operasjoner, og metoder kan man f?r denne tiln?rmingen til ? komme s? n?rme den faktiske l?sningen at forskjellen blir ubetydelig for videre utregninger. Selvf?lgelig m? man alltid ha presisjon i bakhuet, dette vil bli nevnt n?r det blir relevant.
Du har sikkert h?rt om Python og Java, dette er det vi kaller for kodespr?k, som betyr at de fungerer litt forskjellig. Det er alt fra spr?ket, hvordan og hva de aksepterer og m?ten ting skal bli skrevet p?, til hastighet p? hvor lang tid de bruker p? ? fullf?re en kode. Vi vil bruke Python til v?re koder, for det er bedre egnet til den type kode vi vil lage.
Vi har gjort et par viktige antagelser som forenkler systemet og gj?r modellen lettere ? programmere.
- Partiklene i boksen oppf?rer seg som en ideell gass, alts? all kollisjon er elastisk.
- Symmetrien til H2 er en enkelt kule, slik at man ikke m? ta hensyn til posisjonene til hvert av atomene gjennom simuleringen.
- Det er ingen ytre og indre krefter som virker p? systemet, slik som magnetiske felt, tiltrekning mellom partiklene p? grunn av ladningsforskjeller og tyngdekraft osv.
- Partiklene kan ikke kollidere med hverandre, og vil dermed g? rett igjennom hverandre eller okkupere samme posisjon.
- Det skal v?re s? tilfeldig som mulig gitt to modeller, gaussisk fordeling for fartskomponentene, og uniform fordeling for posisjonene.
- Partiklene kan bare kollidere med veggene.
- Antall partikler i boksen skal v?re konstant.
Hensikten bak simuleringen av partiklene er at det blir v?res rakett motor, men hvordan kan dette fungere som en rakett motor, hvordan f?r man en netto kraft som virker i den retningen vi ?nsker?
Her kommer sannsynlighetsfordelingene inn, p? grunn av den gaussiske fordelingen og v?r antagelse vi har en ideell gass s? betyr det at energien i systemet alltid er bevart. Oppdriften vil komme av bevegelsesmengden fra partiklene som treffer veggene, men summen av denne bevegelsesmengden p? hele boksen forblir lik 0 s? lenge det ikke er noe hull som det kan slippes ut fra. Husk at fartskomponentene er tilfeldig fordelt fra den gaussiske modellen. Dersom det blir uklart i hvorfor dette har s? mye ? si s? anbefaler jeg ? g? tilbake til en tidligere post som dekket dette. La oss n? introdusere et hull i bunnen av boksen som partiklene kan slippes ut fra, dette gj?r at ikke all bevegelsesmengden lenger motvirkes i hver sin z-retning. Tapet av bevegelsesmengden ut av hullet vil v?re det samme som vil virke opp i taket, som gir en bevegelse i den retningen. Denne kraften som da virker opp vil v?re det som f?r raketten til ? g? oppover.
Vi har noen parametere som er bestemt p? forh?nd slik som antall partikler og tiden simuleringen dekker. I simuleringen vil vi bruke \(N = 10^5\) partikler, over en tid p? \(t = 10^{-9}\) sekunder. V?re justerbare parametere er temperaturen gitt i kelvin, lengden p? hver side av boksen med enhet \(L = 10^{-6}m\), som er 1000 ganger mindre enn en millimeter, og arealet p? hullet, der hullet er et kvadrat. I tabellen under er v?re startverdier satt opp, noen vil bli endret p? senere.
| N | K | hull |
x |
y | z |
|---|---|---|---|---|---|
| \(10^5\) | 3000 | \(0.25L^2\) | \(1L\) | \(1L\) | \(1L\) |
Vi vet allerede tiden som simuleringen vil dekke,\(10^{-9}\) sekunder, det betyr at vi kan finne et uttrykk for kraften \(F\) ut fra impulsen. Vi kan se dette som en impuls p? grunn av at vi ikke regner ut summen av bevegelsesmengden f?r simuleringen er ferdig, som betyr at vi kan si at bevegelsesmengden virker p? boksen i et kort tidsrom, dermed kan vi skrive kraften om fra impuls slik \(Ft = \Delta p \Rightarrow F = \dfrac {\Delta p}{t}\). \(\Delta p = mv_1 - mv_0\) der \(mv_0 = 0\) som er f?r vi har startet simulasjonen, og \(mv_1 = m\sum v_{total}\) som er total bevegelsesmengde fra alle partiklene som forlot boksen under hele simulasjonen.
Kan du tenke deg hvorfor vi kan sette \(m\) utenfor summeringstegnet \(\sum\)?
Fartskomponentene blir bestemt ut fra en gaussisk fordeling, men posisjonene er bestemt av uniform fordeling. Uniform fordeling betyr at det er like stor sannsynlighet at en partikkel blir plassert i midten av boksen som et hvilket som helst annet sted i boksen. Et eksempel som kan visualisere dette bedre er om vi deler hver side opp i 20 like store deler for s? ? skulle kaste en 20 siders terning, og tallet man f?r bestemmer posisjonen i den gitte retningen. Denne prosessen gj?res s? en gang for hver retning og hvert kast er uavhengig av hverandre, dermed f?r man en uniform fordeling p? posisjonene i boksen. Dette gj?res 3 ganger for hver partikkel og med b?de fartskomponentene og posisjonene p? plass er vi nesten klare til ? starte simuleringen.
Posisjon og fart er p? plass, s? kommer vi til hvordan vi f?r disse til ? bevege seg. P? grunn av noen begrensninger kan man ikke gj?re mer enn en ting om gangen i Python, det betyr at n?r vi skal bevege partiklene s? m? dette gj?res komponent etter komponent, partikkel etter partikkel. Som du kanskje innser s? blir dette en lang prosess for \(N = 10^5\) partikler. M?ten vi f?r de til ? bevege seg er gjennom ? bruke noen av de grunnleggende formlene for fart, strekning og tid. Farten i en gitt retning har vi, og tiden vil v?re det vi kaller for et tidssteg.
S? hva er poenget med et tidssteg? Under de fleste numeriske m?ter ? simulere ting p?, slik som n?, buker man steglengden for ? justere p? presisjonen. Mindre steglengde betyr st?rre presisjon p? bekostning av tiden det tar ? kj?re koden p? grunn av at den m? gj?re samme prosess flere ganger for ? dekke samme tidsramme. Likes? er det for st?rre steglengde, st?rre steg betyr mindre presisjon, men raskere utf?ring av koden. I v?rt tilfelle vil vi bruke tidssteg p? \(\Delta t =dt = 10^{-12} \), med litt enkel regning f?r vi at antall steg er lik 1000. Det betyr at vi skal flytte p? hver partikkel 1000 ganger. Vi kan skrive dette som \(posisjon_1 = posisjon_0 + v \cdot dt\). Vi vil ogs? bruke noe som kalles en "loop" for ? f? dette til ? g? av seg selv, som vil g? igjennom hver partikkel og oppdatere, eller gi, dem en ny posisjon. N?r dette har skjedd og den er ferdig med den siste partikkelen g?r den tilbake til den f?rste partikkelen og gj?r det samme igjen. Dette kalles en iterasjon, en repetitive handling, som gj?res 1000 ganger. Dette kalles for Euler's metode, og vi kunne ha brukt mange andre variasjoner som kunne ha ulik presisjon, men for v?rt tilfelle s? holder det lenge med den standard versjonen. Hadde vi hatt med akselerasjon ville ikke denne fungere, og vi ville ha brukt en variasjon som heter Euler Cromer. Denne kommer vi til litt senere s? vi kan vente med ? forklare den til da.
Det er p? tide ? p?peke hva vi gj?r for ? oppdage om en partikkel har kollidert i en vegg. Veggkollisjon blir oppdaget dersom den nye posisjonen er utenfor en av sidene. Ta for eksempel i x-retning, origo er plassert midt p? x, som betyr at lengden x i positiv og negativ retning er \(\pm \dfrac{1}{2}x \) der \(x = 1L\). P? bilde under er det et eksempel p? en kollisjon med veggen. Her vil den nye posisjonen bli justert slik at det g?r tilbake i motsatt retning, s? partikkelen holder seg inne i boksen, samtidig som fartskomponenten i den retningen vil endre retning. Det er andre m?ter ? gj?re dette p?, som for eksempel ? bare se om det er utenfor og flytte den slik at de er akkurat inne i boksen, inkludert ? endre fortegnet p? farten i den retningen. Denne m?ten ? implementere veggkollisjon er definitivt lettere, men man mister hele lengden som partikkelen tilbakela utenfor boksen. I v?r simulering blir dette ogs? tatt hensyn til, og partikkelen kan v?re s? langt unna som den vil, og fremdeles bli plassert inne i boksen, med fortegnet til fartskomponenten endret med hensyn p? hvilken retning den ville ha fortsatt i. Det betyr at vi har mulighet ? f? partikler med flere veggkollisjoner under samme tidssteg. Dette g?r dessverre p? bekostning av tiden det tar ? kj?re koden p? grunn av de ekstra beregningene som m? gj?res for ? f? dette til ? fungere.
Hva om en partikkel g?r gjennom hullet, hva gj?r vi med dem? For veggkollisjon brukte vi noe vi kaller for reflective boundary condition, for partikler som g?r ut av hullet bruker vi noe vi kaller for periodic boundary condition. Man kan tenke p? det som om man teleporterer partikkelen tilbake inn i boksen fra toppen av med en avstand fra taket lik den avstanden den tilbakela utenfor boksen. Med andre ord, vi trekker fra den nye posisjonen med en lengde lik h?yden av hele boksen. I praksis blir det egentlig ? legge til fordi posisjonen vil alltid ha et negativt fortegn p? grunn av plasseringen av hullet. Dette er en metode som er moderat n?r det kommer til hvordan det vil p?virke farten p? koden. Enklere metoder som ? bare la hullet v?re som en vegg og bare registrere n?r en partikkel treffer innenfor hadde v?rt lettere, men det vil ikke akkurat kunne simulere helt det samme som skjer under p?fyll av drivstoff til motoren.
Det er ogs? bedre enn ? bare plassere dem i en fast z verdi for da mister man noe av tilfeldighetene man satt helt i starten, for da f?r man et slags system som vil motvirke entropien i systemet. Den m?ten vi valgte vil b?de kunne simulere det som skjer i motoren, n?r nye partikler blir sendt inn for ? erstatte det som ble mistet, og bevare tilfeldigheten i systemet p? en bedre m?te. Alt vi trenger ? gj?re er ? telle antall som g?r ut gjennom hullet, samle deres vz og endre p? pz.
Med det er vi nesten klare, bare en ting til og det er oppsamling av dataene. N?r simulasjonen er ferdig vil den returnere en liste med f?lgende verdier og p? den f?lgende formen.
|
|
|---|
S? hva er denne tabellen med mindre tabeller i seg? Dette er en kombinasjon av to m?ter ? lagre verdier p?, den f?rste som inneholder begge tabellene kalles en liste. Den andre typen som hver av de mindre tabellene er kalles en array. "Hva er da forskjellen" lurer du nok sikkert p? n?. En liste kan inneholde mange ulike typer objekter, deriblant andre lister med ting, en array kan bare ha tall verdier, men er i motsetning til lister enklere ? bruke dersom man skal regne ut nye verdier, inkludert mye raskere. N?r vi leser av resultatet er vi hovedsakelig bare interessert i de to siste verdiene, resten er dersom vi trenger ? sammenligne dersom vi skulle trenge ? endre p? verdiene.
- N = antall partikler brukt i simulasjonen
- k = temperatur i boksen
- A = arealet av hullet
- s = seed, en verdi brukt under fordeling av fartskomponentene og posisjonene som gj?r at vi f?r samme fordeling hver gang vi kj?rer koden s? lenge denne verdien er det samme. I v?r simulering brukte vi verdien seed = 75.
- xl = lengden p? boksen langs x-aksen gitt i L = 10-6 m
- yl = lengden p? boksen langs y-aksen gitt i L = 10-6 m
- zl = lengden p? boksen langs z-aksen gitt i L = 10-6 m
- ntot = antall partikler som gikk ut gjennom hullet gjennom hele simulasjonen.
- Vtot = samlet mengde av absolutt verdien av vz fra partiklene som forlot boksen gjennom hele simulasjonen. Trenger bare ? samle vz fordi det bare er en type partikkel som gj?r at ? legge sammen massen ganget med individuelle fartskomponenter blir det samme som ? legge sammen komponentene f?rst, for s? ? gange det med massen.
Brukte disse parameterne da vi kj?rte koden f?rste gang
| N = 105 | k = 3000 | A = 0.25L2 | x = 1L | y = 1L | z = 1L |
|---|
og fikk f?lgende resultat
| ntot = 4.3076 * 104 | Vtot = 2.00004018 * 108 m/s |
|---|
N? gjelder det ? se hvor godt dette stemmer med det vi forventer. For ? sjekke dette bruker formelen for den midlere farten som partiklene vil ha som er \(\langle v \rangle = \int_{0}^{\infty} v P\left(v\right) \, dv\), har ikke tenkt til ? ta utledningen her, men du kan g? hit dersom du er interessert i utledningen. Dette integralet blir s? \(\langle v \rangle = \sqrt{\dfrac{8 kT}{\pi m}}\).
Her legger vi inn verdiene for den analytiske l?sningen og bruker at \(k = 1.38 \cdot 10^{-23} J/K\) som er Boltzmannkonstanten, og massen for H2 blir \(m_{H_2} = 2.018 \cdot 1.66 \cdot 10^{-27}\) kg, setter inn verdiene og f?r \(\langle v \rangle = 5612.68 m/s \). Tar s? ? deler \(\langle v_{numerisk} \rangle = V_{tot}/n_{tot} \) som blir 4643 m/s, dette gir oss en for stor prosent del n?r vi ser p? presisjonen mellom det den ga, og det som er ? forvente.
M?ten vi regner ut det er \(\dfrac{absolutverdi(numerisk - analytisk)}{analytisk}\) i v?rt tilfelle er dette p? 17.3%, som betyr at verdiene er rundt 17.3% for h?ye.
Det betyr at noe ikke fungerte som det skulle, som er mest sannsynlig noe ? gj?re med hvordan vi gikk lagret verdiene, eller andre sm?feil i koden. Mye tyder p? at det er tekniske feil som avrundingsfeil gjort under veggkollisjonene og andre sm?ting som er vanskelig ? fikse.
Noe som er verdt ? merke seg er at 2 versjoner av koden ble laget, eneste forskjell mellom dem var m?ten vi lagret verdiene p?, men de ga ulikt resultat. Vi fortsatte med denne versjonen som hadde endret p? m?ten vi lagret verdiene p?, for det var et fokus p? at det skulle v?re en pc-vennlig kode ? kj?re, den andre kunne fort f?re til at pc'en gikk tom for minne selv om den f?rste var mer presis ut fra ? gj?re en test. Det kan fremdeles ha v?rt like d?rlig presisjon, men presisjonen er ikke det som var hovedfokuset.
Dette er en typisk ting som er vanskelig ? forutse n?r det kommer numeriske l?sninger, det er alltid et forhold mellom presisjon, utf?rings tid og bruk av ressurser for en pc, og denne gangen var det presisjonen som fikk st?yten.
Da har vi det vi trenger for ? skyte opp raketten v?r. Fortsetter i neste blogg-post, eller g? til forrige blogg-post.
Under er det en video for moro skyld som viser boksen med alle de 100 000 partiklene, den ble laget med hjelp av en modifisert versjon av koden, visualisert ved hjelp av et program som heter ovito, og er en tidligere versjon, s? den har ikke implementert hullet i simulasjonen. UIO sin nettside(denne) var ikke helt forn?yd med ? skulle vise videoen s? det kan v?re litt hakkete.