Getting linear results from ODE45 instead of exponential

6 views (last 30 days)
I aam very new to Matlab and have this project to figure out the velocity vs time and displacement vs time of a sphere falling from rest. With my code I am able to get the displacement to work which is r(1) but my velocity vs time plot is linear when it should be exponential. I've tried retyping the whole equation for the code and inputting the squared term in all the ways I could figure out with no success.
function Projecto
clear all
tr=[0 15]; %seconds
initv=[0 0]; %starting point
[t,y]=ode45(@sphere, tr, initv);
plot(t,y(:,1))
ylabel('x (m)')
xlabel('time(s)')
figure
plot(t,y(:,2))
ylabel('velocity (m/s)')
xlabel('time(s)')
function r=sphere(t,y)
mass=167.6;
g=9.81;
rho=1.21; %kilogram per meter cubed
d=.01; %meters
r=zeros(2,1);
r(1)=y(2);
r(2)=(mass*g-rho*g*(pi/6)*d^3-0.5*rho*y(2)*abs(y(2))*(pi/4)*d^2*0.5)/(mass+(pi/12)*d^3*rho
  2 Comments
Ameer Hamza
Ameer Hamza on 30 Sep 2020
Why do you expect it to be exponential? For a sphere falling under constant gravity, the velocity is supposed to increase linearly.
Connor Jack
Connor Jack on 30 Sep 2020
I have to consider drag in the equation so I have a velocity squared term. In the equation r(2) the y(2) is my velocity variable, this is the equation that I started with:
m(dv/dt)=mg-rho*g*(pi/6 * d^3)-0.5rho*v^2*(pi/4 d^2)Cd-(pi/12)d^3*rho*(dv/dt)
rho*g*(pi/6 * d^3) - being the buoyant force
0.5rho*v^2*(pi/4 d^2)Cd - being the drag force
(pi/12)d^3*rho*(dv/dt) - being the force on accelerating body

Sign in to comment.

Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 30 Sep 2020
When you look at the result you have to start looking at how big the different components of the forces are along the trajectory. So write a function where you can calculate the different forces - that way you can judge how big the drag is relative to your force of gravity.
One more thing you can do is to extend the time of your fall - I increased it by a factor of 10 and then you start to see a reduction of the acceleration.
One further thing to try is to modify the mass and density (why are you using both?) to give you a lighter object that should have a lower terminal velocity, when I reduced that by a factor of 10 the terminal velocity can be clearly seen within 150 s.
HTH

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!