Nå som vi har gått over mesteparten av teorien i den forrgie bloggen, kan vi begynne å simulere! Vi har lyst å simulere, for da kan vi veldig lett justere parameterene (konstantene) som temperatur, hvilken gass vi har i boksen, og etter hvert, størrelsen på hullet
For å simulere en gas i en boks, må vi ta utgangspunkt i hvordan partikellene burde oppføre seg. Siden vi antar at det er en ideel gas, antar vi elastiske kollisjoner (se forrgie bloggpost), og trenger dermed kun å implementere kollisjoner med veggene. Når en partikkel treffer en vegg (gitt elastiske kollisjoner), vil vinkelen den kommer ut med være den same som vinkelen den treffer veggen med, på nøyaktig same måte som ett speil. Dette er lettere å forstå hvis vi først ser på det i 2 dimensjoner:
/bilder/momentumforklaring.jpg)
På diagrammet ser vi en kollisjon mellom en partikkel og veggen. Vi ser her at, gitt en boks med vegger parralelt med aksene, hastigheten i y retning vil være uendret. Dette er også logisk, ettersom kraften kun vil virke i x retning. Siden bevegelsesmengde er masse ganget med hastighet, og massen er konstant, kan vi også konkludere med at endringen i x retning er 2px der px er den originale bevegelsesmengden til partikellet i x retning (se gas, gas, gas! bloggposten for forklaring). For 3 dimensjoner gjelder dette også, siden kraften kun virker i en retning av gangen, snur vi kun på partikkelet i en retning.
Vi lagrer posisjon og hastighet som arrays, altså en type liste, med N elementer hvor hver av elementene har 3 elementer til, der N er antall partikeller, og subelementene er komponentene til posisjonen til hver av de. Den vil se slik ut: \([[x_1, y_1, z_1], [x_2, y_2, z_2], ... ] \). Så når vi vil ha ut informasjon om partikkel 1, trekker vi ut den første listen med posisjoner ( \( [x_1, y_1, z_1] \) ). Du kan tenke på det som mange små kuler i en sekk, der hver av kulene har informasjonen til en partikkel. Posisjonen blir generert tilfeldig innenfor boksen, vår framtidig motor, med uniform fordeling. Dette vil si at partiklene har like stor sannsynlighet for å starte hvor som helst i boksen. Når vi har flere titusener av partikeller vil vi ende opp med en fin uniform fordeling gitt store talls lov og CLT. Dette er tilsvarende det vi forventer å se i virkeligheten, ettersom det betyr at gassen er gjevnt fordelt, samt trykket. Vi lagrer også hastighet tilsvarende. Vi bruker formelene vi fant i forrige bloggpost til å tilfeldig generere hastigheter i hver dimensjon for hver partikkel.
Vi definerer et kordinatsystem slik at midten av boksen ligger i origo, altså (0, 0, 0). Dette lar oss sjekke komponentvis om avstanden i x, y, eller z retning er over halvparten av lengden til boksen. I det tilfelle vil partikellen være utenfor boksen, og må bli snudd.
Koden vår er vektorisert, som er mye raskere enn gjentatt iterasjon. Vektorisering i dette tilfelle betyr at vi kan utføre endringer direkte på listen (sekken) og det vil skje på alle elementene (kulene) i listen. Vi kan opptatere alle posisjonene samtidig, ved å gjøre følgende:
\(\vec{X} = \vec{X} + \vec{v} \times dt\)
der \( \vec{X} \) er listen (sekken) med alle posisjonene og \( \vec{v} \) er listen med alle hastighetene. Dette er den normale definisjonen for posisjon og bevegelse, men i vektorform der vi legger til alle posisjonsendringene i alle dimensjonene samtidig. Vi kan også få koden til å kjøre enda raskere ved å vektorisere logikken vi bruker for veggkollisjoner. Numpy har en fantastisk funksjon som lar oss lage en ny array med "True" og "False" verdier for om partikellene er utenfor boksen eller ikke. I tillegg så kan maskiner jobbe veldig fint med True og False verdier, og behandler de som 1 og 0 verdier. Denne nye listen kaller vi for Noutside. Dermed kan vi snu alle partikellene utenfor boksen, samtidig, ved denne logikken
\(\vec{V} = \vec{V} -2 \times \vec{V} \times \vec{Noutside}\)
Legg merke til at når Noutside er 0 i et element, så blir dette til at V = V i den dimensjonen for den partikellen, ettersom man ganger endringen med 0. Når Noutside er 1 i et element, så blir dette V = -V for den dimensjonen for den partikellen. Dermed klarer vi å snu alle retninger på hastighetskomponenetene til partikeller som er utenfor boksen, i de riktige dimensjonene, på denne ene linjen, gitt at vi har laget Noutside først. Dette tok koden vår fra flere minutter kjøretid, til under 2 sekunder!
Vi må lage en boks i bunden av raketten for gassen til å rømme ut av. Dette gjør vi også ved å vektorisere logikken, som tidligere. Dersom posisjonen er innenfor et kvadratisk området i bunden av boksen vår, med størrelse 0.25L2 , så teller vi den vertikale bevegelsesmengden til partikellen og at en partikkel har unslippet. I praksis vil forsvundene partikeller bli erstattet med nye partikeller fra drivstoffet vårt. I koden vår, vil vi simplifisere dette ved å kun la den sprette av åpningen, og fortsette. Dette kan føre til unøyaktigheter, siden nye hastigheter ikke er tilfeldige, men den samme som forrige partikkel, og de begynner i bunnen av boksen istedenfor toppen. For eksempel, kan en partikell sprette rundt i en loop og treffe åpningen på et fast interval, og telle som en ny en hver gang. Når det er sagt, så kjører vi koden vår over et veldig kort tidsinterval (i størrelsesorden med et nanosekund), dermed er det tvilsomt at partikellene vil ha tid til å treffe den flere ganger. I tillegg så bidrar dette med en bit av nøytaktighet til systemet, som er at totalenergien til gassen vil holde seg konstant; den kinetiske energien vil aldri endre seg, ettersom masse og hastighet er konstanter med denne metoden. Dermed, selvom dette ikke er representativt av virkeligheten, så er det usannsynlig at dette vil påvirke resultatene våre.
Deretter setter vi koden vår in i en blokk som lar oss importere den til en annen fil (som heter en "klasse" for de spesielt interesserte). Dette vil la oss prøve ullike parametere for å se hva som funker for å skyte opp raketten når vi vil se nærmere på hva vi trenger for å nå stjernene.
Deretter kan vi også passe på at simulasjonen vår faktisk gjør det den skal, før vi prøver å bygge noe med resultater som kanskje har feil. Fra forrige del, vet vi at gjennomsnitt hastighet i hver komponent av alle partikellene skal være ca. 0. Vi vet også at gjennomsnittlig v² skal være;
\(\langle V² \rangle = \frac{K \times T}{m}\)
Der K er Boltzmann konstant, T er temperatur, og m er massen.
Vi kan få koden vår til å skrive ut gjennomsnittsverdiene til V og V², og sammenligne disse med verdiene vi forventer å se ifølge teorien; disse burde være lignende.
Vi kjører koden med 100,000 partikler, 3500 K, lengde på boksen som 10⁻⁶, massen til hydrogen atomer, tidssteg på 10⁻¹², og totaltid på 10⁻⁹. Når vi kjører den, ser vi vanligvis gjenomsnittshastigheter mellom -15 og 15. Det er hvert å nevne at standardavviket vi ville forventa rundt 0 her gitt ved:
\(\sigma = \sqrt{\frac{3000 \times 1.38 \times 10^{-23}}{1.66 \times 10^{-27}}} \approx 5000 \)
Altså ligger vi veldig godt innenfor hva vi ville forventa å få. Vi kan også sjekke at V² er der vi ville forventa. Her ser vi at vi ender opp med mellom 0.2% og 0.3% avvik fra forventingsverdien de fleste ganger vi kjører den.
Sammenligningen mellom forventingsverdiene til V og V² og det vi får indikerer at simulasjonen vår funker! Nå kan vi kaste dette inn i et annet program for å sjekke hva slags parametere vi trenger for å få oss vekk fra denne steinen.