Landing and Takeoff Delta-V vs TWR and specific impulse

[tavert]

I just did some calculations in Mathematica that are worth sharing. Especially at low thrust-to-weight-ratio (TWR), it turns out to be more efficient to land along a completely horizontal trajectory (pitching up to prevent gravity from reducing your altitude) than the common wisdom of doing a pure-retrograde suicide burn. The following numbers are best-case, assuming you're able to start from a circular prograde equatorial orbit at zero altitude before landing. I'm only considering the bodies without atmospheres here, since drag changes the story completely.

Javascript is disabled. View full album

The reason it takes less delta-V to land with a lower specific impulse (Isp) is that your mass decreases faster. There's no consideration here for engine or fuel mass, just overall TWR at the start of the landing burn. Higher Isp or better mass fraction will naturally give your craft a larger delta-V capacity to work from. The weight is relative to the body in question, at 0 altitude. I am considering the rotation speed of each planet/moon as the final condition, so the infinite-TWR asymptote is the orbital speed at zero altitude minus the sidereal rotation velocity at the equator and zero altitude.

Here are similar plots for takeoff, assuming the same type of constant-altitude trajectory:

Javascript is disabled. View full album

Hope people find this useful and/or interesting. Let me know if you have any questions, or if anything looks off.

Edit: I'm getting my celestial body data from this spreadsheet

Edit again: Sorry for the delay, I finished up the mathematical derivation document, you can find it here. I've also posted the code for landing and takeoff.

Edited July 21, 2013 by tavert
Adding mathematical derivation link

[tavert]

Here is an excellent video by Kosmo-not demonstrating this type of landing trajectory, to show you what I mean:

I should also note that I'm assuming a single stage is used throughout the entire landing, with a constant maximum thrust and the only change in mass being due to expending propellant. If you happen to stage in such a way that your TWR and Isp are exactly the same immediately before and immediately after the staging event, then the original numbers will still work.

Edited July 11, 2013 by tavert

[SchroedingersHat]

Firstly is that TWR standard gravity or gravity on the relevant body?

I just did some calculations in Mathematica that are worth sharing. Especially at low thrust-to-weight-ratio (TWR), it turns out to be more efficient to land along a completely horizontal trajectory (pitching up to prevent gravity from reducing your altitude) than the common wisdom of doing a pure-retrograde suicide burn.

Isn't a zero lift turn still going to be more efficient? I usually go for a zero lift almost-suicide burn. Set my periapsis at a bit past my landing spot then wait a bit until almost 0.5at^2= distance to landing spot and burn pure retrograde.

[tavert]

Firstly is that TWR standard gravity or gravity on the relevant body?

Isn't a zero lift turn still going to be more efficient? I usually go for a zero lift almost-suicide burn. Set my periapsis at a bit past my landing spot then wait a bit until almost 0.5at^2= distance to landing spot and burn pure retrograde.

Gravity on the relevant body. Counterintuitively, zero-lift (pointing retrograde the whole time... no lift for any trajectory in a vacuum) is actually worse. If you allow your altitude to drop, then gravity speeds you up and that ends up more important than if you accepted the steering loss and pitched up to fight gravity. This is more noticeable at low TWR. There are several threads from the past few months that came to this conclusion, both from numerical simulation and in-game tests. Let me dig those up and link to them...

Edit: a few recent (some more relevant than others) discussions to flip through if you're bored

http://forum.kerbalspaceprogram.com/showthread.php/34354-Land-on-Tylo-with-low-Thrust-to-Weight-Ratio

http://forum.kerbalspaceprogram.com/showthread.php/40139-How-to-calculate-optimal-descent-profiles

http://forum.kerbalspaceprogram.com/showthread.php/39721-Tylo-%28Tylo-Lander-Certification%29

http://forum.kerbalspaceprogram.com/showthread.php/39390-How-to-calculate-DeltaV-to-get-to-orbit-or-to-get-from-orbit-to-land

Edited July 18, 2013 by tavert

[SchroedingersHat]

Oh, I get it now. By burning pure horizontal you can accelerate down quicker than you would otherwise (even after steering losses) and spend less time (as gravity loss is a pure function of time and horizontal speed) accumulating gravity loss.

I'm guessing then that the ideal is going to be somewhere in between and would depend on the altitude which you decircularized from.

[tavert]

