I'm simulating a very simplified scenario of a missile falling to earth after being deployed at altitude. Currently I'm evaluating 1D only assuming the missile has already turned and is in a stable configuration downward. Running this case at 50,000 ft I get the following for missile velocity (blue) and terminal velocity (green)

resultant terminal and calculated velocity

Why am I exceeding terminal velocity?? I do see my peak velocity occurs when the terminal velocity and actual velocity are equal, however I would expect it to settle down to terminal velocity eventually? There are only two forces acting on this system, drag and weight, I've rationalized it in my head as the missile has more momentum than the increase in drag is able to counteract, but I'd like to have some verification of this or someone to point me to something more concrete.

Just to note, I've already verified my atmospheric calculations and cd calculations the problem isn't there. Have I misconceptualized something or is the concept of terminal velocity a simplification in how I've represented it? I've also attached the code. Thanks!

Edited to add:

$$\ F = m\alpha = -mg + \frac{1}{2} C_d \rho v^2A $$ $$\ \alpha = -g + \frac{1}{2m} C_d \rho v^2 A $$

where rho(altitude), Cd(mach,altitude), v(t), and I'm holding g constant

assuming a constant acceleration over a very small time step I use the basic kinematic formulations

$$\ v = v_0 + \alpha t $$ $$\ h = h_0 + vt + \frac{1}{2} \alpha t^2 $$

the standard formula for terminal velocity

$$\ F = drag - mg = 0 $$ $$\ v_t = \sqrt{\frac{2(mg)}{C_d \rho A}} $$

where once again Cd is a function of mach and altitude and rho is a function of altitude

The terminal velocity calculation takes into account the changing Cd value due to changes in mach number and with changes in density at altitude.

    releaseAltitude = sweep(iCase)

    %Initial Conditions
    [temp(1), rho(1), p(1), a(1)] = getAtmosphere(releaseAltitude);

    delT = .001; %s
    h(1)= releaseAltitude;
    vel(1)=-a(1)*0.9; %ft/s
    mach(1) = abs(vel(1)/a(1));

    Cd(1) = getCdAtMach(mach(1));
    dragPartial = Cd(1)*0.5*rho(1)*Area;

    dragAccel(1) = (dragPartial/mass)*vel(1)*vel(1);
    acel(1) = -g+dragAccel(1);
    dragForce(1) = dragAccel(1)*mass;
    termVel(1) = sqrt((2.0*weight)/(Cd(1)*rho(1)*Area));

    %Start Loop
    while (h(i) > 0)

        acel(i+1) = -g + dragAccel(i);
        vel(i+1) = vel(i)+acel(i+1)*delT;
        h(i+1) = h(i)+vel(i+1)*delT+.5*acel(i+1)*delT*delT;

        delh = h(i+1)-h(i);

        [temp(i+1), rho(i+1), p(i+1), a(i+1)] = getAtmosphere(h(i+1));

        mach(i+1) = abs(vel(i+1)/a(i+1));
        Cd(i+1) = getCdAtMach(mach(i+1));

        dragPartial = Cd(i+1)*0.5*rho(i+1)*Area;
        dragAccel(i+1) = (dragPartial/mass)*vel(i+1)*vel(i+1);
        dragForce(i+1) = dragAccel(i+1)*mass;

        termVel(i+1) = sqrt((2.0*weight)/(Cd(i+1)*rho(i+1)*Area));

        if(delh > 0)
        i = i+1;

more stuff here plotting etc .....

function [tempNew rhoNew pressureNew aNew] = getAtmosphere(altitude)

R = 1716;
g = 32.174;
gamma = 1.4;
tempSL = 59.018+459.67;    %See Intro to Flight 6th addition pg 112
tempLimit = -69.682+459.67;   %Rankine
rhoSL = 0.002377;   %slug/ft^3
pressureSL = 2116.2; %lb/ft^2

gradientAltitude = 36089; %ft end of gradient layer
run = tempSL - tempLimit; %Rankine
a1 = -run/gradientAltitude;

