Ball bounce - State equations

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!

[Shrey Tripathi]

What does BounceData in your code depict? It hasn't been defined, so please provide the format of the data, or its data type.

Answer 1

[Jon]

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

[Jon]

Here is my modified version of your code. I commented out the portions of your original code that were not needed.

M=0.057;
hold off;
g=9.81;
dt=0.001;
tFinal = 0.4;
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!
numSteps = round(tFinal/dt);
t = zeros(numSteps,1); % preallocate
y = zeros(numSteps,1); % preallocate
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)

[Jon]

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]

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)