# 1D falling object, and terminal velocity

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!

$$\ 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

run = tempSL - tempLimit; %Rankine

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

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

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