Design of Gravity Assist trajectory from Earth to Jupiter

The algorithm that is used is taken from Fundamentals of Astrodynamics and Applications [10]. This algorithm utilizes the bisection method which provide a strong solution for a wide variety of transfer orbits. The formulation of this method begins with the f and g universal variable, defined by the following formulas:

f=1−χ2r1C(z)g=Δt−χ3μC(z)g˙=1−χ3r2C(z)f˙=μr1r2χ(zS(z)−1)$$\begin{array}{} \displaystyle ~f=1-\dfrac{\chi^{2}}{r_{1}}C(z) \\\displaystyle ~g=\Delta t-\dfrac{\chi^{3}}{\sqrt{\mu}}C(z)\\\displaystyle~ \dot{g}=1-\dfrac{\chi^{3}}{r_{2}}C(z)\\\displaystyle \dot{f}=\dfrac{\sqrt{\mu}}{r_{1}r_{2}}\chi(zS(z)-1) \end{array}$$(1)

Where r₁ and r₂ are the magnitudes of the initial and final position vectors,

χ is a universal variable, z is square of the difference in eccentric anomalies, E, at two position (z = (ΔE)²)), C(z) and S(z) are define as:

C(z)=z−sin⁡zz3,S(z)=1−cos⁡zz$$\begin{array}{} \displaystyle C(z)=\dfrac{\sqrt{z}-\sin\sqrt{z}}{\sqrt{z^{3}}}, \quad S(z)=\dfrac{1-\cos\sqrt{z}}{z} \end{array}$$

The f and g expressions in terms of the orbital elements

g=r1r2μPsin⁡Δνf=1−r2P(1−cos⁡Δν)g˙=1−r1P(1−cos⁡Δν)f˙=μPtan⁡Δν2(1−cos⁡ΔνP−1r1−1r2)$$\begin{array}{} \displaystyle ~\,g=\dfrac{r_{1}r_{2}}{\sqrt{\mu P}}\sin\Delta\nu \\ \displaystyle~f=1-\dfrac{r_{2}}{P}(1-\cos\Delta\nu)\\\displaystyle ~\dot{g}=1-\dfrac{r_{1}}{P}(1-\cos\Delta\nu)\\\displaystyle \dot{f}=\sqrt{\dfrac{\mu }{P}}\tan\dfrac{\Delta\nu}{2}(\dfrac{1-\cos\Delta\nu}{P}-\dfrac{1}{r_{1}}-\dfrac{1}{r_{2}}) \end{array}$$(2)

Equating the corresponding equations in the two groups Eqs. (1) and Eqs. (2), we obtain

f=1−r2P(1−cos⁡Δν)=1−χ2r1C(z)$$\begin{array}{} \displaystyle f=1-\dfrac{r_{2}}{P}(1-\cos\Delta\nu)=1-\dfrac{\chi^{2}}{r_{1}}C(z) \end{array}$$(3)

g=r1r2μPsin⁡Δν=Δt−χ3μC(z)$$\begin{array}{} \displaystyle g=\dfrac{r_{1}r_{2}}{\sqrt{\mu P}}\sin\Delta\nu=\Delta t-\dfrac{\chi^{3}}{\sqrt{\mu}}C(z) \end{array}$$(4)

g˙=1−r1P(1−cos⁡Δν)=1−χ3r2C(z)$$\begin{array}{} \displaystyle \dot{g}=1-\frac{r_{1}}{P}(1-\cos\Delta\nu)=1-\dfrac{\chi^{3}}{r_{2}}C(z) \end{array}$$(5)

f˙=μPtan⁡Δν2(1−cos⁡ΔνP−1r1−1r2)=μr1r2χ(zS(z)−1)$$\begin{array}{} \displaystyle \dot{f}=\sqrt{\dfrac{\mu }{P}}\tan\dfrac{\Delta\nu}{2}(\dfrac{1-\cos\Delta\nu}{P}-\dfrac{1}{r_{1}}-\dfrac{1}{r_{2}})=\dfrac{\sqrt{\mu}}{r_{1}r_{2}}\chi(zS(z)-1) \end{array}$$(6)

We get from Eq. (3)

χ=r1r2(1−cos⁡Δν)PC(z)$$\begin{array}{} \displaystyle \label{5} \chi=\sqrt{\dfrac{r_{1}r_{2}(1-\cos\Delta\nu)}{PC(z)}} \end{array}$$(7)

Substituting with χ in Eq. (6) and cancelling h=μp$\begin{array}{} \displaystyle h=\sqrt{\dfrac{\mu}{p}} \end{array}$ from both sides, we obtain after simplification;

