Fuel consumption as a function of atmospheric pressure

[stan-]

Hello,

I noticed that the fuel consumption of my engines drops as my ship rises higher, but the wiki does not state an explicit equation for the relation between consumption and height. I have some empirical evidence to believe that there is some kind of sigmoidal relation, i.e. consumption = a * (1 - S(k * height)) + b, where a and b are engine-specific parameters (derived from what is displayed in the engine description for fuel consumption), and k is a parameter of the game's physics engine.

My question is: am I even remotely right?

Thanks,

- Stan

[PakledHostage]

This question might be better asked in the Gameplay Questions and Tutorials sub-forum, but I'll take a stab at it. I am sure someone will correct whatever I get wrong (and what I get right, for that matter...)

Ans: The Isp of engines in the game varies with atmospheric pressure, but the implementation is bass ackwards. In order to achieve the variable Isp, they've held engine thrust constant over the full range of atmospheric pressures and varied fuel burn instead. In reality, you would expect engine thrust to vary with back pressure (i.e. atmospheric pressure) acting on the nozzle while fuel flow rate would be roughly constant for any given throttle setting. I recall seeing a post once from one of the devs that they were going to fix this, but I don't think it has been done yet.

[K^2]

The relationship is piecewise-linear. The way it works in code is via the following lines in the engine .cfg file.

atmosphereCurve
 {
  	 key = 0 370
 	 key = 1 320
 }

