Raketfysik

Om opgaverne: På denne side finder du en række opgaver som omhandler opsendelsen af raktter fra Jordens overflade, samt navigation i solsystemet. Målet med opgaverne er primært træning af programmering i MATLAB og sekundært fysisk forståelse ved at se samspillet mellem flere effekter der alle påvirker resultatet. Opgaverne et stillet i forbindelse med Videregående Mekanik (Mek2) på Niels Bohr Institutet for årgang 2011.

Målet med en raketopsendelse er at aflevere en nyttelast i den ønskede position på det rette tidspunkt. Satelitter og rumsonders position er påvirket af tyngdefeltet fra Solen og planeterne. En satellit er derfor altid "i kredsløb" om et eller andet. Når man navigerer i solsystemet er det derfor et spørgsmål om at gå fra én bane til en anden på det rette tidspunkt, er tidspunktet forkert ender man måske milioner af kilometer fra sin tiltænkte planet, da den har flyttet sig i mellemtiden.

Starter med en opsendelse fra Jorden, er det første problem at opnå så høj en hastighed at nyttelasten falder "forbi" Jorden og derfor befinder sig i et kredsløb. Et lavt kredsløb få hundrede kilometer over jordoverfladen kræver en hastighed omkring 7.3 km/s. For cirkulære kredsløb er hastigheden givet ved [Y&F, 13.10],

\begin{equation} v_{circ} = \sqrt{ \frac{G m_E}{r} } \end{equation}

1) Simulering af raketligningen

Kraftpåvirkningen på en raket (kaldet "Thrust") er givet ved, \begin{equation} F = - V_{ex} \frac{d M}{d t} = - V_{ex} \ \dot{m} \end{equation} Her er \( V_{ex} \) udstødningshastigheden for reaktionsmassen, og \( \dot{m} = dM/dt \) er den instantane masseændring.

Opgaver (Download MATLAB skabelon)

For en raket kender vi \( V_{ex} = 2500\) m/s og raketten vejer \(1.2 \times 10^6\) kg, hvoraf de \(1.0\times 10^6 \)kg er brændstof, resten er nyttelast. Raketten har et trin som brænder i 205 sekunder.

  1. Bestem \( \dot{m} \) hvis vi antager at forbrænding er konstant over brændetiden.
  2. Benyt det analytiske udtryk \( \Delta V = - V_{ex} \ln \frac{m_0}{m} \) hvor \(m_0\) er den oprindelige masse (nyttelast og brændstof) og \(m\) er den resterende masse efter brændstoffet er brugt, til at bestemme den forventede hastighedsændring.
  3. Plot den analytiske forudsigelse for \(\Delta V\) som funktion af tiden. (Hint: Antag at alt brændstoffet er brugt efter 205 sekunder, og lav en vektor i MATLAB med den den resterende masse som funktion af tid).
  4. Skriv en numerisk simulering i én dimension, der beregner hastigheden og positionen af raketten i løbet af de \(205\) sekunder ud fra thrust-ligningen. Benyt Eulers metode.
  5. Plot \( \Delta V \) som funktion af \(t\) for både det analytiske og numeriske resultat. Indsæt en "legend" samt labels.
  6. Stemmer den numeriske løsning overens med den analytiske? Se hvor stor en påvirkning valget af tidskridt har, ved at ændre \(\Delta t \).
  7. Hvor langt har raketten bevæget sig efter \(205\) s?

2) Multi-trin raketter

Normalt benytter man mere end et rakettrin for at nå kredsløb. Fordelen ved opdelingen er at man kan smide den tomme raket væk og derved spare må massen. Ved udregningerne i opgave 1, er forskellen på den tomme og den fulde raket brændstoffet. Ved flere trin, må vi holde øje med hvad der er nyttelast, hvad der er raketstruktur som kan afkobles når et trin er brændt ud, og hvad der er brændstof i de enkelte trin.
En 2-trins raket består af følgende massebidrag før start: \begin{equation} m_{total} = m^{trin \ 1}_{struktur}+ m^{trin \ 1}_{brændstof} + m^{trin \ 2}_{struktur} + m^{trin \ 2}_{brændstof} + m^{nyttelast} \quad \text{"Start"} \end{equation} Under afbrændingen af det første trin går \( m^{trin \ 1}_{brændstof} \rightarrow 0\) og de resterende led forbliver konstante. Efter første trin er brugt frakobler man den tomme struktur, og den totale masse ved tændingen af trin 2 ser således ud, \begin{equation} m_{total} = m^{trin \ 2}_{struktur} + m^{trin \ 2}_{brændstof} + m^{nyttelast} \quad \text{"Trin 1 afkoblet"} \end{equation} Her er det nu \( m^{trin \ 2}_{brændstof} \) som ændrer sig, indtil de sidste trin er brugt, og strukturen frakobles, og vi står tilbage med nyttelasten som den totale masse, placeret i den tiltænkte bane. \begin{equation} m_{total} = m^{nyttelast} \quad \text{"Trin 2 afkoblet"} \end{equation}

