Application Center - Maplesoft

App Preview:

Hohmann Elliptic Transfer Orbit with Animation

You can switch back to the summary page by clicking here.

Learn about Maple
Download Application




Hohmann Elliptic Transfer Orbit

Dr. A. C. Baroudy
acbaroudy@yahoo.com

Abstract

The main purpose of this article is to show how to use Hohmann elliptic transfer in two situations:
                       a- When one manned spaceship is trying to catch up with an other one
                            on the same circular orbit around Earth.
                       b- When delivering a payload from Earth to a space station on a circular
                            orbit  around Earth using 2-stage rocket .

The way we set up the problem is as follows:
                             Consider two manned spaceships with astronauts Sally & Igor , the latter
                             lagging behind Sally by a given angle = 4.5 degrees while both are on the same
                            circular orbit C2 about Earth. A  2d lower circular orbit  C1 is given.
                            Find the Hohmann elliptic orbit that is tangent to both orbits which allows
                            Sally to maneuver on C1 then to get back to the circular orbit C2 alongside Igor.

Though the math was correct , however the final result we found was not !!
It was somehow tricky to find the culprit!
We have to restate the problem to get the correct answer.
The animation was then set up using the correct data.
The animation is a good teaching help for two reasons:
                       1- it gives a 'hand on' experience for anyone who wants to fully understand it,
                       2- it is a good lesson in Maple programming with many loops of the type 'if..then'.

                       Warning

                       This particular animation is a hog for the CPU memory since data accumulated
                       for plotting reached 20 MB! This is the size of this article when animation is  
                       executed. For this reason and to be able to upload it I left the animation  
                       procedure non executed which drops the size of the article to 350KB.

                       Conclusion

                      If  I can get someone interested in the subject of this article in such away that he or
                      she  would seek further information for learning from other sources,  my efforts
                      would be well rewarded.
                      
                                                                           ______________








Problem #1:

We consider two manned spaceships with astronauts Sally & Igor moving on the same circular orbit of radius  
                                      r2 =  (6378 + 1200) km
around the Earth (Radius = 6378 km) as depicted in Figure-1 (a) &(b).
Igor is lagging behind Sally by 4.5 degrees as seen in Figure-1 (c).
Sally has to adjust its course so that it catches up with Igor. How can she do it?

Problem #2:

A space station S is moving on a circular orbit of radius  r2 =  (6378 + 1200) km around the Earth (Radius = 6378 km). A payload of 500 kg is to be delivered to it from Earth. How can we do it?
Both problems can be solved using the Hohmann elliptic orbit transfer.


                                         Initial settings with the constants

R[Earth] = 6.378*10^6
M[Earth] = 5.9737*10^24
G =constant of gravitation = 6.672590000*10^(-11),
μ =G*M[Earth] = 3.986005088*10^14,
h1 = 450 km, this is the altitude of the smaller circle , its radius becomes:  r[mini] = (6378+450)*10^3 = 6.828*10^6 m.
h2 = 1200 km, this is the altitude of the larger circle, its radius becomes :   r[max] = (6378+1200)*10^3= 7.578*10^6 m.

pi := evalf(Pi); R[Earth] := 6.378*10^6
M[Earth] := 5.9737*10^24; G := 6.672590000*10^(-11); mu := G.M[Earth]; h1 := 450*10^3; h2 := 1200*10^3; r[maxi] := h2+R[Earth]; r[mini] := h1+R[Earth]; r[2] := r[maxi]; r2 := r[maxi]; r[c2] := r[maxi]; r1 := r[mini]; r[1] := r[mini]; r[c1] := r[mini]; v[c2] := sqrt(mu/r[maxi]); v[c1] := sqrt(mu/r[mini]); TC1 := 2*pi*r[mini]/v[c1]; TC2 := 2*pi*r[maxi]/v[c2]; a[E] := (r[maxi]+r[mini])*(1/2); e := (r[maxi]-r[mini])/(r[maxi]+r[mini]); V[perigee] := sqrt(mu*(1+e)/(a[E]*(1-e))); V[apogee] := sqrt(mu*(1-e)/(a[E]*(1+e))); MeanMotion := sqrt(mu/a[E]^3); tau := 2*evalf(pi)/MeanMotion; h := r[mini]*V[perigee]; L := h^2/(G*M[Earth])

3.141592654

 

6378000.000

 

0.5973700000e25

 

0.6672590000e-10

 

0.3986005088e15

 

450000

 

1200000

 

7578000.000

 

6828000.000

 

7578000.000

 

7578000.000

 

7578000.000

 

6828000.000

 

6828000.000

 

6828000.000

 

7252.564901

 

7640.506827

 

5615.018774

 

6565.122672

 

7203000.000

 

0.5206164098e-1

 

7836.872140

 

7061.251382

 

0.1032758557e-2

 

