Math- Calculating delta V for an interplanetary transfer

[Amagi82]

I'm having some difficulty figuring out the math for transferring from orbit around one planet/moon to the orbit around another. I know about all the charts that provide the approximate values, but I'd like to know the mathematical equations necessary to calculate the exact values.

So far, using the Vis-viva equation, along with Hohmann transfer calculations, I've managed to get accurate values from the surface to orbit around a body at a specific altitude, with or without atmosphere, I can calculate orbital velocity at a specific altitude, I can calculate the speed I need to be traveling when I leave the SOI to reach another planet/moon, and I know the speed I will be traveling when I enter the SOI of the target body. What I have not been able to calculate is exactly how much delta V I'll have to expend at a given periapsis altitude to escape at the desired velocity.

Re-arranging the Vis-viva equation, with a as the semi-major axis, I get:

a= 1/((2/r)-(v^2/u)).

But this is giving me incorrect values for a from what I know to be the correct answer. Anybody know the math and willing to walk me through it here? Thanks!

Edited May 5, 2014 by Amagi82

[Supernovy]

If I'm understanding the question correctly, this should help: http://en.wikipedia.org/wiki/Hyperbolic_trajectory#Hyperbolic_excess_velocity

Setting the vis-viva SMA to the same as the one given by that should give you the correct periapsis velocity. The Δv will be this velocity minus circular orbit velocity at the periapsis altitude.

[Yasmy]

Let's go from planet A to planet B, where planet B orbits farther from the sun S.

Let small first letters denote the ship, and capitol first letters denote the planets:

rA = the ship's radius around planet A

RA = the planet's radius around the sun S

vA = the ship's velocity relative to planet A

VA = planet A's velocity relative to sun S = sqrt(muS / RA)

Additionally, let sA be the radius of the SOI of planet A.

1) Start at circular orbit around A:

vA^2 = muA / rA

2) Burn dv1 up to velocity v1. We will determine dv1 and thus v1 later.

v1^2 = (vA + dv1)^2 = muA (2/rA - 1/aA)

or

1/aA = 2/rA - v1^2 / muA

3) Just before you leave planet A's SOI, you have slowed down to:

v2^2 = muA (2/sA - 1/aA)

v2^2 = v1^2 - muA (2/rA - 2/sA)

4) Once you exit A's SOI, (if you exit prograde) you gain planet A's orbital velocity VA:

v3 = v2 + VA

5) Now pick your velocity v3 outside the SOI so that you intercept your target half way around the sun:

So you are on an elliptical orbit around the sun with semi-major axis (RA + RB)/2, with velocity v3 at RA:

v3^2 = muS (2/RA - 2/(RA+RB))

6) You can solve for dv1 now, your transfer burn from planet A.

Find v3.

Solve for v2: v2 = v3 - VA

Solve for v1. v1 = sqrt(v2^2 + muA (2/rA - 2/sA)

dv1 = v1 - vA

7) Now do the reverse at planet B. Solve for the circularization burn at planet B.

Before entering the SOI at planet B you are going v4:

v4^2 = muS (2/RB - 2/(RA+RB))

8) Entering the SOI, you lose B's orbital velocity:

v5 = v4 - VB

9) Inside the SOI, find aB, the semi-major axis of your hyperbolic orbit of planet B:

v5^2 = muB (2/sB - 1/aB)

10) At planet B periapsis, use 1/aB to find your velocity:

v6^2 = muB (2/rB - 1/aB) = v5^2 + muB (2/rB - 2/sB)

11) Circularize at B periapsis:

vB^2 = muB / rB

dvB = v6 - vB

Edited May 1, 2014 by Yasmy

[Dkmdlb]

You ninja'd me. I was just about to post the very same thing....

.

.

.

.

.

.

.

.

.

I swear.

[Yasmy]

Ninja'd again, Dkmdlb!

I'll simplify the above so people can plug in numbers:

ÃŽâ€v₁ = sqrt( (μ_S/R_A) (sqrt(2R_B/(R_A+R_B)) - 1)² + μ_A (2/r_A - 2/s_A) ) - sqrt(μ_A/r_A)

ÃŽâ€v₂ = sqrt( (μ_S/R_B) (sqrt(2R_A/(R_A+R_B)) - 1)² + μ_B (2/r_B - 2/s_B) ) - sqrt(μ_B/r_B)

Pleasingly symmetric, as it should be. What it means is that you have to spend the Hohmann transfer ÃŽâ€v plus the planetary escape ÃŽâ€v burn minus your orbital velocity around the planet at each planet.

