how to simulate elastic collision

I'm trying to simulate elastic collisions for a school project and I cant figure out how to prevent over lap of the boundaries with each circle (boundary collisions). also, I cant figure out circle to circle collision. heres my collisions function sphere has the form [rad x_pos y_pos x_vel y_vel].

function [collision] = elasticCollision(ns,sphere,BC,density)

t = 0; dt =0.001; t_final = 100; u = zeros(ns,1);

while t <= t_final

collision_type = zeros(ns,1);

collision_flag = zeros(ns,1);

x_p1 = zeros(ns,1); y_p1 = zeros(ns,1);

t = t + dt;

for k = 1:1:ns

%creating moving circles

% Boundary test

x_p1(k) = sphere(k,4)*dt + sphere(k,2);

y_p1(k) = sphere(k,5)*dt + sphere(k,3);

if x_p1(k) - sphere(k,1) < BC(1)

collision_type(k) = 1; collision_flag(k) = 1;

elseif x_p1(k) + sphere(k,1) > BC(2)

collision_type(k) = 2; collision_flag(k) = 1;

elseif y_p1(k) - sphere(k,1) < BC(3)

collision_type(k) = 3; collision_flag(k) = 1;

elseif y_p1(k) + sphere(k,1) > BC(4)

collision_type(k) = 4; collision_flag(k) = 1;

else

collision_type(k) = 5; collision_flag(k)= 0;

end

switch(collision_type(k))

case(1)

u(k) = (sphere(k,1)+sphere(k,2))/(sphere(k,4));

case(2)

u(k) = (25 - sphere(k,2) - sphere(k,1))/(sphere(k,4));

case(3)

u(k) = (sphere(k,1)-sphere(k,3))/(sphere(k,5));

case(4)

u(k) = (25 - sphere(k,3) - sphere(k,1))/(sphere(k,5));

case (5)

u(k) = 1;

end

end

if max(collision_flag) ~= 0

[m , k] = min(u);

if min(u) == 1

m = dt;

end

for k = 1:1:ns

if collision_flag(k) == 1

x_p1(k) = sphere(k,4)*m + sphere(k,2);

y_p1(k) = sphere(k,5)*m + sphere(k,3);

sphere(k,2) = x_p1(k);

sphere(k,3) = y_p1(k);

end

end

for k = 1:1:ns

if collision_type(k) == 1 || collision_type(k) ==2

sphere(k,4) = -1 * sphere(k,4);

elseif collision_type(k) == 3 || collision_type(k) ==4

sphere(k,5) = -1 * sphere(k,5);

else

sphere(k,5) = sphere(k,5);

end

end

t = t - dt;

t = t + m;

end

pause(.1)

figure(1)

clf

hold on

for z = 1:1:ns

sphere(z,2) = sphere(z,4) + sphere(z,2);

sphere(z,3) = sphere(z,5) + sphere(z,3);

theta = 0: pi/360 : 2 * pi;

x_point = sphere(z,2) + sphere(z,1) * cos(theta);

y_point = sphere(z,3) + sphere(z,1) * sin(theta);

axis square

plot(x_point, y_point, 'r')

xlim([0 25])

ylim([0 25])

grid on

title(['collision Test t = ' num2str(t) '.']);

fill(x_point,y_point,'r')

end

hold off

end

% Object test

% solving for the fractional time step

% if k > 1

% object_collision = 0;

% for j = 1:1:(k - 1)

% distance = sqrt((x_2(j) - x_1(k))^2 + (y_2(j) - y_1(k))^2);

% x_2(j) = sphere(j,4)*dt + sphere(j,2);

% y_2(j) = sphere(j,5)*dt + sphere(j,3);

% if sphere(k,1) + sphere(j,1) == distance

% object_collision = 1;

% end

% if object_collision ~= 0

% a = sphere(k,2) - sphere(j,2); b = dt * (sphere(k,4) - sphere(j,4));

% c = sphere(k,3) - sphere(j,3); d = dt*(sphere(k,5) - sphere(j,5));

% e = (sphere(k,1) + sphere(j,1))^2;

% A = b^2 + d^2; B = 2*(a*b + c*d); C = a^2 + c^2 - e;

