1D falling object, and terminal velocity

Question

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)

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;

    t(1)=0;
    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
    i=1;
    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;

        t(i+1)=t(i)+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)
            break
        end
        i = i+1;
    end

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);

end

end

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;

end

end

Answer 1

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.

Answer 2

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.