And the game just uses first order interpolation, which means I_SP is 370 at pressure of 0 bar (or bellow, but that's impossible), 320 at pressure of 1 bar or above, and for any pressure x between 0 and 1 bar, it is x*320 + (1-x)*370.

This does, of course, have the general features of a sigmoid, but it's not a smooth curve. Similar code is used for I_SP and thrust of game's jet engines, but these use more key points. Interpolation is still linear, however.

Naturally, knowing thrust and I_SP, you can compute fuel consumption.

Edit: Small correction. For jet engines, thrust varies with velocity, not pressure. But again, using the same key-value interpolation. I_SP varies with pressure, same as rocket engines.

Edited January 28, 2014 by K^2

[stan-]

Thank you both for the answers! What K^2 said was what I was after. Will add that to my computations.

Thanks!

[tavert]

Similar code is used for I_SP and thrust of game's jet engines, but these use more key points. Interpolation is still linear, however.

Interpolation of Isp for jets is not actually linear. You'll notice for the turbojet for example you have

	atmosphereCurve
 	{
		key = 0 1200
		key = 0.3 2500
		key = 1 800
 	}

but the reported Isp in the right-click menu peaks at slightly higher than 2500, around 2525 or so. When there are more than 2 points in the curve it's some sort of spline or polynomial, not piecewise linear.

Also I only recently became aware of the fact that for jets in KSP, the reported specific impulse in the right-click menu is not accurate. The effective Isp, equal to thrust / (g0 * mdot), is currently scaled by the velocityCurve. I suspect that's a bug.

[KerbMav]

Maybe this mod - or the reading of its source code - can help you too?

http://forum.kerbalspaceprogram.com/threads/52882-0-23-Kerbal-Isp-Difficulty-Scaler-v1-3-2-12-17

[K^2]

Interpolation of Isp for jets is not actually linear. You'll notice for the turbojet for example you have

	atmosphereCurve
 	{
		key = 0 1200
		key = 0.3 2500
		key = 1 800
 	}

but the reported Isp in the right-click menu peaks at slightly higher than 2500, around 2525 or so. When there are more than 2 points in the curve it's some sort of spline or polynomial, not piecewise linear.

Also I only recently became aware of the fact that for jets in KSP, the reported specific impulse in the right-click menu is not accurate. The effective Isp, equal to thrust / (g0 * mdot), is currently scaled by the velocityCurve. I suspect that's a bug.

Odd. Polynomial interpolation through these three points peaks at 0.47 bar with value 2,696. If you are saying that the maximum is above 2,500, but nowhere near as high as this, that excludes this particular possibility. On the other hand, something like Bezier can't exceed maximum key value. So that's me out of ideas. I guess I'll go and dig through the actual code now.

Edit: Game uses Unity's AnimationCurve object for interpolation. It does use ClampForever option for end points. So the value at 1 bar is valid for any pressure above 1 bar. Apparently, AnimationCurve can look after tangents, so I tried using interpolating polynomial which has zero slope at the end points as well, but that gave me a maximum value that's even higher - 2,970. And I haven't been able to find more detailed documentation on AnimationCurve to find out exactly what it uses for interpolation.

P.S. Here is the actual code from the ModuleEngines.

private double RequiredPropellantMass(float throttleAmount)
{
    this.realIsp = this.atmosphereCurve.Evaluate((float)base.vessel.staticPressure);
    float single = this.realIsp * this.g;
    this.requestedThrust = Mathf.Lerp(this.minThrust, this.maxThrust, throttleAmount);
    return (double)(this.requestedThrust / single * TimeWarp.deltaTime);
}

Edited January 29, 2014 by K^2

[Idobox]

K^2, that formula can be simplified as 370 - 50x, it's an linear function, not a sigmoid.

It also means that for high enough pressures, the ISP could become negative (meaning negative fuel consumption?). I'm surprised KSP uses such a simple function.

Do you have any idea what it looks like for real engines?

[K^2]

K^2, that formula can be simplified as 370 - 50x, it's an linear function, not a sigmoid.

It also means that for high enough pressures, the ISP could become negative (meaning negative fuel consumption?). I'm surprised KSP uses such a simple function.

I've said that it was piecewise-linear. It only follows the above formula between 0 and 1. Even that might not be exactly right. See bellow. Outside of that range, however, behavior is very simple. It's a constant at whatever value the end point had. So for any pressure above 1 bar, it's going to be 320s.

The actual behavior of the AnimationCurve is more complicated. It seems to be piecewise-cubic, with the 4 fitting parameters being the 2 values and 2 tangents at either end. The problem is that you can define the curve without specifying the tangents. And I have not yet figured out what Unity substitutes for the tangent values in that case.

The cool thing is that even though none of the modules make use of it, KSP does allow you to specify a curve with custom tangents! So instead of using "key pressure value", you can use "key pressure value tangent_in tangent_out". This allows you to define almost any smooth curve you like for atmosphereCurve and the velocityCurve.

Do you have any idea what it looks like for real engines?

I don't, but I might know where to look this up, if you are interested. The thing, though, is that it very much depends on the bell. The actual bell shape depends on a few parameters, and is basically just designed to match the engine. But the point at which the curve of the bell is truncated depends on design operating pressure. The ideal nozzle for a rocket engine narrows at first to make the exhaust speed up. However, once the flow reaches speed of sound, the nozzle needs to start expanding. Hence the bell. The bell then keeps on expanding until either the pressure in the exhaust flow matches external pressure, or it no longer becomes worth it to increase the mass of the engine for the extra impulse you are going to get.

Consequently, an engine designed for operating in atmosphere is going to have a short bell. It will perform a little better in atmosphere than an engine designed for vacuum, but the main benefit is a much lighter engine due to much shorter bell. The drawback is that it will gain almost nothing when operating at lower pressures than the pressure that it's designed for.

Engine designed for vacuum is going to have a very long bell, since there is no external pressure that dictates a cutoff point. It's just as long as they can make it without making an engine too heavy. This engine will have its peak I_SP in vacuum, and will suffer significant penalties when operating in atmosphere.

Edit: I think I worked out what the algorithm for AnimationCurve actually is. I've used the key-value pairs for the jet engine posted earlier in the thread.

atmosphereCurve
{
	key = 0 1200
	key = 0.3 2500
	key = 1 800
}

And I replaced these with key-value-slope triplets.

{{0, 1200, 4333.33}, {0.3, 2500, 952.381}, {1, 800, -2428.57}}

The way I got these is by taking slopes of lines connecting a key to its neighbors and averaging these. For end points, I've averaged them with themselves. Once I had these, I used 3rd order polynomial to connect the keys, matching slopes at both end points. And this is the curve that I got.

Why do I think this might be right? Because it peaks at 0.35 bar with a value of 2,524s. Which is very close to what tavert indicated, and I'm strongly inclined to believe him.

Trying to "average" the end-points as well, assuming slope of 0 beyond end-points, results in a slightly different keys with a lower maximum. So it looks like end-points aren't averaged. That means that if the function is defined by just two keys, it will, indeed, be linear between the keys.

Edited January 30, 2014 by K^2

[tavert]

The cool thing is that even though none of the modules make use of it, KSP does allow you to specify a curve with custom tangents! So instead of using "key pressure value", you can use "key pressure value tangent_in tangent_out". This allows you to define almost any smooth curve you like for atmosphereCurve and the velocityCurve.

The velocityCurve's actually do specify tangents. Ferram also puts some nonzero tangents in some of his engine curves in FAR.

Edit: I think I worked out what the algorithm for AnimationCurve actually is.

Good catch! That looks like the most accurate explanation I've seen yet. I'll run this against numerobis' data here https://github.com/numerobis/KSP-scripts/commit/0b3fb888bb5ec7b10eaadccc6216aa30ad2c604e and see how close it all lines up. He's been wanting to figure this out for a while, if he has already found the solution then his code doesn't reflect it yet.

I was playing around yesterday with Catmull-Rom splines, but that's not quite what this is.

Edited January 30, 2014 by tavert

[K^2]

Perfect fit. Magenta points are from the page you linked to, one with numerobis' data. Blue lines are my piecewise-cubic fits.

Jet Engine:

Turbojet:

[tavert]

You beat me to it, I did the thrust curves first since they have tangents defined. Perfect fits there too:

Javascript is disabled. View full album

Mathematica code:

AnimationCurve[Keys_, Values_, TangentsIn_, TangentsOut_] :=
 Function[x, Piecewise[Join[
    {{Values[[Ordering[Keys][[1]]]], x <= Min[Keys]},
     {Values[[Ordering[Keys][[-1]]]], x >= Max[Keys]}},
    Table[{InterpolatingPolynomial[
       {{Keys[[i]], Values[[i]], TangentsOut[[i]]},
        {Keys[[i + 1]], Values[[i + 1]], TangentsIn[[i + 1]]}}, x],
      Min[Keys[[i ;; i + 1]]] <= x <= Max[Keys[[i ;; i + 1]]]},
     {i, 1, Length[Keys] - 1}]]]]

DefaultCurve[Keys_, Values_] :=
 With[{TangentFirst = (Values[[2]] - Values[[1]])/
     (Keys[[2]] - Keys[[1]]),
   TangentLast = (Values[[-1]] - Values[[-2]])/
     (Keys[[-1]] - Keys[[-2]]),
   TangentMiddle = (Values[[2 ;; -1]] - Values[[1 ;; -2]])/
     (Keys[[2 ;; -1]] - Keys[[1 ;; -2]])},
  AnimationCurve[Keys, Values,
   Join[{TangentFirst}, (TangentMiddle[[1 ;; -2]] +
       TangentMiddle[[2 ;; -1]])/2, {TangentLast}],
   Join[{TangentFirst}, (TangentMiddle[[1 ;; -2]] +
       TangentMiddle[[2 ;; -1]])/2, {TangentLast}]]]

Edited January 31, 2014 by tavert

[K^2]

There is one less mystery in the world. I almost feel sad.

[tavert]

There is one less mystery in the world. I almost feel sad.

Would you feel any differently if a little dialog popped up saying "You have earned 5 science points"?

[Idobox]

I've said that it was piecewise-linear

Sorry, I misunderstood your post.

Thanks for the amazing info

[numerobis]

Ah, great. I'll have to code this up now...

[antbin]

Great work Talvert. Perfectly shows how the RAPIER is dominated by the turbojet under all conditions.

Charts like these should give some insight into game design - what performance maps for different components would make for fun gameplay?

IMHO a final basic jet / turbojet / ramjet (?) / RAPIER game balance design would be most fun if each engine had a different peak thrust speed. That would encourage building craft that had a mix of engines, and switching from one to the other on the way to orbit.

[numerobis]

K2: How do you compute the tangent in and tangent out exactly? You're getting one number per point, whereas I would expect two -- the in tangent and the out tangent.

I updated my scripts to get the thrust curves exactly right.

[K^2]

If the in/out tangents for key are specified, you just use the values that are there. If tangents are not specified, you compute what the in and out tangents would be in pieceweise-linear case. But then you set both of them to average, id est, t_in = t_out = (t_in + t_out)/2. That gives you a smooth curve instead of a piecewise-linear one. Of course, for end points, you only have t_in or t_out, so you just use that value, rather than averaging it with zero or anything else.

I've written a C++ class that handles animation curves. So if you find it easier to reference code, feel free to take a look at it. I wasn't trying to code for clarity, though.

Oh, and if anyone wants to use the whole thing or any parts in their work, consider this an unconditional permission to do so.

curve.h

curve.cpp

Be warned that I have not tested DropKey methods, and there is no clean destructor. (Keys will be lost in memory if you delete parent curve.)

[numerobis]

t_in = t_out is the part I was derping on. Thanks!

[tavert]

One thing I'm not certain on since none of our tests had nonzero tangents provided in the curves, was whether column 3 is the in tangent and column 4 is the out tangent, or the other way around. Until someone points to where the part.cfg data gets turned into an AnimationCurve in the code, or collects some conclusive test data for an engine with some nonzero tangents, I think we're all still guessing on this detail. Doesn't yet matter for any of the stock engines, but worth noting.

Also K^2, where/how did you get that ModuleEngines code snippet from?

[K^2]

void FloatCurve::Add(float time, float value, float inTangent, float outTangent)

...

this.Add(float.Parse(strArrays[0]), float.Parse(strArrays[1]), float.Parse(strArrays[2]), float.Parse(strArrays[3]));

So the order in the config file is key, value, inTangent, outTangent. The in tangent is on the left, while the out tangent is on the right of the key.

[Starwaster]

This was a very informative thread as I'm doing some curve work but I'm still not clear on the tangents.

I thought I had a method worked out for calculating them but when I tried them out using the Mathematica code provided in this thread I got a nasty collection of spikes.

What I did was this basically...

deltaY = P2_y - P1_y

deltaX = P2_x - P1_x

atan2(deltaY, deltaX)

What I've learned from this experience is that apparently I don't know a tangent from a poptart. (actually that's not true... poptarts are frosted)

[K^2]

atan2 will give you the angle of the tangent. You don't need that. You just need the actual slope at each point. Then you construct a 3rd degree polynomial that matches the two end points and the two slopes at either end of the interval between the two keys.

[Starwaster]

atan2 will give you the angle of the tangent. You don't need that. You just need the actual slope at each point. Then you construct a 3rd degree polynomial that matches the two end points and the two slopes at either end of the interval between the two keys.

Thanks, I'll try that. Big question though, will that leave the defined values intact when the curve is evaluated? For instance, if I'm defining a pressureCurve to match Earth and it says that at 5.5 km (key) the pressure will be 0.48 (value) will I still get that value when it's evaluated for 5.5 km? Or will there be unwanted smoothing?