N.B.: I've edited the above to change sqrt(2/r_A - 2/s_A) to (2/r_A - 2/s_A). Likewise for the second equation.

Edited May 5, 2014 by Yasmy

[Yasmy]

The above applies to planets in circular or nearly circular orbits.

It also applies to transfers between two moons orbiting the same planet.

And, of course, if body B has an atmosphere, you can avoid paying most of ÃŽâ€v₂.

[Amagi82]

Ninja'd again, Dkmdlb!

I'll simplify the above so people can plug in numbers:

ÃŽâ€v₁ = sqrt( (μ_S/R_A) (sqrt(2R_B/(R_A+R_B)) - 1)² + μ_A sqrt(2/r_A - 2/s_A) ) - sqrt(μ_A/r_A)

ÃŽâ€v₂ = sqrt( (μ_S/R_B) (sqrt(2R_A/(R_A+R_B)) - 1)² + μ_B sqrt(2/r_B - 2/s_B) ) - sqrt(μ_B/r_B)

Pleasingly symmetric, as it should be. What it means is that you have to spend the Hohmann transfer ÃŽâ€v plus the planetary escape ÃŽâ€v burn minus your orbital velocity around the planet at each planet.

Plugging in values for a Kerbin to Moho trip, I'm getting the following (wrong) value for ÃŽâ€v₁.

http://www.wolframalpha.com/input/?i=%CE%94v1+%3D+sqrt%28+1.17e18%2F13599840256%29+%28sqrt%282*5263138304%2F%2813599840256%2B5263138304%29%29+-+1%29+%2B+3.53e12+sqrt%282%2F700000+-+2%2F84159286%29++-+sqrt%283.53e12%2F700000%29

Edited May 5, 2014 by Amagi82

[Yasmy]

1) See above where I've fixed a mistake. That second sqrt() shouldn't have been there. Sorry about that, and thank you for checking my math. My initial post on how to derive the delta-v was correct (if somewhat simple) though.

2) Careful how you enter into wolframalpha. You made a couple mistakes. Here ya go:

Idealized Kerbin -> Moho transfer burn.

And let me know if you have any more problems.

Edited May 5, 2014 by Yasmy

[LethalDose]

1) See above where I've fixed a mistake. That second sqrt() shouldn't have been there. Sorry about that, and thank you for checking my math. My initial post on how to derive the delta-v was correct (if somewhat simple) though.

2) Careful how you enter into wolframalpha. You made a couple mistakes. Here ya go:

Idealized Kerbin -> Moho transfer burn.

And let me know if you have any more problems.

What assumptions are being made for the "idealization" of the Moho burn? No inclination change and intercepting Moho at it's apoapsis?

Also, don't forget you'll need ~ 2 - 2.5 km/s dV to slow down at your interception.

[Yasmy]

We're not forgetting the insertion burn. That is ÃŽâ€v₂ above. We're just debugging the transfer burn calculation at the moment.

Just for kicks, the above formula for the insertion burn at Moho from Kerbin gives about 2.4k m/s for injection to a 25km circular Moho orbit.

Assumptions/idealizations/approximations:

1) circular orbits

2) zero inclination orbits

3) spacecraft exits/enters planet at the planet's orbital radius around Kerbol, rather than +/- the impact parameter (the distance between your SOI exit line and the planet prograde line)

So, YMMV, by a lot of dV. This was just a quick lesson on how to compute combined escape + Hohmann transfer + insertion.

Edited May 6, 2014 by Yasmy

[Amagi82]

Thank you very much, Yasmy. I've been working on an android app, and the equation is working properly now. For anyone else working with Java programs, the following may be helpful:

static double gravParameter[] = { 1.1723328e18, 1.6860938e11, 8.1717302e12, 8289449.8, 3.5316e12, 6.5138398e10, 1.7658e9, 3.0136321e11,

1.8568369e10, 2.1484489e10, 2.8252800e14, 1.962e12, 2.074815e11, 2.82528e12, 2.4868349e9, 7.2170208e8, 7.4410815e10 };

static double scaleHeight[] = { 0, 0, 7000, 0, 5000, 0, 0, 3000, 0, 0, 10000, 4000, 0, 0, 0, 0, 0 }; // atmospheric scale height

static double surfacePressure[] = { 0, 0, 5, 0, 1, 0, 0, 0.2, 0, 0, 15, 0.8, 0, 0, 0, 0, 0 }; // surface pressure in atmospheres

static double equatorialRadius[] = { 0, 250000, 700000, 13000, 600000, 200000, 60000, 320000, 130000, 138000, 6000000, 500000, 300000,

600000, 65000, 44000, 210000 };

