Ball bounce - State equations

10 views (last 30 days)
Leon Ellis
Leon Ellis on 15 Apr 2022
Edited: Jon on 14 Jun 2023
Good day, I am currently a student assigned with a practical to simulate a ball bouncing over time. This is done by obtaining the spring constant (K) and damping constant (b) of the ball along with the excelleration of the ball from the state equations. State equations:
While in the air: ma=-mg; Thus a=g
While compressing: ma=-mg+K*y+b*Vy; Thus a=-g+(K/m)*y+(b/m)*Vy - where Vy is the velocity in the y direction.
While expanding: ma=-mg+K*y-b*Vy; Thus a=-g+(K/m)*y-(b/m)*Vy - Direction of (b/m)*Vy changes due to the change in direction of velocity.
I tried applying these state equations in code to demonstrate the bounce of a ball in the y direction vs time, but no luck. This is the current code I'm trying to get work:
M=0.057;
hold off;
g=9.81;
dt=0.01;
y2=zeros(1e5,1);
v=0;
flag=0;
b=4.1892;
K=1.5276e+04;
%New k value???????
K2=K;
y(1)=1;
hold on;
i=2;
%BounceData is the tracker data!
for i=2:length(BounceData)
if y(i-1)>0
v=v-g*dt;
y(i)=y(i-1)+v*dt;
flag=0;
plot(dt*(i-1),y(i-1),'r*');
hold on;
end
if (flag==0 && y(i-1)<=0)
v=v+(-g-(b/M)*v-(K2/M)*y(i-1))*dt;
y(i)=y(i-1)+v*dt;
flag=1;
plot(dt*(i-1),y(i-1),'r*');
hold on;
end
if (flag==1 && y(i-1)<=0.01)
v=v+(-g+(b/M)*v+(K2/M)*y(i-1))*dt;
y(i)=y(i-1)+v*dt;
flag=2;
plot(dt*(i-1),y(i-1),'r*');
hold on;
end
if (flag==2 && y(i-1)<=0.01)
v=v+(-g+(b/M)*v+(K2/M)*y(i-1))*dt;
y(i)=y(i-1)+v*dt;
flag=0;
plot(dt*(i-1),y(i-1),'r*');
hold on;
end
xlim([0,4]);
ylim([0,2]);
end
The b and K value have been determined practically. I am looking for a graph like such:
But am getting this:
Advice/help would be much appreciated. Thanks in advance!
  1 Comment
Shrey Tripathi
Shrey Tripathi on 14 Jun 2023
What does BounceData in your code depict? It hasn't been defined, so please provide the format of the data, or its data type.

Sign in to comment.

Answers (1)

Jon
Jon on 14 Jun 2023
Edited: Jon on 14 Jun 2023
First I'm assuming you were told to use simple Euler integration, e.g. y(n+1) = y(n) + dy/dt*tStep, otherwise you could use for example one of MATLAB's ode solver's for this.
Give that you are using the Euler integration, Yyou will need to use a smaller steps size, for example dt = 0.001
You shouldn't have to switch the sign of the damping term based on whether it is compressing or expanding, this will be taken into account by the sign of the velocity. The force generated by the damping will always oppose the current direction of travel.
You should just accumulate your t, and y data in numSteps by 1 vectors and then plot after you complete the simulation.
Note that when I tried making the above changes, I got a decaying bounce as expected, but it has essentialy come to rest after only 0.15 s, and the amplitude of the bounce was on 10-4. You seem to be expecting a longer decay time, and larger amplitude. Maybe check the units, etc on your masses, damping coefficients etc
  3 Comments
Jon
Jon on 14 Jun 2023
You should be able to further clean up the code, I provided to guide you. In particular you no longer need the if on the compressing or expanding, as the equations are the same.
Jon
Jon on 14 Jun 2023
Edited: Jon on 14 Jun 2023
Ooops, I just realized I didn't start with your initial condion of y = 1, here it is corrected. Please ignore my earlier comment about the amplitude and duration
M=0.057;
% % % hold off;
g=9.81;
dt=0.001;
tFinal = 4;
y2=zeros(1e5,1);
v=0;
flag=0;
b=4.1892;
K=1.5276e+04;
%New k value???????
K2=K;
% % % hold on;
% % % i=2;
%BounceData is the tracker data!
numSteps = round(tFinal/dt);
t = zeros(numSteps,1); % preallocate
y = zeros(numSteps,1); % preallocate
y(1)=1;
figure
for i=2:numSteps
t(i) = t(i-1) + dt; % increment time vector
if y(i-1)>0
v=v-g*dt;
y(i)=y(i-1)+v*dt;
flag=0;
% % % plot(dt*(i-1),y(i-1),'r*');
% % % hold on;
end
if (flag==0 && y(i-1)<=0)
v=v+(-g-(b/M)*v-(K2/M)*y(i-1))*dt;
y(i)=y(i-1)+v*dt;
flag=1;
% % % plot(dt*(i-1),y(i-1),'r*');
% % % hold on;
end
if (flag==1 && y(i-1)<=0.01)
% % % v=v+(-g+(b/M)*v+(K2/M)*y(i-1))*dt;
v=v+(-g-(b/M)*v-(K2/M)*y(i-1))*dt;
y(i)=y(i-1)+v*dt;
flag=2;
% % % plot(dt*(i-1),y(i-1),'r*');
% % % hold on;
end
if (flag==2 && y(i-1)<=0.01)
v=v+(-g+(b/M)*v+(K2/M)*y(i-1))*dt;
y(i)=y(i-1)+v*dt;
flag=0;
% % plot(dt*(i-1),y(i-1),'r*');
% % hold on;
end
% % % xlim([0,4]);
% % % ylim([0,2]);
end
plot(t,y)

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!