Fald i atmosfære

\begin{equation} \vec{F}_g = m \vec{g} \end{equation} \begin{equation} \vec{F}_d = \frac{1}{2} \rho \vec{v}^2 C_d A \end{equation}

Analytisk løsning

\begin{equation} v(t) = \sqrt{ \frac{2mg}{\rho A C_d} } \tanh \left(t \sqrt{\frac{g \rho C_d A}{2 m}} \right) \end{equation}
clc;
dt = 0.0001;
g = 9.8; % Tyngde-acceleration m/s^2
a = [0.0 -g 0.0]; % m/s^2
rho = 1; % () Density of fluid

kugle = struct(); % Det faldende objekt
kugle.position = [0, 100 ,0]; % (m) begyndelsesposition
kugle.hastighed = [0, 0, 0]; % (m/s) begyndelseshastighed
kugle.A = 10; % (m^2) Reference areal (tværsnittet af kuglen)
kugle.C_d = 0.47; % Drag coefficient (se http://en.wikipedia.org/wiki/Drag_coefficient)
kugle.masse = 10; %(kg) massen af kuglen

figure(1) % position
clf;
% plot
semilogx(0, kugle.position(2),'o');
xlabel('t')
ylabel('y(t)')
hold on

figure(2) % velocity
clf;
semilogx(0, kugle.hastighed(2), 'or')
xlabel('t')
ylabel('v_y(t)')
hold on

% analytic

t = 0:dt:100000*dt;
v_ana =@(t) -sqrt( (2*kugle.masse.*g)./(rho.*kugle.A.*kugle.C_d) ) .* tanh(t.*sqrt((g.*rho.*kugle.C_d.*kugle.A)./(2.*kugle.masse)));
figure(2)
plot(t, v_ana(t),'-k')

% Numerical
for i=1:100000
    
    % Atmospheric drag
    dv_drag =  dt .* 1/kugle.masse * ( 0.5 * rho .* kugle.hastighed.^2 * kugle.C_d * kugle.A );
    
    % Gravitational attraction
    dv_gravity  = dt .* a;
    
    % Step updates
    kugle.hastighed = kugle.hastighed + dv_gravity + dv_drag; % Gravitional change in velocity    
    
    kugle.position = kugle.position + dt .* kugle.hastighed;
    
    if (mod(i,1000) == 0) || (mod(i,400) == 0 && (i < 10000)) % Precision plotting in the beginning, a bit less points afterwards
        figure(1)
        plot(i*dt, kugle.position(2),'o')
        figure(2)
        plot(i*dt, kugle.hastighed(2), 'or')
    end
end