You don't burn horizontal, you burn at the angle necessary to keep your velocity vector horizontal. But yes the time you spend at full throttle is lower, hence less fuel burnt. There's some in-between trajectory that is the best you can achieve in practice that avoids crashing into the terrain.

Edited July 18, 2013 by tavert

[HoY]

What I'm taking away from this is after 2 the return is very minimal for the ammount of extra engine or lighter payload, so 2.0 is the magic number for me.

And the conversation your having,, that makes sense, and mechjeb's auto land function has been updated to land this way instead of doing suicide burns. I don't know when it changed (sometime after 2.0.2) but the suicide option is gone and it burns for quite a while keeping your vertical and horizontal speed as balanced as it can.

[tavert]

What I'm taking away from this is after 2 the return is very minimal for the ammount of extra engine or lighter payload, so 2.0 is the magic number for me.

You can combine these results with the TWR stats of the different engines and compare payload mass fraction for landing, there will be an optimal TWR and it's fairly low (especially for the LV-N). The slightly lower delta-V cost of landing with a higher TWR is not worth the extra dry weight of the additional engines beyond a point.

And the conversation your having,, that makes sense, and mechjeb's auto land function has been updated to land this way instead of doing suicide burns. I don't know when it changed (sometime after 2.0.2) but the suicide option is gone and it burns for quite a while keeping your vertical and horizontal speed as balanced as it can.

I don't think this is correct. I most often use the "land anywhere" button, which does a pure-retrograde near suicide burn (it's around 80-90% throttle for most of the landing in my experience). IIRC the "land at target" still tracks a velocity setpoint when landing without atmosphere, and points retrograde the whole time. It would be great to see this kind of constant-altitude landing implemented in MechJeb, but it's not there right now. A first step in this direction would be a Translatron-like feature that would control to a desired vertical speed setpoint using pitch instead of throttle.

[fibonatic]

I am looking forward to your mathematical derivation. Since I did try it myself, however I get a first-order nonlinear ordinary differential equation. But I am not that great at solving them XD:

dv / dt = sqrt((F * Isp / (m0 * Isp - t * F * g0))^2 - (mu / r^2 - v^2 / r)^2)

I might try to solve it numerically, since I tried using Wolframalpha and Matlab to solve it, but they weren't able to solve it.

[tavert]

I am looking forward to your mathematical derivation. Since I did try it myself, however I get a first-order nonlinear ordinary differential equation. But I am not that great at solving them XD:

dv / dt = sqrt((F * Isp / (m0 * Isp - t * F * g0))^2 - (mu / r^2 - v^2 / r)^2)

I might try to solve it numerically, since I tried using Wolframalpha and Matlab to solve it, but they weren't able to solve it.

That's the gist of it, yes (though you should be dividing F by g0, not multiplying). Rearranged a bit to be in terms of

TWR0 = F * r^2 / (m0 * mu), then solved numerically using NDSolve in Mathematica. In Matlab you'd use ode45 or similar, and it would be a few more lines of code overall.

I've got a nice LaTeX writeup in progress.

Edited July 19, 2013 by tavert

[fibonatic]

That's the gist of it, yes (though you should be dividing F by g0, not multiplying). Rearranged a bit to be in terms of

TWR0 = F * r^2 / (m0 * mu), then solved numerically using NDSolve in Mathematica. In Matlab you'd use ode45 or similar, and it would be a few more lines of code overall.

I've got a nice LaTeX writeup in progress.

I also realized this while making a Matlab script (I also used a wrong formula for the final delta-v calcualtion, also got another g0 wrong there).

But because the Matlab implementation of ode45 only has the option to choose a time span and not a stop criteria for when the velocity reaches orbital velocity (not 100% sure if events can be used for this) so I wrote my own solver. However I do not know the algorithm used by ode45, so I first used a RK2 method. Also tried a RK4 method, however this seemed to be less accurate than the RK2 compared to the ode45 solver. And also also made them more memory efficient, since you do not need to know the entire curve, but only the final value for time:

function dv = burnBudget(mu, r, v0, Isp, TWR, N)

g = 9.81;

vc = sqrt(mu / r);

f = @(t, u) sqrt((g * TWR * Isp / (Isp - t * TWR))^2 - (mu / r^2 - u^2 / r)^2);

h = Isp * (1 - exp(-vc / (g * Isp))) / (TWR * N);

v = v0;