r1r2(1−cos⁡Δν)P=r1+r2−r1r2sin⁡Δν1−cos⁡Δν(1−zS(z))C(z)$$\begin{array}{} \displaystyle \dfrac{r_{1}r_{2}(1-\cos\Delta\nu)}{P}=r_{1}+r_{2} -\dfrac{\sqrt{r_{1}r_{2}}\sin\Delta\nu}{\sqrt{1-\cos\Delta\nu}}\dfrac{(1-zS(z))}{\sqrt{C(z)}} \end{array}$$(8)

We can write this equation more compactly by defining two auxiliary symbols, A and y as:

A=r1r2sin⁡Δν1−cos⁡Δν,y=r1r2(1−cos⁡Δν)P$$\begin{array}{} \displaystyle A=\dfrac{\sqrt{ r_{1}r_{2}}\sin\Delta\nu}{\sqrt{1-\cos\Delta\nu}}, \quad y=r_{1}r_{2}\dfrac{(1-\cos\Delta\nu)}{P} \end{array}$$

Using these definitions of A and y, Eqs. (7) and (8) may be written more compactly as

χ=yC(z)$$\begin{array}{} \displaystyle \chi=\sqrt{\dfrac{y}{C(z)}} \end{array}$$

y=r1+r2−A(1−zS(z))C(z)$$\begin{array}{} \displaystyle y=r_{1}+r_{2}-A\dfrac{(1-zS(z))}{\sqrt{C(z)}} \end{array}$$

If we now solve for Δt from Eq. (4), we get

μΔt=χ3S(z)+Ay$$\begin{array}{} \displaystyle \sqrt{\mu}\Delta t=\chi^{3}S(z)+A\sqrt{y} \end{array}$$

Using the auxiliary symbols A and y to write Eqs. (3 - 6) in the following simplified expressions:

f=1−yr1,g=Ayμ,g˙=1−yr2$$\begin{array}{} \displaystyle f=1-\dfrac{y}{r_{1}},\quad g=A\sqrt{\dfrac{y}{\mu}},\quad \dot{g}=1-\dfrac{y}{r_{2}} \end{array}$$

Then the solution of Lambert problems yields the following relations:

V1→=1g(r2→−fr1→),V2→=1g(g˙r2→−r1→)$$\begin{array}{} \displaystyle \overrightarrow{V_{1}}=\dfrac{1}{g}(\overrightarrow{r_{2}}-f\overrightarrow{r_{1}}),\quad \overrightarrow{V_{2}}=\dfrac{1}{g}(\dot{g}\overrightarrow{r_{2}}-\overrightarrow{r_{1}}) \end{array}$$

In order to escape the gravitational pull of a planet, the spacecraft must travel a hyperbolic trajectory relative to the planet, arriving at it’s sphere of influence with a relative velocity V_∞1 (hyperbolic excess velocity) greater than zero Fig. (2).

Fig. 2

The heliocentric velocity of S/C V0→$\begin{array}{} \displaystyle \overrightarrow{V_{0}} \end{array}$ at departure from the Earth is the sum of the Earth velocity VE→$\begin{array}{} \displaystyle \overrightarrow{V_{E}} \end{array}$ and the hyperbolic excess velocity V⃗_∞1.

V0→=VE→+V→∞1$$\begin{array}{} \displaystyle \overrightarrow{V_{0}}=\overrightarrow{V_{E}}+\overrightarrow{V}_{\infty 1} \end{array}$$

The latter is assumed to be equal to the spacecraft velocity relative to the Earth. In general it is VE→≫V∞1→$\begin{array}{} \displaystyle \overrightarrow{V_{E}}\gg \overrightarrow{V_{\infty 1}} \end{array}$ , due to the modest capabilities of present space propulsion, so that the maximum angle between VE→andV0→$\begin{array}{} \displaystyle \overrightarrow{V_{E}}~~ \text{and}~~ \overrightarrow{V_{0}} \end{array}$ is quite small. In particular the heliocentric leg will lie in a plane that can assume only a modest inclination away from the ecliptic plane.

The impulse required to be given at the perigee of the hyperbolic orbit to transfer the spacecraft from the parking orbit to the escape hyperbolic orbit is given by

ΔV1→=Vp→−Vc→$$\begin{array}{} \displaystyle \Delta\overrightarrow{V_{1}}=\overrightarrow{V_{p}}-\overrightarrow{V_{c}} \end{array}$$

Vp→$\begin{array}{} \displaystyle \overrightarrow{V_{p}} \end{array}$ is the velocity of the spacecraft at perigee of hyperbolic orbit,

Vc→$\begin{array}{} \displaystyle \overrightarrow{V_{c}} \end{array}$ is the velocity of the spacecraft in the parking orbit.