% u_1 = (-1*B + sqrt(B^2 - 4*A*C))/(2*A);

% u_2 = (-1*B - sqrt(B^2 - 4*A*C))/(2*A);

% if u_1 < 0

% u = u_2;

% elseif u_2 < 0

% u = u_1;

% elseif u_2 < u_1

% u = u_2;

% else

% u = u_1;

% end

% sphere(k,2) = sphere(k,4)*u*dt + sphere(k,2);

% sphere(k,3) = sphere(k,5)*u*dt + sphere(k,3);

% sphere(j,2) = sphere(j,4)*u*dt + sphere(j,2);

% sphere(j,3) = sphere(j,5)*u*dt + sphere(j,3);

% t_1 = atan(sphere(k,5)/sphere(k,4));

% t_2 = atan(sphere(j,5)/sphere(j,4));

% phi = atan((sphere(k,5)-sphere(j,5))/(sphere(k,4)-sphere(j,4)));

% v_1 = sqrt((sphere(k,4))^2 + (sphere(k,5))^2);

% v_2 = sqrt((sphere(j,4))^2 + (sphere(j,5))^2);

% m_1 = density*(4/3)*pi*sphere(k,1)^3;

% m_2 = density*(4/3)*pi*sphere(j,1)^3;

% sphere(k,4)=(v_1*cos(t_1-phi)*(m_1-m_2)+2*m_2*v_2*cos(t_2-phi)/(m_1+m_2))...

% *cos(phi)+v_1*sin(t_1-phi)*cos(phi+(pi/2));

% sphere(k,5)=(v_1*cos(t_1-phi)*(m_1-m_2)+2*m_2*v_2*cos(t_2-phi)/(m_1+m_2))...

% *sin(phi)+v_1*sin(t_1-phi)*sin(phi+(pi/2));

% sphere(j,4)=(v_2*cos(t_2-phi)*(m_1-m_2)+2*m_1*v_1*cos(t_1-phi)/(m_1+m_2))...

% *cos(phi)+v_2*sin(t_2-phi)*cos(phi+(pi/2));

% sphere(j,5)=(v_2*cos(t_2-phi)*(m_1-m_2)+2*m_1*v_1*cos(t_1-phi)/(m_1+m_2))...

% *sin(phi)+v_2*sin(t_2-phi)*sin(phi+(pi/2));

% end

% end

% end

end

Answer 1

[TED MOSBY]

Hi,

To prevent over lap of the boundaries with each circle and for circle to circle collision, you can make the following changes in your code:

Inside the main for k = … loop, right after t = t + dt, add:

% advance the centre once per step
sphere(k,2) = sphere(k,2) + sphere(k,4)*dt;   % x
sphere(k,3) = sphere(k,3) + sphere(k,5)*dt;   % y

Delete the temporary variables x_p1, y_p1, and any second position update later in the loop.

Add this block right after the code above:

% left / right
if sphere(k,2) - sphere(k,1) < BC(1)        % left wall
    sphere(k,2) = BC(1) + sphere(k,1);      % push back inside
    sphere(k,4) = -sphere(k,4);             % flip vx
elseif sphere(k,2) + sphere(k,1) > BC(2)    % right wall
    sphere(k,2) = BC(2) - sphere(k,1);
    sphere(k,4) = -sphere(k,4);
end
% bottom / top
if sphere(k,3) - sphere(k,1) < BC(3)        % bottom wall
    sphere(k,3) = BC(3) + sphere(k,1);
    sphere(k,5) = -sphere(k,5);             % flip vy
elseif sphere(k,3) + sphere(k,1) > BC(4)    % top wall
    sphere(k,3) = BC(4) - sphere(k,1);
    sphere(k,5) = -sphere(k,5);
end

This is done to re-position first so the centre is flush with the wall (no overlap) and to reverse only the velocity component normal to that wall. Together those two lines per wall enforce a perfectly elastic reflection and guarantee discs never leave the box.

In the for z = 1:ns plotting section remove:

sphere(z,2) = sphere(z,4) + sphere(z,2);
sphere(z,3) = sphere(z,5) + sphere(z,3);

Hope this helps!