Dynamics of Vibro-Impact system with impact condition

I am trying to solve the system of differential equations that describes the dynamics of Vibro-Impact system (with impact condition). Namely, dy/dt = x; dx/dt = (A cos(2*pi*F*t))/M; dz/dt = v; dv/dt = 0. Impact occurs at y-z>d+2*w-2*r and impact law, described by the velocities before v(+) and after v(-) an impact, is: v(+) = -m*v(-) +(m+1)*x.

Where y - absolute displacement of the outer structure; x - absolute velocity of the outer structure; z - absolute displacement of the ball; v - absolute velocity of the ball.

A, F, M, d, w, r, m are parameters.

Following Euler's scheme I tried to writ matlab cod, please help me frame the impact condition:

    %Code to solve Dynamics of Impact Pair
clear; clc
    %Parameters
A=0.01; %N
F=2; %Hz
M=124.5; %g
d=42;%mm
w=6;%mm
r=5;%mm
m=1;
    %Inputs
h=0.001;
t_final=55;
N=t_final/h;
    %Initial condition
t(1)=50;
y(1)=-A/(M*(2*pi*F)^2); %abs. displacement of the outer structure M
x(1)=0; %abs. velocity of the outer structure M
z(1)=0; %abs. displacement of the ball m 
v(1)=0; %abs. velocity of the ball m
    %update loop
for i=1:N
    t(i+1)=t(i)+h;
    y(i+1)=y(i)+h*(x(i));
    x(i+1)=x(i)+h*(A*cos(2*pi*F*t(i))/M);
    z(i+1)=z(i)+h*(v(i));
    v(i+1)=v(i);
    if abs(y(i+1)-z(i+1))>(d+2*w-2*r)/2, 
        v(i+1)=-m*v(i)+(m+1)*x(i+1);
        %switch k (i)= t(i+1)
            %k(i+1)=k(i)+h
    end
end
figure(1); clf(1)
plot(t,z)
xlabel ('Time (s)')
ylabel ('abs. displacement M')

if somebody knows how to use ODE23 with impact handled by the EVENT FUNCTION feature in this case the suggestions are also welcome. Thank you.

Answer 1

[Josh Meyer]

If you use ode45 then a chunk of your code becomes moot - you just need the initial conditions y0 and time span [t0 tfinal]. Your equations are encoded into a function that you pass to the solver. Note the use of anonymous functions here to pass the extra parameters to the ode and event functions.

%Inputs
tspan = [50 55];
y0 = [-A/(M*(2*pi*F)^2) 0 0 0];
opts = odeset('Events',@(t,y) impactEvents(t,y,d,w,r));
[t,y,te,ye,ie] = ode23(@(t,y) odefun(t,y,A,F,M,d,w,r,m),tspan,y0,opts);

function yxzv = odefun(t,y,A,F,M,d,w,r,m)
    dydt = y(2);
    dxdt = (A*cos(2*pi*F*t))/M;
    dzdt = y(4);
    dvdt = 0;
    yxzv = [dydt; dxdt; dzdt; dvdt];
end

The impact is controlled by the event function, which tells you when a certain expression is equal to zero. Due to the way your impact condition is expressed (one quantity larger than another), it needs to be rearranged to fit that format (one quantity equal to zero). Here is a quick first pass at what the event function might look like:

function [value,isterminal,direction] = impactEvents(t,y,d,w,r)
    value = y(1)-y(3)-d-2*w+2*r; %value that describes the event
    isterminal = 1; %should the integration stop when an event occurs?
    direction = 1; %From what direction will the value approach zero, from above or below?
end

One snag I noticed immediately is that the way the impact condition is expressed above (which might be wrong), no events ever occur. So you might want to double check the condition, equations, and/or direction.

>> any((y(:,1)-y(:,3))>(d+2*w-2*r))
ans =
  logical
   0

Once you get that working, you'll want to add some more code to restart the integration after an event occurs. The velocity after the impact will be the new initial condition for the velocity, and then the integration proceeds until it finishes or another event occurs. For an example of how to stop/restart based on events type edit ballode - this is an example of a bouncing ball.