if(altitude < 36089)

    tempNew = altitude*a1 + tempSL;
    pressureNew = pressureSL*((tempNew/tempSL)^(-g/(a1*R)));
    rhoNew = rhoSL*((tempNew/tempSL)^-((g/(a1*R))+1));
    aNew = sqrt(gamma*R*tempNew);

elseif (altitude <= 82000)

    temp0 = tempLimit;                                          %temp at end of gradient layer remains constant
    pressure0 = pressureSL*((temp0/tempSL)^(-g/(a1*R)));        %pressure at end of gradient layer
    rho0 = rhoSL*((temp0/tempSL)^-((g/(a1*R))+1));              %density at end of gradient layer

    pressureNew = pressure0*exp(-(g/(R*temp0))*(altitude - gradientAltitude));
    rhoNew = rho0*exp(-(g/(R*temp0))*(altitude - gradientAltitude));
    tempNew = temp0;
    aNew = sqrt(gamma*R*tempNew);



function [Cd] = getCdAtMach(mach)

if(mach <= 1.2)

Cd = -12.0274*(mach^6) + 37.3284*(mach^5) - 39.9017*(mach^4) + 15.9694*(mach^3) + 0.3010*(mach^2) - 1.7736*mach + 0.6229;  %From FUN3D Analysis

elseif(mach > 1.2)

    Cd = 0.755;



  • $\begingroup$ @Eli Yes, I'll edit it now. $\endgroup$ Commented Oct 1, 2019 at 15:49
  • $\begingroup$ What is velocity of missle when it starts its downward motion $\endgroup$
    – Natsfan
    Commented Oct 1, 2019 at 17:38
  • $\begingroup$ @jmh about 1/2 V_terminal. Similar behavior occurs if the starting velocity = 0. $\endgroup$ Commented Oct 1, 2019 at 18:24
  • $\begingroup$ What specifically are your functions for $\rho$ and $C_d$? $\endgroup$ Commented Oct 2, 2019 at 3:24
  • $\begingroup$ @Aaron Stevens I've added those functions in. $\endgroup$ Commented Oct 2, 2019 at 15:50

2 Answers 2


This is just the behavior you'd expect. Terminal velocity is just that velocity where the drag and gravitational forces balance out, yielding zero acceleration--a body going faster than terminal velocity will slow down, and one going slower will speed up. As you've noted, at the one point where missile velocity is equal to terminal velocity, the slope of the missile velocity curve is flat: there is no acceleration. The missile would continue at exactly this velocity forever if terminal velocity didn't change, but it does: as the air gets thicker, drag goes up. Now, missile velocity is faster than terminal velocity, and the missile begins to slow down. It can't continue to match terminal velocity as it slows down, as then there would be no net force to cause it to decelerate. It must maintain a velocity higher than terminal velocity in order to do so.


The terminal velocity applies to a body dropped from rest or dropped with a velocity < terminal velocity. In your case, the missile is already going at mach speed and the Cd is less than if moving from rest. I haven't worked your problem but you run out of time to get below terminal velocity.

  • $\begingroup$ The initial velocity is less than terminal velocity, by almost half. Is this what you are getting at? ellipsix.net/blog/2010/04/calculating-terminal-speed.html $\endgroup$ Commented Oct 1, 2019 at 18:20
  • $\begingroup$ yes. Why are you using the drag coefficient for supersonic missile when missile is subsonic? $\endgroup$
    – Natsfan
    Commented Oct 1, 2019 at 18:52
  • 2
    $\begingroup$ The terminal velocity applies to a body dropped from rest or dropped with a velocity < terminal velocity. This isn't true. $\endgroup$ Commented Oct 2, 2019 at 3:20
  • 1
    $\begingroup$ @AaronStevens Yeah, I can't imagine skydiving would be as popular as it is if that were the case. $\endgroup$
    – JMac
    Commented Oct 2, 2019 at 19:28
  • 1
    $\begingroup$ @JMac I don't know, maybe it would be fun. Just make your cross-sectional area smaller to increase your terminal velocity. Then increase your cross-sectional area to lower the terminal velocity. Since you are now moving faster than your current terminal velocity, according to this answer you can now increase your speed as much as you want. Maybe I just want to go really fast :) $\endgroup$ Commented Oct 2, 2019 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.