% 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
%Two equivalent methods for handling arrays are shown.
g=9.81 %m/s^2 gravitational constant
m=68.1 %kg mass of falling sphere
c=12.5 %kg/s drag coeff. of sphere
dt=0.1 %s time step
nt=60 % number of time steps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% method without for loop %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
maxtime=(nt*dt) %duration of experiment (seconds)
t1=0:dt:maxtime; %generate time array - t1
v1=(g*m/c)*(1-exp(-t1*c/m)); %calculate velocity - v1
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% method using for loop where array %
% elements are explicitly defined %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t2(1)=0; %initialize time
for i=2:nt+1
t2(i)=t2(i-1)+dt; %time array
end
for i=1:nt+1,
v(i)=(g*m/c)*(1-exp(-t2(i)*c/m));
end
% generate plot and label
plot(t2,v,t1,v1,'*')
title('Sphere speed as a function of Time')
xlabel('Time (s)')
ylabel('V (m/s)')
% write data to an ASCII file
fid=fopen('sphere_vel.dat','w');
fprintf(fid,'%7.3f \t %7.3f \n',t1,v1); %c-style formatting
fclose(fid);