Hvor skal vi starte?
Det f?rste vi kan begynne med er ideen bak boksen. Hvordan skal den funke? Hvordan kan vi simulere dette? Vil det funke? Dette er alle sp?rsm?l vi skal f? svart p? nedenfor, og vi kommer til ? bruke Python for ? gjennomf?re det.
Ideen bak hvordan boksen skal funke er at det er masse partikler i boksen som fyker rundt, og som gir fra seg bevegelsesmengde til raketten. S? lenge partikkelen er i boksen og kolliderer med veggene, vil bevegelsesmengden pr. vegg hoppe litt frem og tilbake, men v?re omtrentlig 0 fordi bevegelsesmengden p? den ene veggen "kanselleres" n?r partikkelen kolliderer med veggen p? andre siden. MEN med en gang en partikkel unnslipper ut gjennom hullet p? bunnen av boksen, vil ikke bevegelsesmengden fra partikkelen p? taket bli "kansellert" av kollisjonen partikkelen har med gulvet i bunnen, og raketten "beholder" bevegelsesmengden partikkelen overf?rte til taket av boksen. For de av dere som har fysikken p? plass, s? husker dere at \(\vec{F} = \frac{\Delta \vec{p}}{\Delta t}\) der \(\Delta\vec{p}\) er endringen i bevegelsesmengde og \(\Delta t\) er endring i tid, og dermed vil det virke en kraft p? boksen. Denne kraften trenger vi for ? vite hvor mye kraft raketten f?r fra motoren til slutt, og n?r vi begynner ? telle antall partikler som forlater boksen kan vi regne ut massen som forsvinner ved ? gange opp antall partikler ut med massen til partiklene.
S? hvordan lager vi denne boksen? Vi starter med ? definere alle de n?dvendige parameterne som temperatur, antall partikler og st?rrelse p? boksen. Det er viktig ? f? med at antallet partikler i boksen vil v?re det samme hele tiden, dvs. at n?r en partikkel forsvinner ut, vil en ny partikkel komme inn igjen. Videre vil vi bruke numerisk integrasjon for ? finne hastighet og posisjon til hver partikkel i boksen. Metoden for numerisk integrasjon vi vil bruke her kalles Eulers metode for numerisk integrasjon. Det er ikke den beste eller mest avanserte metoden for integrasjon, men vil v?re godt nok for et kort tidsintervall (i senere bloggposter kommer vi til ? bruke mer n?yaktige og avanserte metoder for numerisk integrasjon). P? denne m?ten vet finner vi hastigheten og posisjonen til partiklene. Mens partiklene virrer rundt m? vi ogs? teste for kollisjoner mellom partikkel og vegg, og snu fortegnet p? en av komponentene i hastigheten s? partikkelen ikke flyr vekk gjennom veggene. Vi blir ogs? n?dt til ? registrere n?r en partikkel forsvinner ut av hullet og antallet partikler som forsvinner ut, slik at vi kan regne ut hvor mye masse vi taper og kraften vi f?r i boksen. N? tar vi hvordan-delen i litt dypere detalj:
Vi setter noen faste parametere n?, som vi kan tweake p? senere. I f?rste omgang bruker vi:
\(\text{sidelengde p? vegg}=10^{-6}m,\:T=3000K,\:\text{antall partikler} = 10^5,\:\text{sidelengde p? hull}=0m\)
Som du ser starter vi uten noe hull i bunnen av boksen. Den legger vi til litt senere. Det er ogs? viktig ? presisere hvordan vi har definert koordinatsystemet i boksen. Boksen er definert slik at origo (punktet (0,0,0)) ligger midt i boksen. Det betyr at dersom sidelengden til boksen er \(L=10^{-6}m\), vil f.eks veggen langs x-aksen i boksen ligge p? koordinatene \(\left(\frac{L}{2},0,0\right)\) og \(\left(-\frac{L}{2},0,0\right)\), og analogt med y og z-aksen.
Det neste steget er ? lage to arrays fylt med null: en som skal inneholde posisjonene til alle partiklene ved et tidspunkt, og en som skal inneholde alle hastighetskomponentene til partiklene. En array er en slags datastruktur, som en kommode med skuffer og inndelinger oppi skuffene s? vi kan samle dataene v?re der, litt som i figur 6. Eneste forskjellen er at vi kun vil ha det for ett tidspunkt om gangen, alts? vil vi ha 100 000 antall skuffer for partikler, 1 inndeling for tid, og 3 inndelinger for hver komponent. Dette er fordi om vi skulle lagret dataen for alle posisjonene og hastighetene til 100000 partiklene for alle de sm? tidsstegene (vi kommer til ? bruke et lite tidssteg \(dt = 10^{-12}\) og g? gjennom 1000 steg med integrasjon) vi tar, ville vi endt opp med 600 000 000 elementer som skulle lagres! Da vil du rett og slett ende opp med et minneproblem p? pc'en din som gj?r at du ikke f?r kj?rt koden at all.
Euler - en ganske allsidig mann
N? som vi har disse arrayene fulle av null, ?nsker vi ? sette en initialposisjon og initialhastighet for hver enkelt partikkel. Hvorfor vil du forh?pentligvis forst? senere, om du allerede ikke vet det. Det finnes mange gode metoder for ? gj?re dette, men vi valgte en enkel og kanskje litt treg m?te ? gj?re det p?: nemlig g? gjennom arrayene og sette dem for hver partikkel én etter én. Men hva skal starthastighetene og startposisjonene v?re? N? f?r vi bruk for statistikken fra tidligere! Vi vet nemlig at en ideell gass har gaussisk distribusjon p? hastighetskomponentene. Vi gir hver enkel partikkel et tilfeldig tall fra gaussisk distribusjon, og vi brukte en random-funksjon i Python som gj?r dette for oss, kalt random.gaussian(mu, sigma), som gir oss et tilfeldig tall fra gaussisk distribusjon. Da sender vi inn \(\mu = 0\) og \(\sigma = \sqrt{\frac{kT}{m}}\) som gjelder for ideell gass. Hva med posisjon? Vel, partiklene har ingen god grunn til ? befinne seg p? et sted i rommet fremfor et annet sted. Derfor fordeler vi partiklene uniformt i rommet, som vil si at det er like sannsynlig ? finne en partikkel et sted i boksen som et annet sted i boksen. Vi bruker samme metode som for hastigheten, og bruker en random.uniform(-L/2, L/2) til ? fordele posisjonene. Da har vi alt vi trenger for ? integrere numerisk med Eulers metode.
Eulers metode er en form for numerisk integrasjon. Det man gj?r er ? regne ut hastighet og posisjon et biiiiiittelite tidssteg \(dt\) etter den forrige hastigheten og posisjonen, gitt at man har et uttrykk for akselerasjonen. Det finnes mange m?ter ? vise hvor Eulers metode kommer fra, men jeg liker ofte ? vise det med definisjonen av den deriverte:
\(\begin{align*} a(t) &= \frac{v(t+dt) - v(t)}{dt} \\ \Rightarrow v(t + dt) &= v(t) + a(t)\cdot dt \\ \\ v(t) &= \frac{r(t+dt) - r(t)}{dt} \\ \Rightarrow r(t+dt) &= r(t) + v(t)\cdot dt \end{align*}\)
Her faller uttrykkene rett ut fra definisjonene! Farten i et lite tidssteg dt etter tiden t er farten i tidspunktet t pluss akselerasjonen i tidspunktet t ganget med det lille tidssteget, og samme med posisjonen. Dette er Eulers metode! Husk at \(dt\) er en infinitesimal st?rrelse, alts? en uendelig liten st?rrelse, s? v?re tall p? datamaskinen blir kun en approksimasjon, og vi m? bruke sm? tidssteg \(\Delta t \). I v?rt tilfelle akselererer ikke partiklene, alts? er \(a(t) = 0\), s? farten slipper vi ? regne ut (husk at bevegelsesmengden er bevart, s? farten til partiklene er de samme hele veien, men bytter fortegn p? en komponent ved kollisjon med vegg. Glemt det allerede? Sjekk det ut her!). Uttrykket vi sitter igjen med ? datamaskinen er da kun \(r(t + \Delta t) = r(t) + v(t)\cdot\Delta t\) som vi trenger ? regne ut. Det vi gj?r da er ? bruke en loop til regne ut tidssteg for tidssteg. Vi bruker et biiitelite intervall \(\Delta t = 10^{-12}\), og regner ut 1000 steg av loopen. Da vil den totale tiden for simulasjonen av boksen v?re \(t_{total} = 10^{-9}\). Sekvensen vil se noe s?nn ut:
\(\begin{align*} r(10^{-12}) &= r(0) + v(0)\cdot10^{-12} \\ r(2\cdot10^{-12}) &= r(10^{-12}) + v(10^{-12})\cdot10^{-12} \\ r(3\cdot10^{-12}) &= r(2\cdot10^{-12}) + v(2\cdot10^{-12})\cdot10^{-12} \\ &... \\ &... \\ r(999\cdot10^{-12}) &= r(998\cdot10^{-12}) + v(998\cdot10^{-12})\cdot10^{-12} \\ r(1000\cdot10^{-12}) &= r(999\cdot10^{-12}) + v(999\cdot10^{-12})\cdot10^{-12} \\ \end{align*}\)
Ogs? er den ferdig etter 1000 steg. Ser du n? hvorfor vi trengte en initialhastighet og initialposisjon ved t=0? Vi bruker vektorisering for ? regne ut posisjonen til alle partiklene samtidig, s? det skal g? litt fortere.
Mellom hver eneste gang vi regner ut en ny posisjon for partiklene, looper vi gjennom alle partiklene for ? sjekke posisjonene deres og evt. bytte fortegnet p? en komponent dersom de har kollidert med en vegg. Vi kunne ogs? brukt vektorisering for ? sjekke veggkollisjoner, som ville v?rt kjappere og mer elegant, men ? g? gjennom alle partiklene en etter en funker det ogs?. M?ten vi tester det p? er ? teste om hvert enkelt koordinat for hver partikkel er st?rre eller mindre enn \(\frac{L}{2}\) eller \(-\frac{L}{2}\). P? den m?ten har vi kontroll p? om koordinatene til partikkelen er i boksen eller ikke. Da m? vi passe ekstra godt p? om \(\Delta t \)-en v?r er liten nok, ellers vil partikkelen fly ut av veggen p? boksen mellom et tidssteg uten at vi fanger det opp. Vi kj?rer simuleringen for \(N=3\:\text{antall partikler}\) og plotter banene til partiklene, bare for ? se om det ser riktig ut:
Dette ser veldig lovende ut! Partiklene spretter fint inne i boksen uten ? forsvinne ut. En annen m?te ? sjekke om det stemmer p? er ? regne ut det analytiske trykket i boksen og det numeriske trykket i boksen, og sammenlikne. Det analytiske uttrykket for trykk av gass i et volum er \(P = \frac{N}{V}kT\) , der \(N\) er antall partikler, \(V \) er volumet gassen befinner seg i, \(k\) er som tidligere Boltzmanns konstant og \(T\) er temperatur i kelvin. Vi kj?rer en simulering av boksen der vi regner ut det numeriske trykket. Vi husker at trykk ogs? er gitt som \(P = \frac{F}{A}\), men vi mangler et uttrykk for kraften. Heldigvis utledet Johan dette i en tidligere post om bevegelsesmengde som du kan lese her. Det tar vi for gitt, og vi vet n? at kraften fra en enkel partikkel er \(F = \frac{2p_z}{\Delta t}\), der \(p_z\) er bevegelsesmengde for partikkelen i z-retning, og \(\Delta t\) er tiden for hele simuleringen (\(10^{-9}s\)). Dette gir oss uttrykket \(P = \frac{2p_z}{\Delta tA}\), der alle nevnte st?rrelser er som for, og \(A\) er arealet p? flaten kraften virker p?. Vi summerer dette sammen for hver gang vi registrerer at en partikkel traff en eller annen vegg, og deler det p? antall sider i boksen for ? f? et fint gjennomsnitt av trykk pr. vegg. Det numeriske trykket blir da \(P=4111.093Pa\). Vi kaster de samme tallene inn i formelen for det analytiske trykket, og f?r \(P=4141.947Pa\). Litt avvik m? vi forvente fra datamaskinen og dataen vi har, og tallene er veldig like. Dermed kaller jeg det suksess! Vi kan n? introdusere hullet i bunnen av esken.
N?r vi introduserer hullet i esken er det tre nye ting ? ta hensyn til: st?rrelsen p? sidelengden av hullet, testing av n?r en partikkel g?r ut av hullet, og hva som skjer med en partikkel n?r den forlater hullet (HUSK! Antallet partikler i boksen skal v?re det samme under hele simuleringen). Vi legger inn at sidelengden p? hullet er halve lengden av boksen, dvs. \(\text{sidelengde hull} = 0.5\cdot10^{-6}\), slik at arealet av hullet er en fjerdedel av arealet p? bunnen. N?r vi skal registrere om en partikkel g?r ut av hullet, bruker vi samme praksis som da vi testet for veggkollisjon. Under testen der vi sjekket om partikkelen treffer gulvet legger vi en ny test, som tester om x og y-koordinatet til partikkelen er mindre enn sidelengden til hullet. Hvis x og y-koordinatet til partikkelen er mindre enn sidelengden til hullet, registrerer vi at en partikkel gikk ut og legger til +1 p? antall partikler som har g?tt ut. Denne starter selvf?lgelig p? 0. I tillegg regner vi ut bevegelsesmengden \(p_z\) som overf?res til taket av boksen. Denne kan vi bruke til ? regne ut total kraft fra boksen senere, ved ? bruke at \(F = \frac{2p_z}{\Delta t}\), der \(\Delta t \) er tiden for hele simulasjonen. Siden \(\frac{2}{\Delta t}\) er felles for kraften fra alle partiklene trenger vi kun ? summere opp bevegelsesmengden, og kaste den inn i uttrykket senere. Vi regner ut \(p_z\) ved ? summere opp bevegelsesmengden i z-retning til en partikkel hver gang den forlater hullet, selv om bevegelsesmengden teknisk sett overf?res n?r partikkelen treffer taket. Men p? denne m?ten slipper vi ? legge til og trekke fra bevegelsesmengde hver gang en partikkel kolliderer. Totalt massetap regner vi ogs? etter at loopingen er ferdig, siden\(\text{totalt massetap} = \text{antall partikler}\cdot \text{massen til partiklene}\). Men siden det kun er for \(10^{-9}s\), m? vi gange det opp med \(10^9\) (eller dele det p? den totale tiden for simuleringen) for ? f? totalt massetap pr. sekund. Da er \(\text{totalt massetap pr. sekund} = \frac{\text{antall partikler}\cdot\text{massen til partiklene}}{\text{total tid}}\). Men hva skjer n? med partikkelen som forlot boksen? Vi gj?r det veldig enkelt, og bare kaster den tilbake i posisjonen \(\left(0,0,\frac{L}{2} - 10^{-12}\right)\) som om ingenting har skjedd. Vi trekker fra en bitteliten del av posisjonen slik at koden ikke skal registrere det som en veggkollisjon og snu fortegnet p? z-komponenten av hastigheten til partikkelen. Vi velger derimot ? beholde hastigheten til partikkelen som forsvant ut, for ? bevare bevegelsesmengden og den gaussiske distribusjonen. Hadde vi valgt ? gi den "nye" partikkelen som kommer inn en nye gaussiske hastighetskomponenter kunne det hende at vi forskyver den gamle distribusjonen, og ?delegger den gamle distribusjonen.
Resultater fra boksen
Vi kj?rer en simulering av boksen med parametere:
\(L=10^{-6}m,\:T=3000K,\:N=10^5,\:l=0.5\cdot10^{-6}m\) , der \(L\) er sidelengden p? boksen, \(T\) er temperaturen i boksen gitt i kelvin, \(N\) er antall partikler i boksen, og \(l\) er sidelengden p? hullet p? bunnen av boksen. Tallene for kraft og massetap pr. sekund pr.boks blir da:
| N (antall partikler) | \(45934999999999 \) |
| thrust (N : Newton) | \(7.329079237148241\cdot10^{-10}N\) |
| Bevegelsesmengde (kgm/s) | \(7.329079237148241\cdot10^{-10}kgm/s\) |
| massetap pr. sekund (kg/s) | \(1.5374972752499938\cdot10^{-13}\frac{kg}{s}\) |
Det som kanskje er verdt ? kommentere p? f?rst her er at total bevegelsesmengde fra boksen ila. et sekund er det samme som thrust-force for raketten. Det h?res kanskje veldig rart ut, men dersom vi har total mengde bevegelsesmengde ila. simuleringstiden p? \(10^{-9}s\) s? vil det v?re ca. \(p = 7.33\cdot10^{-19}kgm/s\). Dersom vi deler dette igjen p? \(10^{-9}s\) vil vi f? total bevegelsesmengde pr. sekund, og vi ser p? dimensjonene ogs? at det g?r fra \(\frac{kgm}{s}\) til \(\frac{kgm}{s^2}\), som er enheten for kraft. Vi vet at arealet p? undersiden av raketten v?r er \(16m^2\), s? dersom én liten boks har sidelengde \(10^{-6}m\) har boksen en underside med areal \(10^{-12}m^2\). Det betyr ogs? da at \(10^{12}\) antall bokser vil ta opp \(1m^2\). Dermed kan vi tillate oss ? bruke s? mye som \(1.6\cdot10^{13}\) antall bokser uten at det tar for mye plass. Bruker vi s? mange bokser for ? lage en rakket vil vi ende opp med en motor som har drivkraft p? ca. \(11.72\cdot10^{3}N\) og massetap \(2.46\frac{kg}{s}\).
Dette er tall med parametere vi kommer til ? endre senere. For eksempel kan det hende vi gj?r boksen mindre slik at vi f?r h?yere antall-tetthet av partikler, som gir oss mer trykk i boksen og derav mer kraft. Med justering av boks-st?rrelse kommer ogs? naturligvis en justering av utslippshullet. Vi kommer til ? justere antall partikler, og antakeligvis putte litt flere partikler inn i boksen. Parametere vi skal pr?ve ? holde ur?rt temperatur, mest fordi forbrenningstemperatur p? hydrogengass faktisk er \(3000K\). H?yere temperatur vil ogs? gi oss mer kraft, men vi pr?ver ? holde den der den er. Vi kommer ogs? til ? beholde et antall bokser slik at rakettmotoren vil ta opp alle de \(16m^2\) vi har til gode p? undersiden av raketten v?r, dermed justeres antallet ut ifra st?rrelsen p? boksen.
Skal vi til verdensrommet?
For ? f? k?l p? dette m? du til neste innlegg.