6083.885982

 

0.5351016297e11

 

                                       Elements of the small circular orbit
r[c1] := h1+R[Earth] = 6.828.10^6 m,
v[c1] := sqrt(mu/r1) = sqrt(3.986005088*10^14/(6.828*10^6))
                                       Elements of the large circular orbit
     r[c2] := h2+R[Earth] = 7.578.10^6 m,
     v[c2]:=                                         Elements of the elliptic orbit
The semi-major axis:
   a[E] := (r[c1]+r[c2])*(1/2)= (r1+r2)/(2)=(1/2)*(6.828*10^6+7.578*10^6) = 7.203000000*10^6m.
Eccentricity:
   
e:=((r[max]-r[mini]))/((r[max]+r[mini]))=((6.828*10^(6)- 7.578*10^(6)))/((6.828*10^(6)+7.578*10^(6)))=0.05206164098.  V[perigee]:=sqrt((2 mu)/()*(1/(r[mini])-1/(2 a)))=sqrt((2 mu)/()*(1/(r[mini])-1/(r[max]+r[mini])))=sqrt((2 mu)/()*((r[max])/(r[mini] (r[max]+r[mini]))))=sqrt((2 mu)/()*((r[max])/(r[mini] 2 a))).
V[perigee] := sqrt(mu*(1+e)/(a[E]*(1-e))) = 7836.872141 m/s.
V[apogee] := sqrt(mu*(1-e)/(a[E]*(1+e))) = V[perigee]*(1-e)/(1+e) and V[perigee]*(1-e)/(1+e) = 7061.251382 m/s.
MeanMotion := sqrt(mu/a[E]^3) = 0.1036808256e-2.
 tau := 2*evalf(pi)/MeanMotion = 2*Pi*sqrt(a[E]^3/mu) and 2*Pi*sqrt(a[E]^3/mu) = 6060.122758 seconds.
 h := r[mini]*V[perigee] = 5.343253258*10^10  kg. m^2/s.                  
L := h^2/(G*M[Earth])
 =                                         Problem #1 in details

                                                 Sally Maneuver

To inject the satellite into elliptic orbit from point B (Figure-1 (a)) on the large circular orbit r2 = r[max] ( apogee of the elliptic orbit)  Sally needs to decrease her velocity from
                                           v[c2] = 7252.564901  to V[apogee] = 7061.251382:
so that the decrement of the velocity is:
                                           `ΔV`[2] := v[c2]-V[apogee] = 191.313518m/s.
To inject the satellite into small circular orbit at point A she has to reduce her velocity by:
                                          `ΔV`[1] := V[perigee]-v[c1] = 196.365312 m/s.
When Sally reaches point A, Igor would be at point i[1] (Figure-1 (a)).
The difference in degrees between Igor at i[1] and Sally at A is angle AOi[1]:
                           AOi[1] := 2*((1/2)*T-(1/2)*tau)*evalf(Pi)/T+(4.5*(1/180))*Pi = .3088248736*rad,
                                             13.19436186                     +   4.5       = 17.694361862degrees.
While on the smaller circular orbit Sally has to make a few revolutions waiting for Igor to get into the right relative position I[2] symmetrical with position i[2] (AOi[2] = AOI[2] and AOI[2] = 13.19436186 degrees) before firing again her spaceship to inject it into the transfer ellipse.
The gap between Sally at A after a few revolutions on the smaller circle and Igor at I[2] is
         the small angle (I[2]Oi[1])=17.694361862+13.19436186=30.88872372 degrees.
The gap between them is in fact the large angle (I[2]Oi[1])= 360 - 30.88872372=329.1112763 degrees.
A 360 degrees revolution on the small circular orbit corresponds to a fraction of a revolution on the large circular orbit for Igor. The difference in degrees is :
                                         360-360*TC1/TC2 = 52.0991642 degrees.
The number of revolutions that Sally has to go through till Igor gets to point I[2] is:
                                       329.1112763/(52.0991642) = 6.317016429 revolutions.
When the time of these revolutions has elapsed:
                                       6.317016429*TC1*(1/60) = 9.852823847 hours,
Sally should start firing at this moment in order to go on elliptic orbit again at A. For this to happen she needs to increase her velocity from  
                                           v[c1] = 7640.506827  to V[perigee] = 7836.872141:
hence the increment of the velocity is:
                                          `ΔV`[1] := V[perigee]-v[c1] = 196.365312 m/s.
She will coast along this elliptic orbit while Igor is moving along the large circular orbit. When their paths are parallel ( hence having the same potential energy) which occurs at B, she has to fire her rocket to increase her forward speed by:
                                         `ΔV`[2] = 191.313518m/s,
and she would sail alongside Igor with no lagging between them.


                                       

                       &nb sp;                What was wrong with problem #1?