Det amerikanske firma SpaceX har udviklet en raket kaldet Falcon 9, som skal benyttes til at forsyne den internationale rumstation fra starten af næste år. I resten af opgaverne bruger vi specifikationerne fra Falcon 9 som reference for at give realistiske værdier. Falcon 9 er en 2-trins raket med følgende specifikationer:

Trin 1Trin 2Nyttelast
\( V_{ex}\)2837 m/s3292 m/s
burntime170 s345 s
\( \dot{m}\)-1972.94 kg/s-141.74 kg/s
\( m_{struktur}\)24700 kg3100 kg
\( m_{brændstof}\)335400 kg48900 kg
\( m_{nyttelast}\)10454 kg
Mængden af nyttelast kan justeres, jo mindre raketten medbringer jo højere fast kan den opnå, og det bestemmer i sidste ende hvilket kredsløb man kan opnå.

Opgaver (Download MATLAB skabelon)

  1. Udvid simuleringen fra opgave 1 med et ekstra trin, som aktiveres efter det første er brændt ud (Hint: du kan vælge flere fremgangsmåder, den nemmeste er at kopiere løkken fra trin 1 og sætte den efter med opdaterede værdier for trin 2).
  2. Bestem den totale hastighedsændring analytisk for ovenstående parametre.
  3. Plot position som funktion af tid
  4. Plot hastighed som funktion af position
  5. Plot hastighed som funktion af tid
  6. Ovenstående simulering tager ikke højde for tyngdekraft, luftmodstand o.lign. vil det give nogen fordel at indføre en "coasting" fase mellem de to trins afbrændinger hvor trin 2 forbliver slukket i en periode?

3) Tyngdekraft

I Tyngdekraften påvirker opsendelsen på to måder. Den ene er som nævnt, at et objekt kræver en specifik hastighed for at være i kredsløb. Den anden påvirkning kaldes for "gravity drag" (tyngdemodstand), og svarer til den kraft der skal modvirkes før at raketten ikke længere er udsat for 1g tyngdeacceleration.
Hvis vi antager at Jorden er flad vil tyngdemodstanden være et ekstra led i hastighedsopdateringen:

            v = v + dt * (F_thrust/m_raket + g)
        

Da tyngdeaccelerationen ikke ændrer sig nævneværdigt de første 100 km over jordoverfladen er det en acceptabel approksimation.
Da målet er at komme i et regulært kredsløb om Jorden, er det dog smartere at skrifte til vektorkoordinater, og bestemme Jordens tyngdekraft på fartøjet,

            F_g = -G .* m_raket .* m_jorden ./(norm(r).^3) .* r; % (N) Kraft
            v = v + dt * (F_thrust/m_raket + F_g/m_raket)
        

Hvor G er Newtons konstant, og r er en vektor med jordens centrum i origo.

Opgaver: (benyt skabelon koden fra del 2. men udvid med vektorer for r og v)

  1. Udvid simuleringen fra del 2 med effekten af tyngdekraft.
  2. (ekstra) I opgave 8.112 i [F&Y] fandt vi det analytiske udtryk for raketligningen med tyngdeacceleration, plot udtrykket sammen med den numeriske løsning. Stemmer de overens?
  3. Plot position som funktion af tid
  4. Plot hastighed som funktion af position
  5. Plot hastighed som funktion af tid

4) Atmosfæren

Hvis Jorden var atmosfæreløs (og helt flad) ville alt der kræves for at opnå kredsløb være en tilpas høj tangentiel hastighed. Et stort problem ved raketopsendelser er at opnå de høje kredsløbshastigheder uden at luftmodstanden får raketten til at smelte. En typisk raketopsendelse starter ved at det første trin bringer raketten op i omkring 50 km højde, her vil hastigheden være omkring 6000 km/t. 50 km oppe er stratosfæren ophørt og luftens massefylde er over 250 gange lavere end ved havoverfladen. Med den tyndere atmosfære er raketten udsat for langt mindre luftmodstand og hastigheden kan forøges betragteligt.

I samme højde udføres der også et "gravity turn" som bare betyder at man drejer raketten så den flyver gradvist mere horisontalt indtil den flyver tangentielt med Jordens overflade.