Clearly the direction of Vp→andVc→$\begin{array}{} \displaystyle \overrightarrow{V_{p}}~~ \text{and}~~ \overrightarrow{V_{c}} \end{array}$ are the same, then ΔV1→$\begin{array}{} \displaystyle \Delta\overrightarrow{V_{1}} \end{array}$ is also in the same direction.

We can obtain V0→$\begin{array}{} \displaystyle \overrightarrow{V_{0}} \end{array}$ from solving Lambert’s problem earth - Mars, see Table (1), we calculate VE→$\begin{array}{} \displaystyle \overrightarrow{V_{E}} \end{array}$ at T₀ to find V⃗_∞1 and it’s magnitude (V_∞1) which by it and a given perigee (r_p1 = R_E + 300) can calculate the hyperbolic trajectory elements. the angular momentum and eccentricity of the hyperbolic orbit can be obtained from the following relations [5]:

h1=rp1V∞12+2μrp1,e1=1+rp1V∞12μ1$$\begin{array}{} \displaystyle h_{1}=r_{p1}\sqrt{V_{\infty 1}^{2}+\dfrac{2\mu}{r_{p1}}}\quad,e_{1}=1+\dfrac{r_{p1}V_{\infty1}^{2}}{\mu_{1}} \end{array}$$

The velocity at the perigee of the hyperbolic orbit is:

VP1=V∞12+2μErp1$$\begin{array}{} \displaystyle V_{P1}=\sqrt{V_{\infty1}^{2}+\dfrac{2\mu_{E}}{r_{p1}}} \end{array}$$

The speed of S/c in its circular parking orbit is given by Vc=μ/rp1$\begin{array}{} \displaystyle V_{c}=\sqrt{\mu/r_{p1}} \end{array}$ Then the Δ V₁ required to put the S/C onto the hyperbolic departure trajectory is:

ΔV1→=Vp1→−Vc→=V∞12+2μErp1−μErp1$$\begin{array}{} \displaystyle \Delta\overrightarrow{V_{1}}=\overrightarrow{V_{p1}}-\overrightarrow{V_{c}}=\sqrt{V_{\infty1}^{2}+\dfrac{2\mu_{E}}{r_{p1}}}-\sqrt{\dfrac{\mu_{E}}{r_{p1}}} \end{array}$$

ΔV1→=Vc2+V∞1Vc2−1$$\begin{array}{} \displaystyle \Delta\overrightarrow{V_{1}}=V_{c}\left( \sqrt{2+\left( \dfrac{V_{\infty1}}{V_{c}}\right)^{2}}-1\right) \end{array}$$

The orientation of the apse line of the hyperbola to the asymptotes of the hyperbolic trajectory measured by the angle, which can be obtained from the relation [5]

β=cos−11e1=cos−111+rp1V∞12μE$$\begin{array}{} \displaystyle \beta=cos^{-1}\left( \dfrac{1}{e_{1}}\right) =\cos^{-1}\left( \dfrac{1}{1+\dfrac{r_{p1}V_{\infty1}^{2}}{\mu_{E}}}\right) \end{array}$$

The results are summarized in the Table (2).

Table 2

Escape from Earth at T₀

V⃗_E0 = {15.113446216468128, -25.49053048920682, 0.0008860019561318039}

V⃗₀ = {18.617165466382446, -28.29950136444239, -1.1521167643597399}

V⃗_∞1 = {3.50372, -2.80897, -1.153}

JD (Julian day No) of T₀ = 2458230.5

V_{(parking orbit)} = 11.8689(km/s)

Now after solving Lambert’s problem for earth - Mars trajectory and Mars, Jupiter trajectory we have V⃗₁ and V⃗₂, see Table (1). Such that:

V⃗₁ is the heliocentric S/C velocity at final position for Earth-Mars Lambert’s algorithm,

V⃗₂ is the heliocentric S/C velocity at initial position for Mars- Jupiter Lambert’s algorithm,

V⃗_M is the velocity of Mars at time T₁.

Then the heliocentric velocity of spacecraft at the SOI of Mars is V⃗₁ and it need to out from the SOI by V⃗₂, to complete the trajectory and finally reach to the SOI of Jupiter, on the other hand the spacecraft enter the SOI velocity relative to mars is V→∞(1)$\begin{array}{} \displaystyle \overrightarrow{V}_{\infty}^{(1)} \end{array}$ and it need to out with velocity relative to Mars is V⃗_outFig. (3).

V→∞(1)=V→1−V→MV→out=V→0−V→M$$\begin{array}{} \displaystyle \overrightarrow{V}_{\infty}^{(1)}= \overrightarrow{V}_{1}- \overrightarrow{V}_{M}\\ \displaystyle\overrightarrow{V}_{out}= \overrightarrow{V}_{0}- \overrightarrow{V}_{M} \end{array}$$

