Eötvös Loránd Tudományegyetem Természettudományi Kar Komplex Rendszerek Fizikája Tanszék
A hintamanőver szimulációja Szakdolgozat Készítette:
Puszta Adrián Fizika BSc
Témavezető:
dr. Csabai István ELTE TTK Komplex Rendszerek Fizikája Tanszék 2009
1
Tartalomjegyzék Bevezetés....................................................................................................................................3 A hintamanőver történeti áttekintése...........................................................................................4 A hintamanőver kinematikája.....................................................................................................6 A három- és többtest-probléma.................................................................................................10 Numerikus differenciálegyenlet-megoldó módszerek..............................................................12 A szimulációs módszer tesztelése.............................................................................................14 Valódi adatok alapján számolt pálya.........................................................................................20 Optimalizációs lehetőségek, újabb numerikus módszerek........................................................25 További fejlesztési lehetőségek.................................................................................................27 Technikai információk..............................................................................................................27 Köszönetnyilvánítás..................................................................................................................27 Irodalomjegyzék........................................................................................................................28 Függelék....................................................................................................................................30
2
Bevezetés Az égitestek megfigyelése mindig is nagy szerepet játszott az emberek életében. Az ókorban a bolygók mozgásához kötöttek vallási szertartásokat, a csillagképekhez mitológiai történeteket társítottak. Később, már a újkor hajnalán, a távcső feltalálásával lehetővé vált, hogy pontosabban is megfigyeljék a bolygók pályáit, és le tudják írni a bolygómozgás törvényeit (Kepler), vagy a köztük ható erőt (Newton). A 20. században viszont már nem elégedtek meg a Földről való megfigyeléssel, közelebb akartak kerülni az égitestekhez. Ezért űrszondákat építettek, hogy azokat a megfelelő helyre eljuttatva információt kaphassanak a bolygók légköréről, mágneses teréről, egyéb tulajdonságaikról. Az űrszondáknak a megfelelő helyre juttatása viszont nem egyszerű feladat. Ugyanis a gravitációs vonzások legyőzéséhez üzemanyagra van szükség. Viszont minél több üzemanyagot visz a szonda, annál nehezebb lesz, és a nagyobb tömeggel manőverezéshez pedig több üzemanyag kell. Ezért nagy segítség volt az űrhajózásnak a hintamanőver felfedezése. A hintamanőver célja, hogy egy bolygótól energiát nyerve megnövelje, vagy épp energiát átadva lecsökkentse az űrszonda sebességét. A manőver szemléltetésére egy vonatos hasonlatot szoktak említeni1. Eszerint egy vonat halad a megfigyelő felé A megfigyelő egy labdát dob a vonatnak
v vonat sebességgel.
v labda sebességgel. A vonatból nézve a labda
v vonat v labda sebességgel jön felé, majd ugyanekkora sebességgel visszapattan. A külső megfigyelő szerint a labda a visszapattanás után labda sebessége tehát
v vonat v labda v vonat sebességgel halad. A
2 v vonat -tal megnőtt. A hintamanővert is hasonlóan lehet elképzelni, a
vonat helyett valamelyik bolygót, a labda helyére a szondát kell helyettesíteni, valamint nem rugalmas ütközés, hanem gravitációs vonzás hat a két test között. A hintamanőver tehát egy olyan módszer az űrhajózásban, aminek segítségével minimális üzemanyag felhasználásával meg lehet növelni az űrszonda sebességét, és így meggyorsítani a célba érését. Tovább lehet növelni a hatékonyságát, ha a manőver közben a hajtóművek bekapcsolásával módosítják a pályát, így még nagyobb sebességre tehet szert a szonda.
1 http://www2.jpl.nasa.gov/basics/grav/primer.html
3
A hintamanőver történeti áttekintése Az
ukrán
származású
Yuri
Kondratyuk, az űrrepülés úttörője írta le először a hintamanővert (1918-ban), mint egy
lehetséges
módszert
űrhajók
gyorsítására. Munkáját csak orosz kollégái értékelték, a többiek viszont nem akartak ezzel a különös ötlettel foglalkozni. Nem sokkal később, 1920-ban Walter Hohmann német mérnök kiszámolta két bolygó közti legkisebb energiájú pályát. Ez a pálya egy ellipszis, ami érinti az indítási bolygó, és a célbolygó pályáját is. Ebben az időben ezt a
1. ábra: A Hohmann-ellipszis
pályát tekintették a leggyorsabb útnak a
célbolygóhoz, és az ezen való végighaladás idejét tekintették az utazás idejének. Még az 1950-es években is úgy gondolták, hogy nagyon hosszú időbe fog telni, mire a Jupiterig, vagy azon túl juttatnak ember készítette űrszondát. Kiszámolható, hogy a Plútóig eljutni ezen az ellipszisen 40-50 évbe kerülne, míg a Neptunuszig is legalább 30 évre lenne szükség. Ez tehát egy meglehetősen lassú módja a Naprendszer feltérképezésének. 1961-ben egy végzős matematikus, Michael A. Minovitch Kondratyuk után újra felfedezte azt a gondolatot, hogy ha megfelelő pályán indítják el a szondát, akkor extra sebességre tehet szert különösebb üzemanyagfelhasználás nélkül, más bolygók gravitációs terének felhasználásával. Addig ugyanis a gravitációs teret egy leküzdeni való akadályként kezelték, aminek leküzdéséhez csak az üzemanyag mennyiségét kell növelni. Minovitch ötlete az volt, hogy a szonda a bolygó mellett elhaladva energiát nyerhetne, miközben a bolygó energiája, és így sebessége elhanyagolható mértékben csökken. A módszer első valódi alkalmazása a Mariner 10 űrszonda volt 1973-ban. A misszió célja a Merkúr vizsgálata volt, annak fizikai adatainak mérése. Útja közben a Vénusz mellett hajtott végre hintamanővert, és így haladt tovább a Merkúr felé. A Merkúr mellett további három alkalommal „hintázott”, eközben felvételeket készített a bolygóról. Végül '75-ben kifogyott az üzemanyagból, és Nap-körüli pályára állt.
4
A Voyager-1 amerikai űrszonda 1977-ben indult. Ez a Jupitert, a Szaturnuszt, és ezek holdjait térképezte fel. A Jupitert használta fel a sebessége növelésére, majd újabb manőverek után elhagyta a Naprendszert. A legújabb adatok szerint (2009. február 1.) 16.247 milliárd km-re jár a Naptól, ezzel a legtávolabbi működő űrszonda lett. A hintamanőverek kihasználása nélkül nem lehetett volna ilyen távolságba eljuttatni. 1989-ben indult el a Galileo szonda, szintén a Jupiterhez. A Galileo a VEEGA pályát követte, a Venus-Earth-Earth Gravity Assist-ot, vagyis egyszer a Vénusz, kétszer a Föld mellett hajtott végre hintamanővert. A hintamanővert legtöbbet kihasználó szondák egyike a MESSENGER, amely egyszer a Föld, kétszer a Vénusz, majd háromszor a Merkúr mellett manőverezett, majd a Merkúr körüli pályára állt.
5
A hintamanőver kinematikája A hintamanővert (angolul gravity assist) fel lehet használni az űrszonda sebességének növelésére és csökkentésére is. Ezekre az esetekre felírhatjuk az impulzusokat, és azok megváltozásával megkaphatjuk a sebesség
változását.
Az
alábbiakban
ennek
mechanizmusát mutatom be ([2] alapján). A sebességeknek 3 indexet adtam: az első J vagy V, ami a Jupiter, ill. Voyager kezdőbetűi (így különböztetem meg a nagy tömegű bolygót a műholdtól); a második S vagy
C,
ami
a
Naphoz
rögzített
(S),
ill.
a
tömegközépponthoz (C) rögzített rendszert jelöli; a 2. ábra: A jelenség leírása harmadik pedig B vagy A, vagyis a kölcsönhatás előtti Naphoz rögzített (S frame) és (Before) ill. utáni (After) sebességvektort jelöli.
tömegközépponti (C frame)
Kezdeti feltételnek olyan helyzetet választottam, rendszerben [2] hogy a bolygó (továbbiakban Jupiter, J) az S rendszerben -x irányban halad, a szonda (továbbiakban Voyager, V) pedig y irányban a hozzájuk tartozó sebességekkel. Tömegközépponti rendszerben (C) a rendszer összimpulzusának nullának kell lennie. Ezért:
v JCB M J =v VCB M V , átrendezve:
behelyettesítve ez az arány:
v JCB M V = , ami a Jupiter és Voyager tömegeit v VCB M J
825 kg =4.38∗10−25 . Mivel ez egy nagyon kicsi szám, 27 1.88∗10 kg
ebből két dolog következik: 1. J sebessége elhanyagolható a TKP sebességéhez képest, így azzal megegyezik v JSB=v TKP ,S
2. a kölcsönhatás2 során nem változik J sebessége (mivel a TKP sebessége sem változik): v JSA=v JSB
2 Kölcsönhatásnak nevezem a műhold bolygó melletti elhaladását
6
Tömegközépponti rendszerben felírható V sebessége, ahogy az a 3. ábrán látható: v VSB − v JSB ⇒ ∣v VCB∣= v 2JSBv 2VSB v VCB = B=arctg
∣v VSB∣ = ∣v JSB∣ A
Ezek után kiszámíthatjuk V sebességét a kölcsönhatás után a Naphoz rögzített rendszerben:
∣v VCA∣ sin A ∣v VSA∣
v VSA=v JSA v VCA , sin =
3. ábra: v VCB számítása [2]
4. ábra: v VSA számítása [2]
A koszinusztétel segítségével ∣v VSA∣ meghatározható:
∣v VSA∣2=∣v JSA∣2∣v VCA∣2−2∣v JSA∣∣v VCA∣cos 180 ˚− A Az eddigi eredményeket behelyettesítve:
∣v VSA∣ =2∣v JSB∣ ∣v VSB∣ −∣v JSB∣∣v JSB∣ ∣v VSB∣ cos 180 ˚−arctg 2
2
Ha ezt az egyenletet leosztjuk
2
2
2
∣v VSB∣ ∣v JSB∣
∣v VSB∣2 -tel, akkor egy gyökvonás után megkapjuk, hogy
hányszorosára változott a sebesség abszolútértéke a manőver során. A jobb átláthatóság
kedvéért a
∣v VSA∣2 2 ∣v VSB∣
arányt P2-tel, a
∣v JSB∣ arányt pedig X-1-nel jelölöm, mert így a ∣v VSB∣
műholdnak a bolygóhoz képesti sebességének függvényében tudom a sebességnövekedést ábrázolni. P 2 X =2 X −21−2 X −1 X −21 cos 180 ˚−arctg X A P(X) függvény ábrázolva azt vehetjük észre, hogy minél nagyobb a sebessége a szondának a bolygóéhoz képest, annál kisebb mértékben növekszik a szonda sebessége. Határértékben 1hez tart az „utána”/„előtte” sebességek hányadosa. (5. ábra) 7
5. ábra: A szonda sebességének relatív növekedése
Bolygó neve
Átlagos keringési
Sebességnövekedés
sebessége [km/s] Merkúr
47.87
9.62
Vénusz
35.02
7.07
Föld
29.78
6.03
Jupiter 13.07 2.79 1. táblázat: A táblázat azt mutatja, hogy a fenti módszerrel (vagyis amikor a bolygó és a szonda sebessége kezdetben merőleges egymásra)
hányszoros
sebességnövekedést
lehet
elérni,
feltételezve, hogy a szonda sebessége kezdetben 10 km/s A 6. ábra egy konkrét esetét mutatja be a hintamanővernek: a Voyager-1 elhaladását a Jupiter mellett. Mint látható, a szonda 16 km/s extra sebességre tett szert a hintamanőver során. Az út során természetesen szükség lehet üzemanyagra kisebb pályamódosításokhoz. A 8
manőver nélkül viszont 5-6 nagyságrenddel több üzemanyagot kéne vinnie a szondának a gravitációs erők leküzdésére. Ezért a manőver mind gazdaságilag, mind energiafelhasználás szempontjából hasznos jelenség.
6. ábra: Hintamanőver a Jupiter mellett
9
A három- és többtest-probléma Egy n testből álló rendszerben egy adott testre ható erő a Newton-féle gravitációs erők n
F i=
eredője:
∑
mi m j r 3ij
j=1, j≠i
rij
Ebből a mozgásegynletek: mi r¨ i=
n
∑
j=1, j≠i
mi m j r 3ij
rij i=1,2,n (1)
n n mi m j Az égi mechanikában szokás bevezetni a következő függvényt: U = ∑ ∑ 2 i =1 j =1, j≠i r ij
Az U függvény segítségével az egyes komponensekre felírható a mozgásegyenlet: m i x¨ i=
∂U ∂U ∂U , mi y¨ i= , mi z¨i = (2). Az U függvény tehát a potenciál ellentettje. Ezek ∂ xi ∂ yi ∂ zi
az egyenletek 3n darab közönséges másodrendű differenciálegyenletet jelentenek, amiből meghatározandó x(t), y(t) és z(t). A differenciálegyenlet-rendszer megoldásának módja az első integrálok megkeresése. Első integrálnak azt a differenciálegyenletet nevezzük, amely f x , y , z , x˙ , y˙ , z˙ , t=const. alakú, vagyis a változóknak csak első deriváltjai szerepelnek bennük. A következő első integrálok egyszerűen levezethetők. Tömegközéppont-integrálok: összegezzük a mozgásegyenleteket! n
n
∑ mi r¨ i= ∑ i=1
n
∑
mi m j
i =1 j=1, j≠i
r 3ij
rij = 0 , mert
rij és
rji kiejti egymást a kettős összegben. n
Így a bal oldalt idő szerint kétszer integrálva a következőt kapjuk:
∑ mi r i=a tb
, ahol
i=1
a
és
b
konstans vektorok. A rendszer tömegközéppontját megadó vektor:
n
∑ mi r i r0=
i=1
m
n
, m=∑ mi
. A fentiek szerint tehát:
m r0= a t b , m r˙0 = a , ami a
i=1
tömegközéppont megmaradásának tétele: a tkp. vagy nyugalomban van, vagy egyenes vonalú egyenletes mozgást végez. Ez az egyenlet a rendszernek 6 első integrálját adja (3 koordináta, és 3 sebességkomponens). 10
Impulzusmomentum-integrál: (1) mindkét oldalát ri -vel vektoriálisan megszorozzuk, majd n
összegzünk:
n
∑ ri × mi r¨i = ∑ i=1
n
∑
i=1 j =1, j≠i
mi m j r 3ij
ri × rij =0 , mert az
ri × rij és
rj × rji
n
vektorok összege 0. Ezt idő szerint integrálva a következőt kapjuk:
∑ ri × mi r˙ij =c
. Az
i=1
egyenlet bal oldalán a rendszer össz impulzusmomentuma áll, ami a jobb oldal szerint állandó. Az impulzusmomentum vektor komponenseit felírva így újabb 3 első integrálhoz jutottunk. Energia-integrál: (1) mindkét oldalát r˙ i -tal skalárisan megszorozzuk, majd összegzünk n
n
∑ mi r¨ i r˙i = ∑ i=1
n
∑
i=1 j=1, j ≠i
mi m j r
3 ij
n
rij r˙i=∑ i=1
∂U ∂U ∂U x˙ y˙ z˙ . ∂ x i i ∂ yi i ∂ z i i
n
T=
Ha
bevezetjük
1 mi r˙ i2 jelölést (ami a kinetikus energia), akkor az egyenlet átírható ∑ 2 i=1
alakba. Ezt idő szerint integrálva:
a
dT dU = dt dt
T −U =h , vagyis a rendszer mechanikai energiája
állandó. Ez egy újabb első integrált ad. A 10 első integrált felhasználva az eredetileg 6n egyenletből álló rendszert át lehet írni (6n-10) differenciálegyenletből álló rendszerré. n=2 esetén egy másodrendű rendszert kapunk, melyet könnyű integrálni. n=3 esetén a rendszer nyolcadrendű lesz, aminek megoldásához további első integrálokra lenne szükség. Azonban Bruns és Poincaré bebizonyították, hogy nem lehet a meglévő első integráloktól független újabbakat felírni. Ez azt jelenti, hogy már három testre sem lehet az egyenleteket analitikusan megoldani. Az elméleti áttekintés után bemutatom azt a numerikus módszert, mellyel a differenciálegyenleteket oldom meg.
11
Numerikus differenciálegyenlet-megoldó módszerek A közönsége differenciálegyenletek megoldása sokszor nem írható fel zárt alakban. Ilyenkor a numerikus módszerekhez kell fordulnunk, hogy megoldhassuk a problémát. A többtest-szimulációkban közönséges differenciálegyenleteket kell megoldani, melyeknek ismert
a
kezdeti
értéke,
ezek
tehát
kezdetiérték-feladatok.
Bevezetem
a
y ' = f t , y és y t 0= y 0 jelöléseket. A megoldási módszer abból áll, hogy y(t) függvény közelítőértékeit számoljuk ki a ti időpontokban ( t i=t 0ih , h a lépésköz). t
A feladat átírható integrál-alakba:
y t = y0 ∫ f t ' , y t dt ' . Az első pontban a t0
t 0h
következő
közelítést
y t 1= y 0 ∫ f t , y t dt≈ y 0 h f t 0, y 0= y 1 .
végezzük:
t0
Általánosítva:
y i1= y ih f t i , y i . Ez az Euler-módszer. Ez a módszer h2 nagyságrendű
hibával közelíti a valódi értéket. A Runge és Kutta által kifejlesztett módszer pontosabb, mint az Euler-módszer, h5 a becslés hibája. A számolás menete a következő: Az
y(t)
függvényt
keressük,
ezt
később
a
bolygók
koordináta-
és
sebességkomponenseivel helyettesítjük. Ennek deriváltja viszont ismert, hiszen egy adott bolygó
elrendeződésnél
ismert
a
koordináta és a sebesség változása. A t n1=t nh
időpontban
y(t)-t
a
következővel közelítjük: h y n1= y n k 12k 22k 3k 4 , 6 ahol
7. ábra: Negyedrendű Runge-Kutta 12
k 1 = f t n , y n h h k 2 = f t n , y n k 1 2 2 h h k 3= f t n , y n k 2 2 2 k 4= f t nh , y n h k 3 Az A kiindulási pontban meghatározzuk a görbe meredekségét (k1). Az A-ból k1 meredekséggel az intervallum feléig lépünk (B). Meghatározzuk a B ponton keresztül menő görbe meredekségét (k2). Az A-ból k2 meredekséggel az intervallum feléig lépünk (C). Meghatározzuk a C ponton keresztül menő görbe meredekségét (k3). Az A-ból k3 meredekséggel az intervallum végéig lépünk (D). Meghatározzuk a D ponton keresztül menő görbe meredekségét (k4). Kiszámoljuk a meredekségek súlyozott összegét, ezzel a meredekséggel lépünk az A pontból az intervallum végére (E). Ez a pont lesz az eredeti görbe becsült értéke az intervallum végpontjában. Ezt a pontot a következő lépésben az új intervallum kezdőpontjának tekintjük. Hogy a Runge-Kutta módszert használom az Euler-módszer helyett, két oka van. Az egyik, hogy a kerekítési hiba az Euler-módszerrel sokkal nagyobb, mint egyéb, fejlettebb módszerekkel. A másik ok, hogy az Euler-módszer hajlamos numerikus instabilitásokra. Természetesen ettől függetlenül meg kell vizsgálni, hogy mennyire megbízható ez a módszer a bolygómozgás tanulmányozásában. Erről a következő fejezetben írtam.
13
A szimulációs módszer tesztelése A szimulációk elvégzése előtt nagyon fontos ellenőrizni, hogy analitikusan kiszámolható pályákat visszakapunk-e a szimulációval. Legelőször egy tömegpont Nap körüli mozgását vizsgáltam. A Nap fixen a (0;0;0) pontban van. Körülötte a kezdőfeltételek miatt körpályán mozog egy test. A test kezdeti koordinátái: (1;0;0), csillagászati egységekben mérve (csillagászati egység helyett a továbbiakban AU-t használok). A körpálya feltétele, hogy a centripetális erő megegyezzen a bolygó és a Nap közti gravitációs erővel: m
M m M Nap v2 = Nap2 v 0= x0 x0 4x 0
amibe a kezdeti feltételeket behelyettesítve kapjuk: v 0 =6.28
AU év
Az elméleti helynek a pozitív x tengelyhez mért szögét az ívhossznak a kerülethez képesti arányából számolom: =
vt vt 2= . A koordináták: 2 x 0 x0
x= x 0 cos
vt vt ; y= x 0 sin x0 x0
A pontok távolságát Pitagorasz-tétellel, a köztük lévő szöget pedig szintén az ívhossz és a kerület arányából számoltam. A 9. ábra mutatja a szögeltérés változását az időben. Az általam futtatott szimulációk hossza néhány hét, ez alatt az eltérés az elméleti és a számolt pálya közt 10-6 AU, ami gyakorlatilag nullának tekinthető. A másik kérdés a numerikus módszer idő invarianciájának megléte. A rendszert t=0ból fejlesztve t=t0-ig, majd innen vissza t=0-ig az eredeti kiinduló állapotba kell visszajutnia a rendszernek. Ezt úgy tesztelem, hogy egy nagy tömegű, mozgó bolygó mellett elküldök egy kis tömegű testet, majd a pálya egy tetszőlegesen kiválasztott pontjában megfordítom a sebességeket, és úgy küldöm vissza. Azt várom, hogy az oda és vissza irányú pálya megegyezzen. Azt vizsgáltam, hogy egy tetszőleges koordináta (legyen ez a 11. ábra vízszintes tengelye, továbbiakban x) mikor minimális (az ábrán látható, hogy van minimuma). Ezt a pontot P1-nek, illetve P-1-nek nevezem, attól függően, hogy melyik irányba fut a szimuláció. A visszafordulás idejét változtatva megkaphatjuk a visszatérés pontosságának a futási időtől való függését. Vagyis azt, hogy hosszabb ideig futtatva növekszik-e az eltérése az
14
oda és vissza pályának. Ennek eredménye a 12. ábra mutatja. Az látható, hogy a futási időtől lényegében nem függ P1 és P-1 távolsága, folyamatosan 0.00012 AU körül van, ami nagyjából 18 ezer km.
8. ábra: A körpálya tesztelése
9. ábra: A szögbeli eltérés időbeli változása
15
10. ábra: A szögbeli eltérés egy kinagyított részlete
11. ábra: Nagy tömegű bolygó melletti elhaladás (kinagyítva a manőver)
16
12. ábra: A minimum x helyek távolsága A mozgás másik legfontosabb tulajdonsága, hogy időben hogyan változik a sebesség. Ebben az esetben a szonda sebessége csökkent a kezdeti sebességéhez viszonyítva. Ez a két bolygóból álló elrendezés zárt rendszert alkot. Az összenergiát ábrázolva tehát egy konstans egyenest kell kapnunk. Az összenergia a szonda és a bolygó helyzeti és mozgási energiájának összegéből áll. A 13. ábra mutatja a szonda és a bolygó energiáját, valamint az összenergia 10-szeresét, hogy elkülönüljön a bolygó energiájától. Látható, hogy bár a szonda energiája csökken, ez elhanyagolható a bolygó energiájához képest. Az összenergia pedig konstans marad.
13. ábra: A rendszer energiája 17
A kezdeti paraméterek megváltoztatásával az eddigiektől teljesen eltérő eredményt is kaphatunk, ha a szonda a másik irányból kerüli meg a bolygót. A 15. ábrán a szonda kezdeti sebességének y irányú komponensét 1.5 AU/év-ről 1 AU/év-re változtattam. Ennek hatására a szonda nem jobbról, hanem balról kerüli meg a bolygót. A sebessége pedig nem csökken, hanem nő.
14. ábra: A sebesség időbeli változása
15. ábra: A megváltoztatott sebességkomponens hatására megváltozik a megkerülés iránya
18
16. ábra: A megváltoztatott sebességkomponens hatására a sebesség növekedett a manőver során A tesztelés eredménye a következő: a program képes visszaadni analitikusan számolható eredményeket, valamint bonyolultabb pályákat is képes számolni. Az idő invariancia megléte biztosítja, hogy ha kényelmi szempontok miatt az érkezés helyétől időben visszafelé akarom elindítani a számolást, akkor ugyanazt kapom, mintha a kezdőpontból indítottam volna időben előrefelé a szimulációt. Ezen eredmények alapján úgy ítélem meg, hogy a program alkalmas többtest-szimulációk elvégzésére. Ezért a továbbiakban is ezt a programot, illetve ennek némileg módosított változatát használom.
19
Valódi adatok alapján számolt pálya A Pioneer-11 űrszonda pályáját szimuláltam. Azért választottam a Pioneert, mert arról lehetett tudni, hogy a hintamanőver közben pályájába beavatkozás nem történt, így a programom képes szimulálni a mozgását. A pálya számítás előkészítéséhez szükségem volt a bolygók, és a Pioneer pályájának adataira. Ezeket a NASA honlapjáról3 tudtam letölteni. Az innen elérhető adatok nem teljes mértékben mérésekből származnak, hanem egy szimulációs program eredményei, de feltehetőleg összevetik időnként mérési adatokkal is. A továbbiakban ezekre mint valódi adatokra hivatkozom. Mivel a pálya számítás eredménye elég érzékeny a kezdeti paraméterek értékére, ezért a kiindulási pontot a manőver helyéhez közel választottam. Így a szimuláció kezdete 1974. november 28., és innen követtem 80 földi napon keresztül. A számolás során figyelembe vettem a Nap, a Föld, a Mars, a Jupiter és a Szaturnusz hatását. A többi bolygót azért hagytam ki a számításból, mert a Jupiter mellett ezek hatása elenyészően kicsi. A figyelembe vett bolygók mozgását nem számítottam, hanem a pályájuk adatait beolvastam, és interpoláltam az aktuális időpontra. A Pioneer mozgását viszont a programom számolta. Az alábbi ábrákon láthatók a futtatás eredményei.
17. ábra: A Pioneer-11 hintamanőverei az x-y síkra vetítve (valódi pályák) 3 http://ssd.jpl.nasa.gov/horizons.cgi
20
18. ábra: A Pioneer-11 hintamanőverei 3D-ben (1973. április 7. - 1980. december 31.)
19. ábra: A manőver kinagyított részlete 21
20. ábra: A manőver kinagyított részlete 3D-ben
21. ábra: A teljes szimuláció
22
A 21. ábrán látható a szimuláció tovább futtatása. Itt már a kezdeti feltételek, és egyéb számítási hibák miatt a számított pálya jelentősen eltér a valóditól. Emiatt nem találkozik a második hintamanőver helyén a Szaturnusszal, így maga a manőver sem jön létre. A kezdeti paraméterek (sebességkomponensek) módosításával sikerült ugyan a manőver helyéhez közelebb küldeni a szondát, viszont az odavezető úton, mivel az hosszabb lett, időkésést szenvedett, és nem a megfelelő időpontban találkozott a Szaturnusszal. A szimulációs programnak tehát a három térkoordináta mellett az időbeliséget is nagyon pontosan kell számolnia ahhoz, hogy visszakapjam a valódi pályát. A szonda mozgását nem csak a pályája jellemzi, hanem a sebessége is, aminek vizsgálata szintén nagyon fontos, hiszen az a manőver célja, hogy megváltoztassa a szonda sebességét, konkrétan ebben az esetben növelje. A 22. ábra elég jó egyezést mutat a valódi és a számolt sebességek közt. A kinagyított ábrán viszont az látszik, hogy a manőver során a legnagyobb sebességet nem egyidőben, hanem körülbelül 22 óra eltéréssel éri el a szonda a kétféle adatsor szerint. Ennek az időbeli eltérésnek később jelentős következménye az, hogy a Pioneer nem találkozik a Szaturnusszal, így a valódihoz képest egy teljesen másik pályán megy tovább. A manőver hatására az űrszonda sebessége kb. kétszeresére növekszik (9.5 km/s → 19 km/s). Érdekes megfigyelni, hogy a Pioneer kilövésekor volt a leggyorsabb, az útja többi részén, csak a kezdeti sebességénél kisebb sebességgel haladt.
22. ábra: A sebesség változása a futtatás ideje alatt 23
23. ábra: Kinagyított részlet az előző ábráról
24. ábra: A lépésköz változása a lépések során A szimulációhoz az adaptív Runge-Kutta módszert használtam. Most ennek az adaptív tulajdonsága érdekes, ugyanis a lépések során változik azok mérete. A lépés méretét akkor kell lecsökkenteni, ha a szonda olyan helyen tartózkodik, ahol nagyon erősen változik a gravitáció. Ez azért szükséges, mert a Runge-Kutta módszer akkor számol jól, ha a különböző 24
helyen vett deriváltak közel egyformák. Ez pedig a gyorsan változó gravitációs térben csak kis lépések esetén teljesül. Az algoritmusból következően a Jupiter mellett (200-300. lépés) 4 nagyságrenddel kisebbek a lépések (kb. fél perc), mint attól távol (10 földi nap).
25
Optimalizációs lehetőségek, újabb numerikus módszerek Az utóbbi időben több új numerikus módszert fejlesztettek ki bolygómozgás szimulációjára. Az egyik legfontosabb a szimplektikus integrátor (angolul sympectic integrator) módszere. Erre azért volt szükség, mert kellett egy megbízható módszer, ami a rendszer összenergiáját megtartja. Ennek a fejlesztése az 1980-as évek elején kezdődött: Ruth (1983) explicit módszert dolgozott ki a szétválasztható Hamilton-függvényű rendszerekre; Feng és Quin (1987) ugyanerre implicit módszert fejlesztett; Sanz-Serna (1988) és Lasagni (1988) pedig a Runge-Kutta módszert dolgozták át, hogy szimplektikusan viselkedjenek. Az alábbiakban röviden bemutatom, hogy hogyan működik a szimplektikus integrátorok módszere. A rendszer Hamilton-függvénye a kinetikus és potenciális energia összege: N N p 2i mj H =∑ −G ∑ mi ∑ , továbbá tudjuk, hogy i =1 2 mi i=1 j=i1 r ij N
dx i ∂ H dpi ∂H = , =− , ahol xi dt ∂ pi dt ∂ xi
az i-edik koordináta, pi az i-edik impulzus komponens.
q-val jelölve az általánosított
koordináták n-esét, felírható a következő egyenlet: n n dq ∂ q dx i ∂ q dp i ∂ q dH ∂q dH =∑ =∑ = F q , ahol dt i =1 ∂ x i dt ∂ pi dt i=1 ∂ x i dt ∂ pi dt
általános megoldása az egyenletnek: ahol
F
q t =e
F
egy operátor. Az
2 F 2 q t −=1 F qt − , 2
q t− egy kis időlépéssel korábbi állapota a rendszernek. A módszer lényege, hogy
a Hamilton-operátort felbontjuk egy összegre, melynek tagjai lehetőleg analitikusan H =H A H B
számolhatók. Például a Nap körül keringő testekre így írható fel: N
H A= ∑ i=1
2
p i G m Nap mi − , 2m i r i , Nap
N
H B=−G ∑
N
∑
i=1 j=i1
mi m j r ij
HA a Nap hatását írja le a körötte keringő testekre, HB pedig a testek egymásra hatását. A q(t)re vonatkozó egyenlet analitikusan megoldható HA behelyettesítésével, így pontos eredmény kapunk. Ugyanide HB-t helyettesítve csak numerikusan tudunk számolni. Ez azonban nem probléma, mivel a legtöbb esetben
H A≫ H B H B = H A , így elég az exponenciális első
két, legfeljebb három tagját figyelembe venni a számoláskor. Így
O 3 mértékű hibát
vétünk. Ez az eljárás akkor jó, ha nincsenek közeli elhaladások a bolygók közt, különben nem 26
teljesül a
H A≫ H B feltétel.
Az N-test szimulációk módszerei folyamatosan fejlődnek. A [4] szerint a közeljövőben tovább fognak fejleszteni egy olyan technikát, amivel nagy számú objektumot lehet a „dinamikai súrlódás” (dynamical friction) figyelembevételével szimulálni. A „dinamikai súrlódás” azt jelenti, hogy egy kis részecskékből álló anyagfelhőn keresztülhaladva a bolygó/ űrszonda a felhőnek energiát és impulzusmomentumot ad át. Súrlódásnak tehát az energia disszipációja miatt hívják a jelenséget. Ezt a módszert főleg galaxisok ütközésekor érdemes használni, mivel a galaxisban lévő csillagok egy „felhőnek” tekinthetők, amelyben mozog a másik galaxis (és viszont). Egy még összetettebb módszer azt az esetet kezeli, amikor nagyon sok, kis részecske koronggá áll össze, és ebben sűrűséghullámok alakulnak ki. A nagyszámú részecske szimulálását úgy érik el, hogy a teret fa-struktúrába rendezik, így összevonva az egymáshoz közeli részecskék hatását a távolabbiakra. A fa-struktúra miatt csak kell megoldani, nem
N log N egyenletet
N 2 -et, ahogy hagyományosan kéne. Ez pedig nagyon nagy számú
részecskénél jelentős csökkenést jelent. A módszert például a Naprendszer kialakulásának modellezésére lehet használni. A ma elfogadott egyik keletkezési elmélet szerint a Naprendszer egy forgó porfelhőből állt össze különálló égitestekké. A porfelhőt fel lehet fastruktúrába osztani, így gyorsítva a számolást. A hintamanővereket általában nem egy, hanem egymás után több bolygó mellett is kihasználják. A manőverek időpontjai előre le vannak rögzítve, eszerint kell megtervezni a pályát. Pályatervezéskor egy újabb szempont, hogy a szonda képes meghajtással módosítani a pályát. A meghajtás kétféle lehet: ritka, impulzusszerű és erősebb, vagy folyamatos, de gyengébb. Ezen túl vannak sebességi korlátok, valamint az, hogy az utazás végén egy adott méretű tértartományba érkezzen a szonda. A pálya így nagyon sok paramétertől függ. A paraméterek optimális megválasztása tehát egy sokdimenziós térbeli optimum keresésének felel meg, azzal kiegészítve, hogy több Lambert-problémát is meg kell oldani. Lambertproblémának nevezik azt, amikor egy gravitációs térben mozgó test pályáját kell úgy meghatározni, hogy az átmenjen a tér két meghatározott pontján. A pálya meghatározásának menete: 1. Meghatározni
x =T 0, T 1, T 2, ... ,T N 1 vektort, ahol
pedig az i-edik manőver időpontja (összesen N+1 darab). 27
T 0 az indulás ideje,
Ti
2. Megoldani N+1 darab Lamber-problémát. 3. Ellenőrizni a kapott pálya megvalósíthatóságát. 4. Kiszámolni a szonda tömegét minden lehetséges pályára, és ennek alapján kiválasztani a legalkalmasabbat. Optimalizációs eljárások közül nagyon sokfélét használnak, ezek részletezésétől most eltekintenék.
További fejlesztési lehetőségek Az eddigiekben azt sikerült megvalósítanom, hogy egy szondát elindítva, és magára hagyva számoltam a pályáját. Egy fejlesztési lehetőség a pályába való beavatkozás számítása. A beavatkozással (hajtómű bekapcsolás hosszabb-rövidebb időre) lehetséges megváltoztatni a pályát úgy, hogy még nagyobb legyen a manőver hatékonysága. Ehhez a fejlesztéshez szorosan kapcsolódik az optimum keresés feladata a pályára vonatkozó feltételek mellett (adott időben adott helyen, sebesség limit, bolygótól való megfelelő távolság tartása, üzemanyag felhasználás optimalizálása).
Technikai információk E dolgozat elkészítéséhez az alábbi eszközöket használtam: •
szövegszerkesztéshez:
OpenOffice 2.4
•
ábrák elkészítéséhez:
Gnuplot
•
program futási környezet:
Linux, Pentium PC
•
a számításokhoz szükséges programokat magam fejlesztettem C++ nyelven
•
fölhasználtam a publikusan elérhető RK adaptív modult
Köszönetnyilvánítás Ezúton szeretnék köszönetet mondani témavezetőmnek, dr. Csabai Istvánnak a szakdolgozatom elkészítésében nyújtott segítségéért, és útmutató tanácsaiért.
28
Irodalomjegyzék
1. Budó Ágoston. Mechanika. Budapest, Tankönyvkiadó, 1965. 2. Albert A. Bartlett, Charles W. Hord. “The slingshot effect: explanation and analogies.” The Physics Teacher 1985 november: 466-473 3. Érdi Bálint. Égi mechanika. Budapest, Nemzeti Tankönyvkiadó, 1996. 4. A. Morbidelli. “Modern Integrations of Solar System Dynamics.” Annu. Rev. Earth Planet. Sci. 30 (2002): 89-112 5. R. C. Johnson. “The Slingshot Effect.” (2003) 6. John D. Anderson, James K. Campbell, John E. Ekelund, Jordan Ellis, James F. Jordan. “Anomalous Orbital-Energy Changes Observed during Spacecraft Flybys of Earth.” Phys. Rev. Letters 100.091102 (2008): 1-4 7. J. E. Chambers. “A hybrid symplectic integrator that permits close encounters between massive bodies.” Mon. Not. R. Astron. Soc. 304 (1999): 793-799 8. Martin J. Duncan. “A multiple time step symplectic algorithm for integrating close encounters.” The Astronomical Journal 116 (1998): 2067-2077 9. Xin Wu, Tian-Yi Huang, Hong Zhang, Xiao-Sheng Wan. “A note on the algorithms of symplectic integrators.” Astrophysics and Space Sciences 283 (2003): 53-65 10. Vacheslav Vasilievitch Emel'yanenko. “A method of symplectic integrations with adaptive time-step for individual Hamiltonians in the planetary N-body problem.” Celestial Mech. Dyn. Astr. 98 (2007): 191-202 11. V. M. Becerra, D. R. Myatt, S. J. Nasuto, J. M. Bishop, D. Izzo. “An Efficient Pruning Technique for the Global Optimisation of Multiple Gravity Assist Trajectories.” Proceedings of GO (2005): 1-7 12. C. Kohlhase. “The Voyager Neptune Travel Guide.” NASA Jet Propulsion Laboratory (1989) 13. Roger E. Diehl. “Gravitational assist.” NASA Jet Propulsion Laboratory 14. A. F. B. A. Prado. “The Dynamics of the Gravity-assisted Maneuver.” Associação Brasileira de Ciências Mecânicas (1995) 15. C. R. H. Solórzano, A. A. Sukhanov, A. F. B. A. Prado. “A study of trajectories to the
29
Neptune system using gravity assist.” Advances in Space Research 40 (2007): 125-133 16. M. Vasile, P. De Pascale. “On the Preliminary Design of Multiple Gravity-Assist Trajectories.” Journal of Spacecraft and Rockets 43: 794-805 17. Chit Hong Yam, T. Troy McConaghy, K. Joseph Chen, James M. Longuski. “Design of Low-Thrust Gravity-Assist Trajectories to the Outer Planets.” 55th International Astronautical Congress (2004) 18. A Gravity Assist Primer
2009. május 17. 19. Gravity assist 2009. május 17. 20. Pioneer anomaly 2009. május 17. 21. N-test szimuláció 2009. május 17. 22. The
mystery
of
the
fly-by
anomaly
2009. május 17. 23. Runge-Kutta
eljárás
2009.
május 17. 24. Hintamanőver
leírás
2009. május 17. 25. Hintamanőver leírás 2009. május 17. 26. A pályaadatok helye 2009. május 17. 27. Runge-Kutta módszer leírás 28. Richard Fitzpatrick. Computational Physics: An introductory course. The University of Texas 29. D. Izzo, V. M. Becerra, D. R. Myatt, S. J. Nasuto, J. M. Bishop. “Search space pruning and global optimisation of multiple gravity assist spacecraft trajectories .” J. Glob. Optim. 38 (2007): 283-296 30. Dinamikai súrlódás leírása 2009. május 23. 31. J. N. Bronstejn, K. A. Szemengyajev, G. Musiol, H. Mühlig. Matematikai kézikönyv. Budapest, Typotex Kiadó, 2000.
30
Függelék A programkód: // // // // // // //
Pioneer-11 pályáját számoló program készítette: Puszta Adrián 2009-ben Fizika Bsc, III. évfolyam G = 1 m -> m * 1.982621412e-25 t -> t * 3.17098e-10 l -> l * 6.68458e-12
#include #include #include #include using namespace std; #include "vector.hpp" #include "odeint.hpp" using namespace cpl; const double pi = 4 * atan(1.0); const double G = 1.0/10000.0; const const const const const const
double double double double double double
mJup=376.4998062; mSat=112.7040968; m=5.134989458e-23; mNap=394363.2251; mMars=0.1272545553; mEarth=1.184457684;
// // // // // //
Jupiter Szaturnusz Pioneer Nap Mars Fold
Vector coords(12); double double double double double
linInterp(double, double, double, double); beolvasSaturn(double, int); beolvasJupiter(double, int); beolvasMars(double, int); beolvasEarth(double, int);
// Derivative vector for Newton's law of gravitation Vector f(const Vector& x) { double t = x[0], r_x = x[1], r_y = x[2], r_z = x[3], v_x = x[4], v_y = x[5], v_z = x[6]; double rSquaredJup = (r_x-coords[0])*(r_x-coords[0]) + (r_y-coords[1])*(r_y-coords[1]) + (r_z-coords[2])*(r_z-coords[2]); double rCubedJup = rSquaredJup * sqrt(rSquaredJup); double rSquaredSat = (r_x-coords[3])*(r_x-coords[3]) + (r_y-coords[4])*(r_y-coords[4]) + (r_z-coords[5])*(r_z-coords[5]); double rCubedSat = rSquaredSat * sqrt(rSquaredSat); double rSquaredMars = (r_x-coords[6])*(r_x-coords[6]) + (r_y-coords[7])*(r_y-coords[7]) + (r_z-coords[8])*(r_z-coords[8]); double rCubedMars = rSquaredMars * sqrt(rSquaredMars); double rSquaredSun = r_x*r_x + r_y*r_y + r_z*r_z; double rCubedSun = rSquaredSun * sqrt(rSquaredSun); double rSquaredEarth = (r_x-coords[9])*(r_x-coords[9]) + (r_y-coords[10])*(r_y-coords[10]) + (r_z-coords[11])*(r_z-coords[11]); double rCubedEarth = rSquaredEarth * sqrt(rSquaredEarth);
31
Vector f(7); f[0] = 1; f[1] = v_x; f[2] = v_y; f[3] = v_z; f[4] = - G* mJup * (r_x-coords[0]) / rCubedJup - G* mSat * (r_x-coords[3]) - G* mNap* r_x/ rCubedSun - G* mMars * (r_x-coords[6]) / rCubedMars - G* mEarth * coords[9]) / rCubedEarth; f[5] = - G* mJup * (r_y-coords[1]) / rCubedJup - G* mSat * (r_y-coords[4]) - G* mNap* r_y/ rCubedSun - G* mMars * (r_x-coords[7]) / rCubedMars - G* mEarth * coords[10]) / rCubedEarth; f[6] = - G* mJup * (r_z-coords[2]) / rCubedJup - G* mSat * (r_z-coords[5]) - G* mNap* r_z/ rCubedSun - G* mMars * (r_x-coords[8]) / rCubedMars - G* mEarth * coords[11]) / rCubedEarth;
/ rCubedSat (r_x/ rCubedSat (r_x/ rCubedSat (r_x-
return f; } double beolvasSaturn (double t, int n){ idopont; n koordinata (x:0, y:1, z:2) double beszam [4]; double betmp [4]; double ido; double tmp;
// beolvassa a file-bol a koordinatakat; t
ifstream be; be.open("Saturn_traj.data"); for (int j=0; j<4; ++j) be >> betmp[j]; while ((betmp[0] <= t) && !be.eof()) { for (int j=0; j<4; ++j){ be >> betmp[j]; } if (betmp[0]<=t){ for (int j=0; j<4; ++j){ beszam[j] = betmp[j]; } }else{ tmp = betmp[n+1]; ido = betmp[0]; } }; be.close(); return linInterp(ido, t, tmp, beszam[n+1]); } double beolvasJupiter (double t, int n){ idopont; n koordinata (x:0, y:1, z:2) double beszam [4]; double betmp [4]; double ido; double tmp;
// beolvassa a file-bol a koordinatakat; t
ifstream be; be.open("Jupiter_traj.data"); for (int j=0; j<4; ++j) be >> betmp[j]; while ((betmp[0] < t) && !be.eof()) { for (int j=0; j<4; ++j){ be >> betmp[j]; }
32
if (betmp[0]<=t){ for (int j=0; j<4; ++j){ beszam[j] = betmp[j]; } }else{ tmp = betmp[n+1]; ido = betmp[0]; } }; be.close(); return linInterp(ido, t, tmp, beszam[n+1]); } double beolvasMars (double t, int n){ idopont; n koordinata (x:0, y:1, z:2) double beszam [4]; double betmp [4]; double ido; double tmp;
// beolvassa a file-bol a koordinatakat; t
ifstream be; be.open("Mars_traj.data"); for (int j=0; j<4; ++j) be >> betmp[j]; while ((betmp[0] <= t) && !be.eof()) { for (int j=0; j<4; ++j){ be >> betmp[j]; } if (betmp[0]<=t){ for (int j=0; j<4; ++j){ beszam[j] = betmp[j]; } }else{ tmp = betmp[n+1]; ido = betmp[0]; } }; be.close(); return linInterp(ido, t, tmp, beszam[n+1]); } double beolvasEarth (double t, int n){ idopont; n koordinata (x:0, y:1, z:2) double beszam [4]; double betmp [4]; double ido; double tmp;
// beolvassa a file-bol a koordinatakat; t
ifstream be; be.open("Earth_traj.data"); for (int j=0; j<4; ++j) be >> betmp[j]; while ((betmp[0] <= t) && !be.eof()) { for (int j=0; j<4; ++j){ be >> betmp[j]; } if (betmp[0]
33
tmp = betmp[n+1]; ido = betmp[0]; } }; be.close(); return linInterp(ido, t, tmp, beszam[n+1]); } double linInterp (double tbegin, double t, double e, double b){ hely1, hely2 return b+(t-tbegin)*(e-b)*365.0; } int main() { double dt=0.0025, accuracy=1e-10; Vector x0(7); // kezdoertekek x0[0] = 600/365.0; x0[1] x0[2] x0[3] x0[4] x0[5] x0[6]
= = = = = =
4.87525; -0.769709; -0.111106; 0.00528*365; 0.003104*365*1.0045; 0.000204*365;
// idopont1, idopont2,
// PONTOSSAG
// t // // // // // //
x y z vx vy vz
Vector x = x0; int steps = 0; ofstream dataFile("k3n.data"); ofstream step("stepsize.data"); x = x0; steps = 0; int NN=0; double T = 3000/365.0; // SZAMOLASI HATAR double dt_max = 0, dt_min = 100; cout << "\n Integrating with adaptive step size" << endl; do { coords[0] = beolvasJupiter(x[0], 0); coords[1] = beolvasJupiter(x[0], 1); coords[2] = beolvasJupiter(x[0], 2); coords[3] = beolvasSaturn(x[0], 0); coords[4] = beolvasSaturn(x[0], 1); coords[5] = beolvasSaturn(x[0], 2); coords[6] = beolvasMars(x[0], 0); coords[7] = beolvasMars(x[0], 1); coords[8] = beolvasMars(x[0], 2); coords[6] = beolvasEarth(x[0], 0); coords[7] = beolvasEarth(x[0], 1); coords[8] = beolvasEarth(x[0], 2); for (int i = 0; i < 7; i++) dataFile << x[i] << '\t'; dataFile << coords[0] << "\t" << coords[1] << "\t" << coords[2]; dataFile << '\n'; double t_save = x[0]; adaptiveRK4Step(x, dt, accuracy, f);
// Runge-Kutta meghívása
34
step << x[0] << "\t" << x[0]-t_save << endl; double step_size = x[0] - t_save; ++steps; if (step_size < dt_min) dt_min = step_size; if (step_size > dt_max) dt_max = step_size; } while (fabs(x[0]) < T); cout << " number of adaptive steps = " << steps << endl; cout << " step size: min = " << dt_min << " max = " << dt_max << endl; cout << " data in file k3N.data" << endl; dataFile.close(); step.close(); }
35
NYILATKOZAT
Név: Puszta Adrián ELTE Természettudományi Kar, szak: Fizika BSc ETR azonosító: PUAOAAT.ELTE Szakdolgozat címe: A hintamanőver szimulációja
A szakdolgozat szerzőjeként fegyelmi felelősségem tudatában kijelentem, hogy a dolgozatom önálló munkám eredménye, saját szellemi termékem, abban a hivatkozások és idézések standard szabályait következetesen alkalmaztam, mások által írt részeket a megfelelő idézés nélkül nem használtam fel.
Budapest, 2009. május
_______________________________ a hallgató aláírása
