[1.1.2][kOS] Air drag calculation attempt

[schedar]

Hi,

So because there is no real answer so far on how to do a perfect ascend I did an attempt to calculate Air drag in kOS as a first step to finding method for perfect ascend. Basically I just assumed that:

where

 

So basically if we knew Potential and Real acceleration we could calculate AirDragAcceleration. We could see then if drag resistance is to high.

Below is a script in kOS that wil show thrust force (Ft), and accelerations Ar(real), Ap(potential) and Ad (drag). Be aware that it will calculate properly only when you are ascending in the atmosphere. Also turning your rocket will probably make readings inaccurate (I'm working on it). It's not going to be 100% correct but i think it's good start.

 

BTW. I'm using SSRSS mod that's why it's using Earth. You can use Kerbin if you don't have this mod.

 

Spoiler

CLEARSCREEN.
//Standard Earth Gravity Acceleration
set c_Earth_G0 to (earth:mass*constant:g/earth:radius^2).

until false {
//Calculate actual thrust
    LIST engines in l_Engines.
    set v_Ft to 0.
    FOR eng IN l_Engines {
        IF eng:IGNITION set v_Ft to v_Ft + eng:THRUST.
    }
    //Thrust force Ft in kN
    print "Ft:"+round(v_Ft,2) at(50,2).
    
//Earth gravitional acceleration with altitude
set v_Earth_Gh to c_Earth_G0*(earth:radius/(earth:radius+ship:altitude))^2.
//Earth gravitional force with altitude
set v_Earth_Fg to ship:mass*v_Earth_Gh.

//Potential acceleration without air drag on Earth
set v_Earth_Ap to (v_Ft-v_Earth_Fg)/ship:mass.
    print "Ap:"+round(v_Earth_Ap,2) at(50,5).
    
//Calculate dV and dT
set T_temp to sessiontime.
set V_temp_vector to SHIP:VELOCITY:SURFACE.
wait 0.1.
set delta_T to sessiontime - T_temp.
set delta_V_vector to SHIP:VELOCITY:SURFACE - V_temp_vector.
set delta_V_vector_mag to delta_V_vector:MAG.

//Real acceleration
//in direction of velocity
set v_Ar to delta_V_vector_mag/delta_T.
print "Ar:"+round(v_Ar,2) at(50,4).

//Airdrag Ad = Ap - Ar
Set v_Earth_Ad to v_Earth_Ap - v_Ar.
print "Ad:"+round(v_Earth_Ad,2) at(50,6).
    
}

 

Let me know what you think.

Edited June 3, 2016 by schedar
screenshot

[Adelaar]

A few questions... The first is; why don't you tap into the datastream of KSP itself, rather than to resort to calculations rond by a script. I'm no modder, but I do know how to visualize the drag by turning on the Aero GUI via Alt-F12. 

Secondly, why use drag for an optimal ascentprofile, drag isn't the only factor, the slingshot effect is too... 

Also, why not Mechjeb? I don't use it myself  but it does provide ascentprofile.

[stenole]

Your reasoning is correct, but I would try to keep things as vectors to make sure you're not adding errors into your calculations. Vectors will make your code cleaner because a lot of the variables that you extract in kOS will already be vectors and you have access to simple vector functions.

"set v_Earth_Ap to (v_Ft-v_Earth_Fg)/ship:mass. "

This line of code assumes that your thrust is always pointing opposite of gravity which means it will be incorrect if your rocket is not going straight up. If you had thrust and gravity as vectors instead of scalar numbers this would make it easier to apply newton's 2nd law of motion.

If you intend to use scalars and keep gravity going in the -1 direction, be aware that the angular speed (horizontal speed) will affect the gravity value because you are in a rotating reference frame; Think Centripetal force

[schedar]

On 6/6/2016 at 0:16 AM, Adelaar said:

A few questions... The first is; why don't you tap into the datastream of KSP itself, rather than to resort to calculations rond by a script. I'm no modder, but I do know how to visualize the drag by turning on the Aero GUI via Alt-F12. 

Secondly, why use drag for an optimal ascentprofile, drag isn't the only factor, the slingshot effect is too... 

Also, why not Mechjeb? I don't use it myself  but it does provide ascentprofile.

The idea is to use kOS instead of Mechjeb in order to have full control over vessel during ascent.

Drag is of course one of the elements I want to use as a parameter when ascending. I want to keep it optimal especially in lower parts of the atmosphere.

I'd love to extract the air drag numbers from datastream of KSP but I don't know how.

On 6/6/2016 at 1:30 AM, stenole said:

Your reasoning is correct, but I would try to keep things as vectors to make sure you're not adding errors into your calculations. Vectors will make your code cleaner because a lot of the variables that you extract in kOS will already be vectors and you have access to simple vector functions.

"set v_Earth_Ap to (v_Ft-v_Earth_Fg)/ship:mass. "

This line of code assumes that your thrust is always pointing opposite of gravity which means it will be incorrect if your rocket is not going straight up. If you had thrust and gravity as vectors instead of scalar numbers this would make it easier to apply newton's 2nd law of motion.

If you intend to use scalars and keep gravity going in the -1 direction, be aware that the angular speed (horizontal speed) will affect the gravity value because you are in a rotating reference frame; Think Centripetal force

Yes, you are very correct about this. I simply got lost in all those vectors, directions and rotations in kOS so I started simply with scalars and added some vectors where I could. Well it was my first attempt.

[surge]

Those equations will only work if you're going straight up, because gravity only acts in one direction whilst the thrust and drag can be in any.  Pretty useless for actual ascents.

A more effective equation for drag is 1/2 *q *Cd *A, where:

• q = dynamic pressure
• Cd = coefficient of drag
• A = surface area

(from here)

The problem is we don't know (A *Cd) anymore in 1.0 without doing experiments.  Cd is based on the 'drag boxes' you may have heard about, and no one knows how to calculate them yet AFAIK.

General dimensions of the ship aren't even available in kOS for some inexplicable reason. A is just frontal surface area, though, so width^2 in the editor.

Edited July 9, 2016 by surge
equation was wrong

[surge]

If you want ascent guidance try the following (xnode just does the maneuvre node).

The key seems to be vectored thrust pitch-over angle/altitude at 0-lift start, not drag... you can make nearly any rocket get into space with this by adjusting that value. Some of my more 'pointy' rockets use poang as low as 50, but if you're trying to launch a draggy, horrible pile of junk, you may need as much as 85-80. And even then it probably wont get into orbit - it's just not possible.

 

Spoiler

 

parameter talt, tinc, poang.
// Target Altitude, Target Inclination, Pitchover Angle

clearscreen.
print "== Launch Ascent Guidance System (LAGS) ==".

set talt to talt *1000. // target orbital perigee
set safetyspeed to 50. // speed at which safety turn ends.
set safestart to 20. // Altitude at which safety turn starts.
set pitchhold to 10. // halt zero-lift turn at this pitch
set debugmode to true.

set critq to 15. // maxQ Not implemented yet.
// fudge factor - how 'pointy' rocket is
set dragf to 0.4. // Not implemented yet

// ===END CONFIG OPTIONS===

// release controls early in case we get interrupted.
SET SHIP:CONTROL:PILOTMAINTHROTTLE TO 0.

lock LAc to 90 -vang(up:vector, ship:velocity:surface).
lock dynp to SHIP:DYNAMICPRESSURE *constant:AtmToKPa.

function printscreen {
// were going to use a 15 line screen total.
// Debug messages start on line 20
print " " at(0,1).
print "Target altitude: " +round(talt/1000, 1) +"km   " at(0,2).
print "Target inclination: " +round(tinc, 1) +"o   " at(25,2).
// Vector start, Vector end
// Launch mass, twr
print "Apogee: " +round(SHIP:OBT:APOAPSIS/1000,1) +"km     "at(0,6).
print "Perigee: " +round(SHIP:OBT:PERIAPSIS/1000,1) +"km    "at(25,6).
print "Time to apogee: " +round(ETA:APOAPSIS) +"s   "at(0,7).
print "Ecc: " +round(SHIP:OBT:ECCENTRICITY,2) +"   " at(0,8).
print "Inc: " +round(SHIP:OBT:INCLINATION,2)  +"   " at(25,8).
print "Dynamic pressure: " +round(dynp, 2) +"/" +round(critq, 2)
    +"kPa     " at(0,9).
print "Flight path angle: " +round(LAc, 1) +"   " at(0,10).
print " " at(0,11).
print " " at(0,12).
print " " at(0,13).
print " " at(0,14).
}

wait 0.5.
printscreen().
print "STATUS: calculating...                    " at(0,20).

function clamped {
// returns stage number of clamps any
    local stagecmp to STAGE:number.
    local highclamp to -1.
    
    if SHIP:STATUS <> "PRELAUNCH" and SHIP:STATUS <> "LANDED"
        return -1.

    local clamps to ship:partsnamed("launchClamp1").
    for c in clamps {
        if (c:stage > highclamp) set highclamp to c:stage.
    }
    return highclamp.
}

set coremotors to list().
set boostmotors to list().
function detect_motors {
// parse and sort motors in current stage into
// globals coremotors and boostmotors
    
    local el to list().
    list engines in el.
    
    set coremotors to list().
    set boostmotors to list().
    
    local stagehigh to STAGE:number.
    local stagelow to clamped().
    if (SHIP:STATUS = "PRELAUNCH") and (stagelow = -1) {
        set stagelow to stagehigh -1.
    } else if (stagelow = -1)
        set stagelow to stagehigh.

    for e in el {
        // we have to detect if stage above has launch
        // clamps AND if we are in launch ready state
        if (e:stage >= stagelow) and (e:stage <= stagehigh)
        {
            // Check if the engine is
            // on centre stack.
            if (e:tag = "boostMotor") {

                boostmotors:add(e).
            } else {
                // everything else should be cores.
                coremotors:add(e).
            }
        }
    }
    return coremotors:LENGTH +boostmotors:LENGTH.
}

function setcorethrottle {
    parameter lvl.
    for e in coremotors {
        set e:THRUSTLIMIT to lvl *100.
    }
}

function getcorethrottle {
    if (coremotors:length = 0) return 0.
    return coremotors[0]:THRUSTLIMIT /100.
}

function getctwr {
    // returns current TWR.
    // *sigh* :MAXTHRUST doesnt work unless the engine is
    // staged. This SHOULD return the TWR upon launch, so
    // we only have to calculate it once.
    local gf to (SHIP:BODY:MU
        /(SHIP:BODY:RADIUS +SHIP:ALTITUDE)^2).
    local thrust to 0.

    list engines in el.
    for e in el
        if e:IGNITION set thrust to thrust +e:THRUST.

//    return SHIP:AVAILABLETHRUST /(SHIP:MASS*gf).
    return thrust /(SHIP:MASS*gf).
    
}

function inc2hdg {
// does the wierd maths when converting from an orbital inclination
// to a heading to fly for it.
    parameter inc.
    set rt to 180 -inc -90.
    if (rt < 0) set rt to 180 -rt.
    return rt.
}
set thdg to inc2hdg(tinc).

set east to up *R(0, -90, 0).
function getvhdg {
    local r1 to vcrs(up:vector,
        ship:velocity:orbit):normalized.
    local vhdg to vang(r1, east:vector).
    if (r1:y) > 0.0 set vhdg to 360 -vhdg.
    return vhdg.
}

set laststage to STAGE:NUMBER.
function stagectl {
    if not STAGE:READY return false.
    if THROTTLE <= 0 return false.

    // detect_motors() takes long and only needs to be
    // called once after staging has finished
    if (laststage <> STAGE:NUMBER) {
        detect_motors().
        set laststage to STAGE:NUMBER.
        return false.
    }
    
    // also check for no motors
    if (coremotors:length = 0) and (boostmotors:length = 0){
        lock throttle to 1.
        stage.
        return true.
    }

    for e in boostmotors {
        if (e:maxthrust <= 0) {
            lock throttle to 1.
            stage.
            return true.
        }
    }

    if (coremotors:length > 0) {
        // check that boost motors are gone.
        // but still need to do the fuel check thing.
        if (boostmotors:length = 0) setcorethrottle(1).

        local corethrust to 0.
        for e in coremotors {
            set corethrust to corethrust +e:maxthrust.
        }
        if (corethrust = 0) {
            stage.
            return true.
        }
    }
    return false.
}

function hdgerr {
    // Calculates the heading offset for the given orbital
    // inclination.
    parameter targethdg.
    
    // ffs. SHIP:OBT:INCLINATION doesnt do negatives?
    set shipinc to SHIP:OBT:INCLINATION.
    local shiphdg to getvhdg().
    if (shiphdg > 90 and shiphdg < 270) set shipinc to -shipinc.
    return targethdg -inc2hdg(shipinc).
}

detect_motors().
if boostmotors:LENGTH = 0 {
    print "WARNING NO BOOST MOTORS DETECTED!" at (0,14).
//    exit.
}

//
// ===GUIDANCE CALCULATION===
//

set poalt to SHIP:BODY:ATM:HEIGHT/7.

// FIXME: This needs to settle faster.
SET thpid to PIDLOOP(0.1, 0.8, 0.15, 0, 1).
set thpid:setpoint to 0.0.

SET pchpid to PIDLOOP(0.5, 0.05, 0.5, 0, 20).
set pchpid:setpoint to 15.

//set STEERINGMANAGER:PITCHPID:KD to 0.5.
//set STEERINGMANAGER:YAWPID:KD to 0.5.

FROM {local countdown is 5.}
UNTIL countdown = 0 STEP {SET countdown to countdown - 1.} DO {
    PRINT "STATUS: pre-launch... T-"
    + countdown +"s            " at(0,20).
    WAIT 1.
}

print "STATUS: Launch...                    " at(0,20).
lock throttle to 1.
lock steering to up.
// Since engine:MAXTHRUST is broken, we have to actually fire
// up the engines (possibly launch if no clamps) to get a value.
set launchtwr to 0.
set sa to SHIP:ALTITUDE -SHIP:GEOPOSITION:TERRAINHEIGHT.
until SHIP:ALTITUDE -SHIP:GEOPOSITION:TERRAINHEIGHT-sa > safestart {

    set launchtwr to getctwr().

    local dTWR to getctwr() -1.3.
    setcorethrottle(thpid:update(TIME:SECONDS*10, dTWR)).
    
    print "Launch TWR: " +round(launchtwr, 2) +"     " at(25,3).
    print "Launch Mass: " +round(SHIP:MASS, 2) +"    " at(0,3).
    stagectl().
    
    printscreen().
    wait 0.1.
}
unset sa.
//set poang to 90-(-100*(1/dragf*0.2)*(launchtwr -1.8)).
print "Zero-lift start: " +round(poalt) +"m   " at(0,4).
print "Vectoring angle: " +round(poang, 1) +"o    "at(25,4).
unset startalt.

print "STATUS: safety turn...                    " at(0,20).
lock steering to heading(thdg, 85).
set rcs to true.
until (SHIP:VELOCITY:SURFACE:MAG >= safetyspeed) {
    stagectl().

    printscreen().
    wait 0.5.
}

print "STATUS: vector guidance...                    " at(0,20).
lock pitcherr to SHIP:ALTITUDE *(90-poang)/poalt.
lock steering to heading(thdg, 90 -pitcherr).
until (SHIP:ALTITUDE > poalt)
{
    stagectl().

    printscreen().
    wait 0.5.
}

print "STATUS: zero-lift turn...                    " at(0,20).
lock steering to heading(thdg, LAc).
// FIXME: do throttle control
lock steering to heading(thdg +hdgerr(thdg), LAc).
until (LAc <= pitchhold) or (SHIP:APOAPSIS >= talt) {
    stagectl().

    printscreen().
    wait 0.5.
}

print "STATUS: pitch hold...                    " at(0,20).
lock steering to heading(thdg +hdgerr(thdg), pitchhold).
until (SHIP:ALTITUDE >= SHIP:BODY:ATM:HEIGHT) {
    stagectl().
    lock throttle to thpid:update(TIME:SECONDS,
        SHIP:APOAPSIS -talt).

    local pitchang to pchpid:update(TIME:SECONDS*10, LAc).
    lock steering to heading(thdg +hdgerr(thdg), pitchang).

    printscreen().
    wait 0.5.
}
lock throttle to 0.
unlock steering.
unlock LAc.
unlock dynp.
RCS off.

print "STATUS: Orbital Insertion...                 " at(0,20).
set apR to APOAPSIS +BODY:RADIUS.
set dV to sqrt(BODY:MU /apR)
    -sqrt(BODY:MU* (2/apR -1/OBT:SEMIMAJORAXIS)).
set circNode to node(TIME:SECONDS +ETA:APOAPSIS, 0, 0, dV).
add circNode.
// we could check the node is correct by comparing
// circNode:ORBIT:APOAPSIS and circNode:ORBIT:PERIAPSIS
// to make sure they are a) nearly equal, and b) close to talt

run xnode(1, 0).
print "STATUS: Program ended.                " at(0,20).

 

 

 

Edited July 9, 2016 by surge