The problem that we just solved doesn't give the correct answer for Sally to be at point A when Igor is at point I[2] on his orbit.
If the number of revolutions were an integral number  then Sally would've been back into A when Igor got at point I[2] on his orbit. However the number of revolutions was not an integral number it was found to be
                                   6.317016429 = 6 revolutions  +  .317016429 of a revolution
                                                         =  6 revolutions + 114.1259144 degrees.
Hence when Igor was at the right position at  I[2] , Sally was not at A . She has already passed this point by 114 degrees.
Firing her forward rocket won't get her at point B: her new point of contact with the large circular orbit will be ahead of Igor by 114 degrees!
This error is the result of the constraint  that comes from setting the radius of the smaller circular orbit to a fixed value ( 450 km) in advance. We should leave this radius as unknown, say (x), then solve our problem till we get to the expression that gives the number of revolutions as function of (x). We then set this expression as equal to an integral number, say 6, then with the help of  Maple to 'fsolve' it for (x).
Doing so we get                                   h[1] = 412463.0090 meters.
Now if restating the problem by specifying that
                                                r[mini] = 6378*10^6+412463.0090 meters,
then the number of revolutions that Sally has to go through on the smaller circular orbit till Igor gets to point I[2] on his orbit is:
                                                 5.9999999 ≈ 6 revolutions.
Now Sally can be sure that rendezvous with Igor will occur at point B.

                                         Finding the time Δt during which the engine should be fired
Knowing
1- the velocity change each time Sally has to pass form one orbit to another in m/s,
2- the thrust that her rocket (or retrorocket) can deliver in Newtons ,
3- the mass of the spaceship in kg,
Sally will use the impulse - linear momentum change  to find the necessary time during which she has to fire her rocket.
From this relation :
                                                  f*`Δt` = m*`Δv`,
we get :                                    `Δt` = `mΔ`/f.

                                        Problem #2

Here we use two-stage rocket to raise the payload to the apogee (point B in Figure(1)(a)) of an initial elliptic orbit which is supposed to be tangent at B to the circular orbit of the station. Then at this point the second stage is ignited to increase its velocity to that of the large circular orbit to put the payload on it. The timing for both stages is very important as we found above in order that the rendezvous with the space station occurs.

with(plottools); with(plots); h1 := 4.12463*10^5; r[mini] := h1+R[Earth]; r[2] := r[maxi]; r2 := r[maxi]; r[c2] := r[maxi]; r1 := r[mini]; r[1] := r[mini]; r[c1] := r[mini]

 

412463.0000

6790463.000

7578000.000

7578000.000

7578000.000

6790463.000

6790463.000

6790463.000

(2)

 

a[E] := (r[maxi]+r[mini])*(1/2); e := (r[maxi]-r[mini])/(r[maxi]+r[mini]); V[perigee] := sqrt(mu*(1+e)/(a[E]*(1-e))); V[apogee] := sqrt(mu*(1-e)/(a[E]*(1+e))); h := r[mini]*V[perigee]; L := h^2/(G*M[Earth]); v[c2] := sqrt(mu/r[maxi]); TC2 := 2*pi*r[maxi]/v[c2]; v[c1] := sqrt(mu/r[mini]); TC1 := 2*pi*r[mini]/v[c1]; tau := 2*pi/sqrt(mu/a[E]^3)

 

7184231.500

0.5481010739e-1

7868.761316

7051.007204

0.5343253257e11

7162649.003

7252.564901

6565.122672

7661.595702

5568.779536

6060.122758

(3)

 

Kepler equation to get:
                                           E =  eccentric anomaly, 
then we get θ from the relation :
                                         cos(E) = (e+cos(theta))/(1+e*cos(theta)) where θ = true anomaly.

Kepler := E-e*sin(E)-t*sqrt(G*M[Earth]/a[E]^3) = 0

E-0.5481010739e-1*sin(E)-0.1036808256e-2*t = 0

(4)

If in Kepler's equation we replace ( t ) by a small interval Δt = period/100:
                                                              
                                        `Δt` = (1/100)*P and (1/100)*P = (2*Pi*(1/100))*sqrt(a[E]^3/(G*M[Earth])) 
we get E then angle θ and finally r (the position of the satellite) after this interval of time. We then get these data for (1/100)*P*2 , (1/100)*P*3, (1/100)*P*4,..., (1/100)*P*n,...etc till one complete period where n = 100. Plotting these positions gives the elliptic orbit. With this substitution the equation becomes:
                 E-0.5481010740e-1*sin(E)-2*pi*((1/100)*n) = 0.
This is a transcendental equation that is solved using Newton's iterative method.

E-0.5481010740e-1*sin(E)-2*pi*((1/100)*n)

Plot of the three orbits: the elliptic, the small circular orbit & the large circular orbit