We take the direction of V→∞(2)$\begin{array}{} \displaystyle \overrightarrow{V}_{\infty}^{(2)} \end{array}$ (the out velocity of spacecraft in hyperbolic orbit fly by) in the same direction of V⃗_out then the turn angle of the hyperbolic orbit flyby can be calculate from this relation [5],

cosδ=V→∞(1).V→∞(2)V∞2=V→∞(1).V→outV∞V→out$$\begin{array}{} \displaystyle cos\delta=\dfrac{\overrightarrow{V}_{\infty}^{(1)}.\overrightarrow{V}_{\infty}^{(2)}}{V_{\infty}^{2}}=\dfrac{\overrightarrow{V}_{\infty}^{(1)}.\overrightarrow{V}_{out}}{V_{\infty}\left|\overrightarrow{V}_{out}\right| } \end{array}$$

Fig. 3

We know that

V→∞(1)=V→∞(2)$$\begin{array}{} \displaystyle \left|\overrightarrow{V}_{\infty}^{(1)}\right| =\left| \overrightarrow{V}_{\infty}^{(2)}\right| \end{array}$$

Now we can calculate the hyperbolic orbital elements using the relation [5],

e2=1sin⁡(δ2),rp2=μMV∞22(e2−1),Δ=rp1+2μMrpV∞22$$\begin{array}{} \displaystyle e_{2}=\dfrac{1}{\sin(\frac{\delta}{2})},\quad r_{p2}=\dfrac{\mu_{M}}{V_{\infty 2}^{2}}(e_{2}-1),\quad\Delta =r_{p}\sqrt{1+\dfrac{2\mu_{M} }{r_{p}V_{\infty2}^{2}}} \end{array}$$

Then we can calculate the perigee of the hyperbolic fly by Mars, now the spacecraft out the SOI of Mars with the velocity V→∞(2)$\begin{array}{} \displaystyle \overrightarrow{V}_{\infty }^{(2)} \end{array}$ relative to Mars with magnitude is V_∞2 and in the same direction of V⃗_out, to make the velocity of spacecraft relative to Mars is V⃗_out we give it Δ V₂ in the same direction of V→∞(2)$\begin{array}{} \displaystyle \overrightarrow{V}_{\infty }^{(2)} \end{array}$,

ΔV2=V→out−V∞2$$\begin{array}{} \displaystyle \Delta V_{2}= \left| \overrightarrow{V}_{out}\right|-V_{\infty 2} \end{array}$$

After that the spacecraft out from the SOI of Mars with heliocentric velocity is V⃗₂ = V⃗_M + V⃗_out

which by it can complete it’s trajectory to reach the SOI of Jupiter. The results are summarized in Table (3).

Table 3

Fly by Mars at T₁

V⃗_in = V⃗₁ = {21.72542643031754, 13.844699399653631, 0.013528429195021072}

V⃗_out = V⃗₂ = {30.823073404118684, 4.934176592416298, -0.8994633203622276}

V⃗_M = {16.48592538172737, 20.64382791494634, 0.02774276866812586}

JD (Julian day No) of T₁ = 2458230.5

A spacecraft arrives at the sphere of influence of the Jupiter with a hyperbolic excess velocity V⃗_∞3 relative to Jupiter, where

V→∞3=V→3−V→J$$\begin{array}{} \displaystyle \overrightarrow{V}_{\infty3}=\overrightarrow{V}_{3} -\overrightarrow{V}_{J} \end{array}$$

Such that:

V⃗₃ is the heliocentric S/C velocity at final position for Mars- Jupiter lambert algorithm. see Table (1), or it is the heliocentric S/C velocity at SOI of Jupiter at T₂

V⃗_J is the heliocentric velocity of Jupiter at T₂.

the goal’s mission is landing on Jupiter To achieve this goal we make the perigee of the hyperbola equal to the Jupiter radius (r_p3 = R_J = 69911 km) and the velocity of S/C relative Jupiter equal zero at perigee. we give S/C the third impulse in the perigee of the hyperbola equal to the it’s velocity and in the opposite direction. The results are summarized in Table (4),

Table 4

Capture from Jupiter at T₂

V⃗₃ = {-5.888427200446674, 3.2105733597573787, 0.22569580208145187}

V⃗_J = {9.130213923868334, 9.803076409849977, -0.2456201552966279}

JD (Julian day No) of T₂ = 0.21253784951961086

V⃗_∞3 = {-15.0186, -6.5925, 0.471316}

ΔV→3=−V→p3$$\begin{array}{} \displaystyle \Delta\overrightarrow{V}_{3}=-\overrightarrow{V}_{p3} \end{array}$$