k = 1;

while v < vc

fi = f(h * (k - 1), v);

vp = v;

v = v + h / 2 *(fi + f(h * k, v + h * fi));

k = k + 1;

end

dh = h * (vc - vp) / (v - vp);

for l = 1 : 5

dh = 2 * (vc - vp) / (f(h * k, vp) + f(h * k + dh, vc));

end

tf = h * (k - 1) + dh;

dv = g * Isp * log(Isp / (Isp - tf * TWR));

end

By the way I used another definition for: TWR = F / (m * g0).

[tavert]

FYI, I finished up the mathematical derivation document, you can find it here. I've also posted the code for landing and takeoff. These links are also at the end of the original post.

Finding the time for a final condition is exactly the kind of thing events are designed to do. You can also easily calculate the time at which mass goes to zero and use that (or some large fraction of it) as a worst-case maximum time, then interpolate on the result to find the final-condition crossing time. That's what I was doing at first, before I found that Mathematica 9 added events to NDSolve, and using them was less prone to convergence warnings than FindRoot on the InterpolatingFunction output of NDSolve.

Your TWR convention is only particularly useful for Kerbin, but I guess you can just do the scaling for different bodies at a higher level if you want. RK2's not very accurate, I prefer using a variable-step integrator for simple stuff like this. It will be a different story when I start incorporating dynamics into an optimization algorithm (it's nice to have a known finite dimension of your problem before you start), but this is just run-of-the-mill ODE integration here.

Edited July 21, 2013 by tavert

[fibonatic]

FYI, I finished up the mathematical derivation document, you can find it here. I've also posted the code for landing and takeoff. These links are also at the end of the original post.

Finding the time for a final condition is exactly the kind of thing events are designed to do. You can also easily calculate the time at which mass goes to zero and use that (or some large fraction of it) as a worst-case maximum time, then interpolate on the result to find the final-condition crossing time. That's what I was doing at first, before I found that Mathematica 9 added events to NDSolve, and using them was less prone to convergence warnings than FindRoot on the InterpolatingFunction output of NDSolve.

Your TWR convention is only particularly useful for Kerbin, but I guess you can just do the scaling for different bodies at a higher level if you want. RK2's not very accurate, I prefer using a variable-step integrator for simple stuff like this. It will be a different story when I start incorporating dynamics into an optimization algorithm (it's nice to have a known finite dimension of your problem before you start), but this is just run-of-the-mill ODE integration here.

Nice documents (however I find the mathematica code hard to read, partially because I have hardly any experience with it). My derivation was a bit shorter since I started in cylindrical coordinates, and just made one balancing equation in the radial direction, which yields an expression for the angle in which to burn and using sin(arccos(x))=sqrt(1-x^2) I was able to find the final equation for dv/dt.

I think I have found a better alternative for the numerical solution, since I noticed that both for take-off and landing the solution seems to be an exponential growth in the region of interest, which causes the solution to be unstable. Note that this is only true for landing when the Isp is not extremely high. Therefore it will always be hard to get good accuracy, since the error also experiences an exponential growth. But by changing the unknown from velocity to time you will be able to turn it into a stable numerical simulation, which needs significant less steps, and also has the advantage that you can end exactly at your desired velocity:

function dv = burnBudgetTakeoff(mu, r, v, Isp, TWR, N)

g = 9.81;

vc = sqrt(mu / r);

h = (vc - v) / N;

f = @(t, u) h / sqrt((g * TWR * Isp / (Isp - t * TWR))^2 - (mu / r^2 - u^2 / r)^2);

t = 0;

for k = 1 : N

fi = f(t, v);

v = v + h;

t = t + (fi + f(t + fi, v)) / 2;

end

dv = g * Isp * log(Isp / (Isp - t * TWR));

end

By the way I used F / (m * g0) for TWR, since I thought this would be easier to calculate for a user, since you are able to find the mass and thrust of your ship quite easily in the vanilla game.

[tavert]

Please ask if there are any specific things I could clarify about the Mathematica code, other than just unfamiliar syntax. In case it isn't clear, note that f[a_, b_] := a + b is the Mathematica syntax (or one of several, anyway) for defining f as a function of a and b. NDSolve is the numerical differential equation solver, which gets called as NDSolve[{x'[t] == ..., x[0] == xi, WhenEvent[x[t] == xf, "StopIntegration"]}, x, {t, 0, tmax}]. The rest of the code is derived values and plotting, pretty much.

