How to simulate an atmospheric straight-up ascent?

[AlexisBV]

Hello all,

I am trying to make a performance simulation spreadsheet for a very simple straight-up atmospheric launch, but for now I'm unable to match what I get in KSP, as my results are pessimistic by almost 4% for a sample case. I know 4% doesn't seem like much, but considering this should be a simple case, it doesn't given me much confidence to push into more complex stuff in its current state.

I've gone through it several times and even recoded it in Matlab, only to end up with the same results as the spreadsheet, which as I mentioned are different from flight test results. Can you help me find why?

Basically, the test case is a Mk1 pod, FL-T400 tank and a LV-N engine, plus 4 sensors (temp, baro, gravity, acceleration), for a total initial mass of 5320kg. The scenario is a straight-up launch from the KSC, starting at 77m altitude as indicated in game for this ship.

We then proceed with a full-throttle launch, with SAS on, pointing up (no lateral movement). In flight, we run out of fuel at about 16530m, then coast up to a max altitude of 20553m. No parachutes after that, sorry Jeb.

But from my spreadsheet, we should run out of fuel at about 15930m, and coast up to 19822m only.

Now I'd understand if my results were better than in KSP, as SAS isn't perfect and can induce some steering losses. But it's the other way around here… so I seem to be overestimating the losses in the spreadsheet, but I don't really see how.

Note I'm using an iteration step of 0.1 seconds, though it doesn't seem to make a huge difference - increasing the time resolution gives better results, but even at 1000Hz (0.001s time step), the apoapsis is only 19870m, still 682m short of KSP.

Anyway, here's the matlab code, let me know if you can find what's missing…

Thanks!

format long g

% --- Section 1 - Initializations ---
%
% In this section, environmental conditions and rocket parameters are set.
%

%Simulation time
endtime=500;
dt=0.1; %time step [s]

%Ship mass [kg]
Mass_Payload=820;  %Command Pod Mk1+4 sensors
Mass_FuelTankDry=250; %FL-T400
Mass_FuelTankFull=2250; %FL-T400
Mass_Engine=2250; %LV-N Atomic
FuelOnBoard=Mass_FuelTankFull-Mass_FuelTankDry;

%Thrust parameters
%LV-N Atomic
EngineMaxFg=60000; %[N]
EngineIspAtm=220;
EngineIspVac=800;
Throttle=1; %1 = 100%

%Environment
rho0=1.2230948554874;
Cd=0.2;
atmosScale=5000;
g0=9.81;
G=0.0000000000667384;
KerbinMass=5.2915793E+22;
KerbinRadius=600000;
GravParameter=G*KerbinMass;

%Variable initializations
altitude=zeros(1,endtime/dt);
airspeed=zeros(1,endtime/dt);
time=0:dt:endtime;
DryMass = Mass_Payload + Mass_Engine + Mass_FuelTankDry;
LostDeltaV = 0;
GroundElevation=77;

%Initial conditions
altitude(1)=77; %[m]
airspeed(1)=0; %[m/s]

%Potential DeltaV calculation
PotDeltaVatm = log((DryMass+FuelOnBoard)/DryMass) * EngineIspAtm * g0;
PotDeltaVvac = log((DryMass+FuelOnBoard)/DryMass) * EngineIspVac * g0;