En simpel model for en atmosfæres massefylde er givet ved,

            function a = rho( h )
            % Adiabatic atmosfære
            % http://farside.ph.utexas.edu/teaching/sm1/lectures/node56.html
            rho0 = 1.225; % kg/m^4-1 - massefylden i havniveau
            gamma = 1.4; % 
            h0 = 29400; %km

            a = rho0 * (1 - (gamma-1)./gamma .* h./h0).^(1/(gamma-1));
            a(h > 29400) = 0;
            end    
        

Hvor \(h\) er højden i meter over havoverfladen. Modellen er ikke realistisk i højder over ~30 km, hvor massefylden er eksponentielt aftagende. Derfor retunerer funktionen 0 for h > 29400 m.

Luftmodstand

Vi har allerede i Mekanik 1 set på luftmodstanden, som approksimativt kan beskrives ved, \begin{equation} F_{luft} = \frac{1}{2} \rho(z) v^2 C_d A \end{equation} Hvor \(F_{luft}\) er luftmodstanden i flyveretningen, \(\rho(z)\) er massefylden som funktion af højden over havoverfladen, v er farten i flyveretningen, \(C_d\) er en dimensionsløs dragkoefficient og \(A\) er rakettens tværsnit i flyveretningen.

            F_luft = @(v,A,rho,C_D) (C_D .* 0.5 .* rho .* A .* v.^2);
        

Hvis vi antager at raketten altid er vendt i flyveretningen så tværsnittet \(A\) er konstant, og farten svarer til \( |\vec{v}|\), kan vi skrive et vektoriseret udtryk ved at multiplicere med hastighedens enhedsvektor,

            F_drag = - v/norm(v) .* F_luft(norm(v) , A, rho(norm(r) - R), C_D);
        

Her antager vi at positionsvektoren r har formen \( \langle x,y,z \rangle \) m og R er Jordens radius i meter.

Opgaver

  1. Udvid simuleringen fra del 3 med en kraftpåvirkning fra luftmodstand
  2. Plot position som funktion af tid
  3. Plot hastighed som funktion af position
  4. Plot hastighed som funktion af tid
  5. (Ekstra) Ved at dreje hastighedsvektoren en smule per tidsskridt kan rakettens retning ændres. Indfør en drejning således at raketten flyver tangentielt på Jorden efter den har nået 50 km.
    Den anonyme funktion rotxy(theta, vec) roterer en 3D vektor i xy planet med en vinkel theta,
    	    rotxy = @(theta, vec) [cosd(theta), -sind(theta), 0; ...
                                    sind(theta),  cosd(theta), 0; ...
                                    0,            0,           1] * vec';
            
    Benyt denne funktion, på hastighedsvektoren når betingelsen norm(r) > 50000 er opfyldt. (Hint: et ægte gravity turn udnytter tyngdekraften til at trække rakettens retningsvektor mod Jorden, med små ændringer i vinklen, undersøg hvordan!)
  6. (Ekstra) Hvis man affyrer raketten mod øst kan man udnytte Jordens egenrotation til at få op mod 460 m/s yderligere omløbshastighed. Det faktisk bidrag afhænger af afstanden til ækvator, med følgende udtryk, \begin{equation} \Delta v_{rotation} = \Omega \ R_{jord} cosd(L) \ cosd(\beta) \quad \text{hvor } \quad\Omega = 7.26 \times 10^{-5} rad/s \end{equation} \(L\) er bredegraden og \(\beta\) er vinklen fra fra vandret affyring mod øst.
    Vi antager at afsendelsen foregår fra Cape Caneveral, hvor \(L=28 deg\). Plot \(v_{rotation}\) som funktion af \(\beta\).
  7. (Ekstra) Indfør at raketten flyver mod øst i simulationen.

5) Kredsløb om Jorden

Hvis raketten har haft tilstrækkelig energi til at levere raketten i en højde hvor banehastighden matcher rakettens vil vi nu være i et kredsløb.
Selve simuleringen af rakettens bane er allerede indbygget i del 3, da vi bruger det fulde vektoriserede udtryk for tyngdekraften. Derfor kan vi forsætte simuleringen som hidtil, men uden effekten fra thrust og luftmodstand, og vi er tilbage med effekten af tyngdekraften alene,

        while t < t_stop
            t = t + dt;
            
            F_g = -G .* m_raket .* m_jorden ./(norm(r).^3) .* r; % (N) Kraft
            r = r + v * dt;
            v = v + dt * (F_g/m_raket);
        end
        

Typisk er banen efter opsendelsen stærkt elliptisk, og for at undgå at ramme Jorden på vej mod perigeum laver man flere korrektioner efter man har opnået størstedelen af kredsløbshastigheden. Det klassiske eksempel på denne type baneændring kaldes for et Hohmann transfer.