There's no instability for any of typical cases, since none of the bodies in KSP are large enough to require extreme mass ratios to land on or takeoff from. There is a singularity as mass goes to zero, but that's way up at crazy delta-V's that we never need here (and can't reach with the mass ratio of KSP's fuel tanks anyway). Landing on Gilly with initial TWR < 1 leads to numerical trouble, since it's only feasible by slowing down until thrust counters the net difference between gravity and centripetal acceleration, and burning nearly vertically until you shed enough mass to have TWR >= 1.

Your approach of taking the reciprocal of the differential equation and integrating dt/dv instead of dv/dt does work, and you don't have to worry about events since you have fixed start and end velocities. I tried changing my code to do this in Mathematica, and integrating time as a function of velocity wound up over 2x slower in the end for takeoff, though it was a little bit faster for landing... odd. The results are close enough to be indistinguishable between the two methods (with default settings for NDSolve, anyway), and there's still no closed-form solution, as best I can tell.

Edited July 24, 2013 by tavert

[fibonatic]

When using the dt/dv method instead dv/dt, did you also compare to accuracy, for instance by comparing both results to a calculation using one of the methods in which you force the step sizes to be much smaller.

Because when I tested this with my RK2 solver, there was a significant difference. Because for the dt/dv method 1E2 steps would already yield a relative error of 1E-10, while for the dv/dt method I still only got a relative error of 1E-7 when using 1E7 steps.

This can be explained by the fact that the solution is either convergent or divergent, as can be seen in the image below, for which the explicit Euler method is used (so not very accurate solver).

For a convergent analytical solution the local error of the numerical solution does not get amplified by the following numerical steps, since the analytical solution of the same ODE, but with a small disturbance in the initial value will still convert to the same end value. However for a divergent analytical solution the local error of the numerical solution does get amplified by the following numerical steps, since a small disturbance in the initial value will eventually result in a totally different value.

PS: Your Mathematica code with be a little bit more readable if you would add some blank lines and maybe use some extra variable, for instance for the orbital and ground velocity. Because especially the first part is a very big block of code, which can be hard to read.

Edit: I just did a quick test for dy/dt = a * y (which is y(t) = y0 * exp(a * t)) in MATLAB and the reason I stated does not seems to hold very well, so I do not know why dt/dv seems to solve much more accurate for me than dv/dt.

Edited July 24, 2013 by fibonatic

[tavert]

NDSolve is using a quite advanced variable-step solver with error control and (IIRC) very tight default tolerances. It refines the step size wherever it needs to, whether that's at the start for a decaying solution or at the end for a growing solution. In your example above, remember that if the y-axis on the right is important, then it's the x-axis on the left that's important, and you'll be quite far off there too even though the y-axis on the left still converges to where it should eventually - it just takes a much different amount of time to do so, and that's what you actually cared about in this situation.

The differential equation residual from the default NDSolve solution is on the order 1e-6 m/s^2 for the v'(t) method, or 1e-7 s^2/m for the t'(v) method. Really only an issue for fixed-step methods if you use a step size that's already too large. And Mathematica has a whole set of solution methods that it automatically chooses from, I didn't ask which it used - it may have gone for an implicit method or something collocation-based, or something else entirely if it thought the problem was stiff or otherwise troublesome.

It does look nicer when viewed in Mathematica than it came out in the pdf, the line breaks were usually in better places and it does some syntax coloring. I couldn't get it to save to pdf both preserving screen appearance and in landscape, which I needed to avoid cutting off the figure legends, and I didn't want people to have to download Wolfram's notebook viewer program just to look at a few lines of code. I can have a go at formatting the Mathematica code to look even closer to LaTeX-style pretty-printing, just didn't do that to start with since it makes porting the code to and from other languages much more painful. And seeing code that uses equation-formatted fractions and Greek characters may confuse some people more than the equations themselves.

Edited July 24, 2013 by tavert

[Kerolyov]

Here is an excellent video by Kosmo-not demonstrating this type of landing trajectory, to show you what I mean:

I should also note that I'm assuming a single stage is used throughout the entire landing, with a constant maximum thrust and the only change in mass being due to expending propellant. If you happen to stage in such a way that your TWR and Isp are exactly the same immediately before and immediately after the staging event, then the original numbers will still work.