% --- Section 2 - Simulation ---
%
% In this section, ascent performance is simulated.
%
% Three forces are calculated:
% - Gross thrust: fixed as long as there's fuel
% - Drag: in KSP, based on mass
% - Weight
%
% 
for i=1:(endtime/dt+1)
    Mass=Mass_Payload+Mass_FuelTankDry+Mass_Engine+FuelOnBoard;
    g=GravParameter/((KerbinRadius+altitude(i))^2);

    % --- Weight ---
    Weight=-Mass*g; %Weight always points down (negative)


    rho=exp(-altitude(i)/atmosScale) * rho0;

    % --- Gross thrust ---
    %Spec impulse interpolates linearly with rho from IspAtm to IspVac
    IspRatio=rho/rho0;
    Isp=IspRatio*EngineIspAtm+(1-IspRatio)*EngineIspVac;

    %Max thrust that can be generated with the fuel on board:
    Fg = min(Throttle*EngineMaxFg,FuelOnBoard*Isp*g0/dt);

    % --- Aerodynamic drag ---
    Drag=0.5*rho*(airspeed(i)^2)*Cd*Mass*0.008;
    if airspeed(i)>0
        Drag=-Drag; %Drag points down (negative) if airspeed is positive
    end

    accel=(Fg+Weight+Drag)/Mass;
    altitude(i+1)=altitude(i)+airspeed(i)*dt+(0.5*accel)*dt^2;
    if altitude(i+1)>=GroundElevation
        airspeed(i+1)=airspeed(i)+accel*dt;
    else
        airspeed(i+1)=0;
        altitude(i+1)=0;
    end

    %Fuel burn, cannot be more than the fuel on board
    FuelBurn=min(FuelOnBoard,Fg/(Isp*g0)*dt);

    %Potential DeltaV for this iteration (to calculate deltaV losses to
    %gravity and drag)
    if FuelBurn>0
        PotDeltaV = log((DryMass+FuelOnBoard)/ ...
            (DryMass+FuelOnBoard-FuelBurn))*Isp*g0;
        LostDeltaV = LostDeltaV + PotDeltaV - accel*dt;
    end

    %Fuel consumed
    FuelOnBoard = FuelOnBoard - FuelBurn;
end

%Trim last array element beyond time scope
altitude(i)=[];
airspeed(i)=[];

plot(time,altitude)
title(['Max alt: ',num2str(max(altitude))])

Edited December 9, 2013 by AlexisBV
Issue resolved

[Taki117]

I believe your spreadsheet is assuming the same specific impulse for your engine as you ascend. this isn't the case. As you get higher the atmosphere thins and the specific impulse of your engine increases towards it's vaccum value. I believe there is a point in atmo where your actual specific impule ans your vaccum specific impulse are the same.

[AlexisBV]

Hello Taki,

The spreadsheet/Matlab model does take into account the variable Isp (see line 78 in the code, under "Gross Thrust"), and I've checked it against the flight test. It matches nicely... so it's not that.

I've also double checked calculated gravity and baro pressure against flight test data (hence the sensors on my test rocket), and they match beautifully... so it has to be something else. But what?

Edited November 29, 2013 by AlexisBV

[Taki117]

AH! taking a closer look at your code your rocket equation is off. its supposed to be the natural log (ln), not the log.

http://en.m.wikipedia.org/wiki/Tsiolkovsky_rocket_equation

[AlexisBV]

Not that either... in Matlab, the natural log (base e) is "log()". The (base 10) log in Matlab would be spelled "log10()". In the Excel spreadsheet, I use LN().

Either way, that wouldn't change much as the delta-V calculation is not actually used to obtain acceleration, only as a provision for a planned feature in the code. It also serves as a good sanity-check for validating the acceleration code in vacuum.

