Nested for loop not functioning correctly for 2D particle collisions
显示 更早的评论
Hello, I have very basic level knowledge of matlab and am trying to code 2D particle collisions with boundary walls and other particles. I am pretty sure that I have done the wall collisions correctly, because with just the wall collisions the graphs make sense. But when I add my particle collision code, the data gets really weird (ex: one particle has an x position of -0.0113 x 10^4 when it shouldn't be able to travel out of the boundary.)
I think it has something to do with the nested for loop, as when I adjust it, it starts giving me data that makes more sense, but I am unsure of what exactly is going wrong. If someone could take a look at this code and give suggestions as to what I could fix, that would be greatly appreciated!
(I apolgize if this is long for an ask, I just have no other way of getting feedback right now)
Initial Data (shouldn't have anything wrong here):
% number of particles
N = 3;
% set limits
ground = 0;
ceiling = 50;
left_wall = 0;
right_wall = 50;
collisions = 0;
c = zeros(1,1001);
% data for euler's method:
t = linspace(0,100,1001); % time from t=0 to t=100 with stepsize of 0.1s
% initial velocities for all particles
Vx = unifrnd(-10,10,N,1)
Vy = unifrnd(-10,10,N,1)
% initial positions for all particles
xi = unifrnd(5,35,N,1);
yi = unifrnd(5,35,N,1);
x_f = zeros(N,1001);
y_f = zeros(N,1001);
x_f(:,1) = xi;
y_f(:,1) = yi;
% c = (x1 - x2).^2 + (y1 - y2).^2; % condition for collision of p1 and p2
r1 = 8; % radius of p1
r2 = 8; % radius of p2
rad = (r1 + r2).^2; % radius of collision between any particle
Collisions:
figure
hold on;
for j = 1:numel(t)-1 % time
for i = 1:1:N % Nth particle
% wall collision for any particle
if x_f(i,j) <= left_wall || x_f(i,j) >= right_wall
Vx(i,1) = Vx(i,1) .* -1;
end
if y_f(i,j) >= ceiling || y_f(i,j) <= ground
Vy(i,1) = Vy(i,1).*-1;
end
% Particle-Particle Collision (PROBLEM AREA)
for K = 1:N
for P = K+1:N % I think this is what is going wrong!!!!!
c(j) = (x_f(K,j) - x_f(P,j)) .^2 + (y_f(K,j) - y_f(P,j)).^2; % this part also seems a bit incorrect
if c(j) <= rad
x1_coord = [x_f(K,j),y_f(K,j)];
x2_coord = [x_f(P,j), y_f(P,j)];
Vel_p1 = [Vx(K,1),Vy(K,1)];
Vel_p2 = [Vx(P,1), Vy(P,1)];
% changing velocity
Vel_p1 = Vel_p1 - ( ( (dot( (Vel_p1 - Vel_p2),(x1_coord - x2_coord) )) / (norm(x1_coord - x2_coord)) )* (x1_coord - x2_coord) );
Vel_p2 = Vel_p2 - ( ( (dot( (Vel_p2 - Vel_p2),(x2_coord - x1_coord) )) / (norm(x2_coord - x1_coord)) )* (x2_coord - x1_coord) );
% extracting Vx and Vy
Vx(K,1) = Vel_p1(1,1); Vy(K,1) = Vel_p1(1,2);
Vx(P,1) = Vel_p2(1,1); Vy(P,1) = Vel_p2(1,2);
collisions = collisions + 1;
end
end
end
x_f(i,j+1) = Vx(i,1).*0.1 + xi(i,1);
xi(i,1) = x_f(i,j+1);
y_f(i,j+1) = Vy(i,1).*0.1 + yi(i,1);
yi(i,1) = y_f(i,j+1);
end
end
What I'm using to plot:
figure(1)
for r = 1:1:N
plot(x_f(r,:),y_f(r,:),'r')
hold on
end
5 个评论
Torsten
2024-12-14
Maybe you could explain in your own words what you are trying to do in the code and how you try to achieve this. It's always difficult or - if the code is errornous - even impossible to reconstruct this from the code.
Anna
2024-12-14
What is your aim in each time step ? Determine x- and y-velocities for each particle that avoid a collision with any one of the other particles ?
What if a particle will collide with more than only one particle in the next time step ? Will the adjustments
% changing velocity
Vel_p1 = Vel_p1 - ( ( (dot( (Vel_p1 - Vel_p2),(x1_coord - x2_coord) )) / (norm(x1_coord - x2_coord)) )* (x1_coord - x2_coord) );
Vel_p2 = Vel_p2 - ( ( (dot( (Vel_p2 - Vel_p2),(x2_coord - x1_coord) )) / (norm(x2_coord - x1_coord)) )* (x2_coord - x1_coord) );
hinder collision with all particles satisfying c(j) <= rad ?
Anna
2024-12-14
Torsten
2024-12-15
At least the loop that detects wall collisions for the particles should be placed after the double loop that detects particle-particle collisions. Otherwise, necessary reversions in direction are simply overwritten in the particle-particle collision loop.
回答(0 个)
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