While this works a treat for low TWR landing on flat terrain I'm wondering how to apply this landing method to low gravity mountainous moons like Minmus especially landing on the plains, as impossible to get in low due to mountains in the way. Also the low gravity means its very likely that a lander will have high TWR unless you have Apollo style separate lander with tiny engines, which seems overkill as in this game separate lander/command modules usually together weigh more than an all in one.

So my latest Minmus lander has 3 man command pod and LV909 and has TWR of 6, so is there any benefit to doing anything other than a pure retro landing? If there is a better way than retro how do I assess what point to start landing burn and how do I estimate the angle to retro to start burning at?

[tavert]

Come in as horizontal as you can, don't put your apoapsis any higher than it needs to be to avoid terrain and don't put your periapsis too far below the surface.

Well you could save mass by using fewer/smaller engines - using a Rockomax 48-7S would give you a TWR of somewhere around 2.4-2.5, with only slightly worse vacuum Isp but saving 0.4 tons of dry mass. Your transfer injection burn may take a little longer though.

For higher TWR, the landing burn is so short that the difference between a pure-retrograde landing and a constant-altitude landing gets smaller and smaller. They're in fact identical in the limit of infinite TWR.

[Kerolyov]

Come in as horizontal as you can, don't put your apoapsis any higher than it needs to be to avoid terrain and don't put your periapsis too far below the surface.

Well you could save mass by using fewer/smaller engines - using a Rockomax 48-7S would give you a TWR of somewhere around 2.4-2.5, with only slightly worse vacuum Isp but saving 0.4 tons of dry mass. Your transfer injection burn may take a little longer though.

For higher TWR, the landing burn is so short that the difference between a pure-retrograde landing and a constant-altitude landing gets smaller and smaller. They're in fact identical in the limit of infinite TWR.

Thanks.

Problem is the plains have high mountains each end around the equator so you find yourself at 5km above the plains with only a short distance of maybe 30km on the plains before the next mountain so I think the best that can be done is to descend steeper and maintain a constant descent rate bleeding off horizontal speed as you descend with aim to reach zero altitude at a controlled rate just as horizontal velocity goes to zero.

Using Rockomax 48-7S does help but does lose some efficinecy with the TMI and plane change burns as your have to burn longer each side of the burn point.

[Kulebron]

So this graph practically means that higher Isp leads to spending more delta-v than with a lower Isp. (Of course it may still weigh less than with a low-isp engine.) Correct?

[tavert]

Yes. If you're landing on Tylo with a single stage, here's what you can expect the payload fraction to be,

meaning ((final mass) - (engine mass) - (propellant mass used in landing) * 9/8) / (initial mass):

If you want to land and take back off with the same stage, your payload fraction will be quite a bit lower (this is hard to do), and the best engine to use is LV-N at an under-1 starting TWR (it will go higher than 1 as you burn fuel):

Things change a bit as you stage, your TWR changes when you stage, and it can change either up or down depending how many engines and fuel tanks you drop.

[Superluminaut]

Hi nice data!

But, you clearly show a high twr is cheaper, but then show a video of a low twr landing? What?

Also whats the final mass of your model? It may be interesting to apply your data to various payloads, like landing fuel is x% of total mass.

[tavert]

That video is an extremely low TWR, yes. But it demonstrates the type of trajectory that I'm calculating for. The lower the TWR, the more difference the type of trajectory makes.

The ratio of final mass to initial mass can easily be calculated from delta-V and Isp using the rocket equation. The payload fraction data for Tylo I'm showing here is quite simple to derive from the initial delta-V numbers, you just need to plug in the engine thrust-to-weight data and the fuel tank mass fraction (9 units wet per 1 unit dry, for all but the smallest stock fuel tanks).

You need to consider not just fuel mass (and dry tank mass), but also engine mass. Higher TWR means you need less delta-V to land, but you need to carry more mass in engines.

Edited November 1, 2013 by tavert

[Superluminaut]

So essentially the perfect landing is producing a large directional explosion the instance you make contact with the ground, arresting all you surface relative velocity. Its a very Kerbal approach.

If I can figure out a good way of piloting such a landing, the mark 55 may actually have a legitimate purpose in my space program.

[mic_e]

So, how did your LaTeX writeup go?

Also, an other interesting graph would be payload/fuel mass fraction, for multiple-use ferrys.