static double siderealRotationVelocity[] = { 0, 1.2982, 54.636, 2.8909, 174.53, 9.0416, 9.3315, 30.688, 12.467, 24.916, 1047.2, 59.297,

17.789, 17.789, 0.75005, 0.30653, 1.2982 };

static double semiMajorAxis[] = { 0, 5263138304L, 9832684544L, 31500000, 13599840256L, 12000000, 47000000, 20726155264L, 3200000,

40839348203L, 68773560320L, 27184000, 43152000, 68500000, 128500000, 179890000, 90118820000L };

static double eccentricity[] = { 0, 0.2, 0.01, 0.55, 0, 0, 0, 0.05, 0.03, 0.14, 0.05, 0, 0, 0, 0.24, 0.17, 0.26 };

static double sphereOfInfluence[] = { 0, 9646663, 85109365, 126123, 84159286, 2429559, 2429559, 47921949, 1049599, 32832840,

2455985200L, 3723646, 2406401, 10856518, 1221061, 1042139, 119100000 };

static int minOrbit[] = { 0, 6817, 96708, 6400, 69078, 7061, 5725, 41447, 12725, 5670, 138155, 55262, 7976, 12695, 21758, 5590, 3874 };

static int highPoint[] = { 0, 6817, 7526, 6400, 6761, 7061, 5725, 8264, 12725, 5670, 0, 5600, 7976, 12695, 21758, 5590, 3874 };

// This method calculates the exact deltaV required to get to orbit around any body in KSP

public static double getToOrbit(int planet, double takeOffAltitude, double orbitClearance) {

double r1 = equatorialRadius[planet] + takeOffAltitude;

double r2 = equatorialRadius[planet] + minOrbit[planet] + orbitClearance;

double gravity = (gravParameter[planet]) / (r1 * r1); // force of gravity at a given altitude

double pressure = surfacePressure[planet] * Math.exp(-takeOffAltitude / scaleHeight[planet]); // pressure at a given altitude

double atmosphericDensity = 1.2230948554874 * pressure;

double terminalVelocity = Math.sqrt((1250 * gravParameter[planet]) / (r1 * r1 * atmosphericDensity));

double deltaV;

// Includes calculations for atmospheric losses only when atmosphere is present

if (surfacePressure[planet] != 0) {

deltaV = (Math.sqrt(gravParameter[planet] / r1) * Math.sqrt((2 * r2) / (r1 + r2))

+ (Math.sqrt(gravParameter[planet] / r2) * (1 - (Math.sqrt((2 * r1) / (r1 + r2))))) - siderealRotationVelocity[planet])

+ ((4 * gravity * scaleHeight[planet]) / terminalVelocity);

} else {

deltaV = Math.sqrt(gravParameter[planet] / r1) * Math.sqrt((2 * r2) / (r1 + r2))

+ (Math.sqrt(gravParameter[planet] / r2) * (1 - (Math.sqrt((2 * r1) / (r1 + r2))))) - siderealRotationVelocity[planet];

}

// Returns deltaV rounded to the nearest tenth.

return Math.round(deltaV * 10.0) / 10.0;

}

// Minimum delta V from orbit around planetStart necessary to intersect planetEnd's orbit, not including inclination changes

public static double getDeltaVInjectionBurn(int planetStart, int planetEnd, double altitudeStart) {

return Math.sqrt((gravParameter[0] / semiMajorAxis[planetStart])

* Math.pow(Math.sqrt((2 * semiMajorAxis[planetEnd] / (semiMajorAxis[planetStart] + semiMajorAxis[planetEnd]))) - 1, 2)

+ gravParameter[planetStart]

* ((2 / (equatorialRadius[planetStart] + altitudeStart)) - (2 / sphereOfInfluence[planetStart])))

- Math.sqrt(gravParameter[planetStart] / (equatorialRadius[planetStart] + altitudeStart));

}

// Minimum delta V necessary for capture and circularization at planetEnd, assuming no aerobraking

public static double getDeltaVInsertionBurn(int planetStart, int planetEnd, double altitudeEnd) {

return Math.sqrt((gravParameter[0] / semiMajorAxis[planetEnd])

* Math.pow(Math.sqrt((2 * semiMajorAxis[planetStart] / (semiMajorAxis[planetStart] + semiMajorAxis[planetEnd]))) - 1, 2)

+ gravParameter[planetEnd] * ((2 / (equatorialRadius[planetEnd] + altitudeEnd)) - (2 / sphereOfInfluence[planetEnd])))

- Math.sqrt(gravParameter[planetEnd] / (equatorialRadius[planetEnd] + altitudeEnd));

}