plot([(1/1000)*r[mini]*cos(`ϕ`), (1/1000)*r[mini]*sin(`ϕ`), `ϕ` = 0 .. 2*pi], color = red); C1 := %; plot([(1/1000)*r[maxi]*cos(`ϕ`), (1/1000)*r[maxi]*sin(`ϕ`), `ϕ` = 0 .. 2*pi], color = gold); C2 := %; plot([0.1e-2*L*cos(`ϕ`)/(1+e*cos(`ϕ`+(1/2)*pi)), 0.1e-2*L*sin(`ϕ`)/(1+e*cos(`ϕ`+(1/2)*pi)), `ϕ` = 0 .. 2*pi], color = blue); EL := %; display(C1, C2, EL)

 

                                       This is the procedure for the animation.

This procedure has so many points to plot that the total size of the file went from 300 kilobytes to 20 megabytes! For this reason I left this procedure non executed in order to reduce the size of the file to a mere 300kB for easy uploading.
To see the animation just remove (#) in the last two commands then execute them. 
Be patient since the loading of the data for plotting takes at least 20 to 30 seconds on my computer.
The animation starts when Sally, being at B and Igor lagging behind her by 4.5 degrees, fires her spaceship retrorocket to move into elliptic orbit apogee.
At A she moves into smaller orbit by reducing her speed ( using retrorocket again) to that of the smaller circular orbit (red).
She makes exactly 6 revolutions (count them) on the small circular orbit to be at A when Igor is at the right point I[2]on his larger circular orbit (gold).
At A she uses her forward rocket to put her spaceship into elliptic orbit (blue) again going to point B which she reaches at the same time as Igor.
At B she uses again her forward rocket to get on the larger circular orbit alongside Igor.
I use different colors for the spaceships ( blue for Sally & red for Igor) as well as for the 3 different orbits (blue for Sally on elliptic then red on the smaller circuar orbit & gold for Igor orbit). Enjoy it.

r1; r2; TC1; TC2; tau; e; h; L

6790463.000

 

7578000.000

 

5568.779536

 

6565.122672

 

6060.122758

 

0.5481010739e-1

 

0.5343253257e11

 

7162649.003

(5)

sat := proc (n) local i, E1, E2, r, j, k, theta, GEL, GC1, GC2; global pt1, D2, D1, D3; j := n; k := n-650; E1 := 0.6283185308e-1*j; if 50 <= n and n <= 650 then j := 50 end if; if 650 < n then j := 50+k end if; for i to 50 do E2 := E1+(-1)*0.5481010740e-1*sin(E1)+(-1)*0.6283185308e-1*j; if abs(E2) < 1/1000000 then break end if; E1 := E1-E2/(1-e*cos(E1)) end do; theta := arccos((e-cos(E1))/(e*cos(E1)-1)); if pi < E1 and E1 < 2*pi then theta := 2*pi-theta end if; r := L/(1+e*cos(theta+pi)); pt1 := [r*cos(theta), r*sin(theta)]; if 50 <= n then GEL := disk([-(1/1000)*pt1[2], (1/1000)*pt1[1]], 1, color = blue) else GEL := disk([-(1/1000)*pt1[2], (1/1000)*pt1[1]], 150, color = blue) end if; if 650 < n then GEL := disk([-(1/1000)*pt1[2], (1/1000)*pt1[1]], 150, color = blue) end if; if 50 <= n then GC1 := disk([(1/1000)*r1*cos((1/50)*Pi*n+(1/2)*pi), (1/1000)*r1*sin((1/50)*Pi*n+(1/2)*pi)], 150, color = blue) else GC1 := disk([0, -(1/1000)*r1], 1, color = red) end if; if 650 < n then GC1 := disk([0, -(1/1000)*r1], 1, color = red) end if; if n <= 50 then GC2 := disk([(1/1000)*r2*cos((1/50)*Pi*n*tau/TC2+(1/2)*pi-0.7853981635e-1), (1/1000)*r2*sin((1/50)*Pi*n*tau/TC2+(1/2)*pi-0.7853981635e-1)], 150, color = red) else GC2 := disk([(1/1000)*r2*cos((1/50)*Pi*(n-50)*TC1/TC2+(3/2)*pi+(-1)*18.34589276*pi/180), (1/1000)*r2*sin((1/50)*Pi*(n-50)*TC1/TC2+(3/2)*pi+(-1)*18.34589276*pi/180)], 150, color = red) end if; if 650 < n then GC2 := disk([(1/1000)*r2*cos((1/50)*Pi*(n-650)*tau/TC2+(3/2)*pi+13.84589277*pi/180), (1/1000)*r2*sin((1/50)*Pi*n*tau/TC2+(3/2)*pi+13.84589277*pi/180)], 150, color = red) end if; display(GEL, GC1, GC2, C1, C2, EL) end proc

n := seq(k, k = 0 .. 700)

NULL

NULL

``