Opgaver:

  1. Forsæt simuleringen efter begge rakettrin er udbrændt og afkoblet (Hint tilføj evt. endnu en løkke som den vist ovenfor, med m_raket, t, r og v startende med værdierne fra slutningen af trin 2 simuleringen).
  2. (Avanceret) Lad os antage at rumfartøjet allerede er i kredsløb med startbetingelserne: \(\vec{r} = \langle 0, R_{jord} + 400000, 0 \rangle \) m, \(\vec{v} = \langle v_{circ}(R_{jord} + 400000), 0, 0 \rangle \) m/s. Benyt ligningerne på Wikipedia og i kildekoden på kredssløbssimulatoren til at beregne hvor meget \(\Delta V\) en baneændring fra 4000 km højde til 300.000 km vil kræve.
  3. (Ekstra, bonus, hardcore, brug ikke din læseferie på det her...) Indfør et kommandomodul der alt efter bane højde, tid og/eller position laver korrektioner i \(\Delta V \) med Hohmann transfers (Hint: se kildekoden fra kredssløbssimulatoren her). Hvis vi antager at nyttelasten har raketter til banekorrektion og et begrænset \(\Delta V\) budget, benyt en serie transfers til at få nyttelasten i et cirkulært kredsløb om Jorden, efter opsendelsen.

6) Solsystemet (førsteårsprojekter)

Denne del er mest af alt en perspektivering, og bør ikke påbegyndes uden godt med fritid, eller som et decideret førsteårsprojekt.

For et rumfartøj i kredsløb om Jorden taler man om et "reduceret 2-legeme problem", da fartøjets påvirkning af Jordens bane er ubetydelig, det er kun Jordens tyngdepåvirkning som påvirker fartøjets retning. For Jord-Månesystemet eller Solsystemet som helhed er det anderledes kompliceret. En mission til Jupiter eller Merkur kræver præcis beregning af planeternes indbyrdes position. Det nytter ikke noget at man kan lave den rigtige Hohmann manøvre der vil bringe en rumsonde ud til Jupiter, hvis planeten befinder sig på den anden side af solen når rumsonden ankommer!

I kildekoden til del 6 findes startbetingelserne for hele solsystemet, og en Euler simulation af planeternes bevægelse, som planeterne stod i july 2008. Samlet svarer det til et 9-legeme problem.

Opgaver (Download MATLAB skabelon)
  1. Kør koden og se hvad der sker ved forskellige dt antagelser. N-legeme problemer som dette er meget følsome overfor startbetingelser og den numeriske præcision. Bare små unøjagtigheder kan akkumuleres og simuleringen kan enten gå i stå, eller legemerne kan styrte mod hinanden for så igen at skyde ud og aldrig vende tilbage!
  2. Opgrader den numeriske integrator til Runge-Kutta. En måde at redde simuleringen fra katastrofale effekter er ved at benytte små dt værdier. Desværre vil simuleringen gå smertefuldt langsomt. I stedet kan man benytte en forbedre integrationsmetode som kan give bedre forudsigelser end Eulers metode, med større tidsskridt. Metoden kaldet "Runge-Kutta O(4)" er beskrevet i sektionen om differentialligninger og er en af de mest benyttede teknikker i videnskaben.
  3. Navigation mellem planeterne. Kombiner del 5 med del 6 ved at indføre et rumfartøj hvis position afhænger af kraftpåvirkningen fra alle planeterne (og solen). Fartøjet behøver ikke at påvirke solsystemet (Denne konfiguration kaldes et 9+1 problem). Start med samme begyndelsesbetingelse som i opgave 5.2, og design en mission til en anden planet.
    • Hvor mange Hohmann transfers skal der til?
    • Hvor stort et \(\Delta v\) budget skal der til? - Hvor meget brændstof?
    • Hvordan bestemmer man hvornår en manøvre skal udføres?
    • Flere problematikker
  4. En verden uden Neptun. Lav en simulering med begyndelsesbetingelserne fra raket_del7.m men hvor Neptun er taget ud. Lav den tilsvarende simulering med Neptun. Er det muligt at forudsige positionen af Neptun ud fra perturbationerne af Uranus?
  5. N-legeme simulering generelt. N-legeme simulering er et aktivt forskningsområde, specielt galaksedannelse og andre astrofysiske fænomener så som simulering af mørkt stof bygger på samme type simulering. Prøv at eksperimentere med N-legeme koden ved at lave en samling "partikler" med forskellig masse, position og hastigheder som startbetingelse og se hvordan de udvikler sig ved deres indbyrdes interaktion. (Advarsel: koden i linket er baseret på Euler metoden, benyt Runge-Kutta for at opnå bedre ydeevne ved \(N_{partikler} > 2\)).