% DRAG EXAMPLE - Newton's 2nd Law
%ME 2040 - Matlab script to calculate the velocity of a falling
%sphere subject to aerodynamic forces as a function of time.
%The differential equation that is modeled here is:
% dv/dt=g-(c/m)v
% v is the velocity of the sphere
% t is time
% c drag coef.
% g gravitational accelleration cost.
% m mass of sphere
%In this example, the analytical solution solution, a numerical
%solution using Euler's method and the error between the two
%solutions are calculated and then plotted.
%
clc
clear
g=9.8 %m/s^2 gravitational constant
m=70 %kg mass of falling sphere 65
c=12.5 %kg/s drag coeff. of sphere
dt=2 %s time step
nt=24 % number of time steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% calculate analytical solution v1 %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxtime=(nt*dt) %duration of experiment (seconds)
t1=0:dt:maxtime; %generate time array
v_analytic =(g*m/c)*(1-exp(-t1*c/m)); %calculate velocity analytically
v_terminal(1:length(v_analytic))=m*g/c; %calculate terminal velocity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%generate time array for Euler method %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t2(1)=0; %initialize time
for i=2:nt+1
t2(i)=t2(i-1)+dt; %time array
end
% implement Euler's method for falling sphere
veuler(1)=0. %initialize velocity
for i=2:nt+1,
veuler(i)=veuler(i-1)+(g-veuler(i-1)*c/m)*(t2(i)-t2(i-1));
end
% calculate absolute error - difference between the analytical solution
% abd Euler Solution
for i=1:nt+1
error(i) = abs(v_analytic(i)-veuler(i));
end
% generate plot and labels
plot(t1,v_terminal,'-',t1,v_analytic,t2,veuler,'*',t2,error*10)
title('Sphere speed as a function of Time')
xlabel('Time (s)')
ylabel('V (m/s)')
key2 = ['Euler Method dt =', num2str(dt)]; %convert dt into a character string
%and concatonate with text to be
%used in legend
legend('Terminal Velocity','Analytic Solution',key2,'Absolute Error*10',0)