In fact, that gives me an idea, I'll try this code/spreadsheet on a Mun scenario in KSP, if it matches, then I'll know the issue lies with the drag calculation (which isn't physically representative anyway)...

Thanks!

[cmoody]

AH! taking a closer look at your code your rocket equation is off. its supposed to be the natural log (ln), not the log.

http://en.m.wikipedia.org/wiki/Tsiolkovsky_rocket_equation

The log function in MatLab in the natural log. (Source: GEEN 1300 General Computing for Engineers)

I believe your spreadsheet is assuming the same specific impulse for your engine as you ascend. this isn't the case. As you get higher the atmosphere thins and the specific impulse of your engine increases towards it's vaccum value. I believe there is a point in atom where your actual specific impule and your vacuum specific impulse are the same.

This is your problem. Make your specific impulse a linear function (which it won't be in actuality) from sea level to vacuum with the endpoints the values given in the part info box. You'll see that your equations is much more accurate. Also note by Newton's equations that gravity is NOT constant with altitude. While this is not terribly important, you'll find that your numbers are still coming up short.

[mhoram]

The drag has to be measured by the relative airspeed.

And since the atmosphere moves with the planet, this must be accounted for in the calculation.

Have a look at the function get_drag_force in the module Satellite.pm of my ascent path optimizer.

These threads might be interesting for you:

http://forum.kerbalspaceprogram.com/threads/46194-I-need-someone-help-me-do-some-math-for-launch-optimization

http://forum.kerbalspaceprogram.com/threads/6664

[ferram4]

There's another issue; you're assuming that Isp changes linearly with atmospheric pressure, when it actually follows a Bezier Curve that has slopes of 0 at vacuum and one atmosphere. It's a small error, but it might account for the problem.

[AlexisBV]

Hmm, interesting point about the air moving with the planet. However I think that would only make things worse as relative airspeed (and drag) would then be higher.

Ferram4, the isp values calculated currently match with what I see in-game (it varies linearly with pressure)...

The program also calculates pressure and gravity matching what I see in KSP. I'm on my phone now, I'll post a table later when I get home.

Is there some way (add-on?) we can see what the current total drag force on the ship is?

[mhoram]

"MechJeb -> Vessel -> Atmospheric drag" shows the current drag in m/s^2

All references I found also stated that in KSP Isp scales linear with pressure.

Edit:

Another aspect is that a rocket also has a horizontal velocity. This increases the altitude of the rocket during each simulation step, because the planet is a spere. I think that this effect can be summarized as centrifugal force.

Edited November 29, 2013 by mhoram

[Specialist290]

There's another issue; you're assuming that Isp changes linearly with atmospheric pressure, when it actually follows a Bezier Curve that has slopes of 0 at vacuum and one atmosphere. It's a small error, but it might account for the problem.

Ferram4, the isp values calculated currently match with what I see in-game (it varies linearly with pressure)...

I'm willing to back AlexisBV up on this, given that it matches my own experiences as well (at least with rockets -- I'm not sure if jet engines work differently in that regard).

Is there some way (add-on?) we can see what the current total drag force on the ship is?

Kerbal Engineer Redux and MechJeb both output drag effects on your craft. I believe MJ outputs it as acceleration while KER outputs it as force, so that's something to be mindful of.

[numerobis]

The curves are linear if there are two points, and not linear if there are more than two points; so a turbojet gets more than its supposed maximum 2500 Isp. But for rockets, it's linear.

g0 is 9.82 in KSP. That probably explains almost everything.

Another effect: you are going 174.5 m/s orbit horizontal speed when you fly with 0 m/s surface horizontal speed. So gravity is a little tiny bit less than what you simulate.

Finally, I believe the gravitational parameter is more accurate than G * mass. Grab those values from the wiki: Kerbin is 3.5316e12.

Edited November 30, 2013 by numerobis

[Meithan]

The discrepancy is small that it's hard to assess the cause. I've been developing an aerobraking simulator mod, and drag is where I still get slight discrepancies with MechJeb values and game results. Since drag depends on the square of velocity, small errors in velocity will translate into noticeable errors in drag, which in turn alter the calculated velocity, and so on.

You might try switching to a higher-order numerical solver. If I understand your code correctly, you're using a first-order Euler method right now. It's not the best of methods. A simple option is the midpoint method. It's quite simple: first use the Euler method to estimate values for airspeed and altitude half a timestep ahead, and then use those "midpoint" values to obtain the new airspeed and altitude. You'll need to compute acceleration twice so you might want to define a function that does just that based on altitude and airspeed values.

[tavert]

Here's Mathematica's opinion, NDSolve uses a fancy adaptive method with very tight default tolerances (probably more accurate than KSP's own integration):

https://dl.dropboxusercontent.com/u/8244638/Vertical%20ascent.pdf (pardon the lack of plot labels)

Couple hundred meters off where you said burnout should happen, but matches apoapsis to within 25 meters. I'd put that small discrepancy down to the launchpad not being exactly on the equator.

Edited December 1, 2013 by tavert

[Camacha]

In fact, that gives me an idea, I'll try this code/spreadsheet on a Mun scenario in KSP, if it matches, then I'll know the issue lies with the drag calculation (which isn't physically representative anyway)...

You might want to use HyperEdit to manipulate the circumstances in KSP (change the atmosphere, gravity or else), as it helps to isolate contributing factors.

[AlexisBV]

I think I got it, based on tavert's example, it seems like it's the centripetal factor that was missing from my calculations. (I was only using thrust, gravity and drag, without correcting for the rotating planet reference). I added it to my spreadsheet, and while still not exact, it is way closer (about 0.1%) to what we see in game. I can live with that