clc
clear
clc
clear
%% //--------------------------------Parameters-----------------------------------//
%radius
r = 5;
%mass
m = 1;
%length of boundary
l = 100;
%volume fraction
V = 0.1;
%number of particles
n = floor((V * l^3)/(r^3 * pi * (4/3)));
%initial speed
k = 0;
velocity = rand(n,3)*k;
%initial acceleration
acc = zeros(n,3);
%stiffness coeficient
k1 = 8;
%attenuation coefficient of speed
k2 = 0.05;
%step size
step = 0.1;
%threshold value of speed
th_v = 0.1;
%energy
Ecc = [];
Tcc = [];
T = 0;
%gif
flag = 0;
%initial positon of particles
pos(1,1) = (l - 2*r)*rand + r;
pos(1,2) = (l - 2*r)*rand + r;
pos(1,3) = (l - 2*r)*rand + r;
for i = 2 : n
pos(i,1) = (l - 2*r)*rand + r;
pos(i,2) = (l - 2*r)*rand + r;
pos(i,3) = (l - 2*r)*rand + r;
end
%% //--------------------------------Simulation process-----------------------------------//
tic
drawcount=0;
while 1
judge = 0;
for i = 1 : n %acceleration
acc(i,1) = 0;
acc(i,2) = 0;
acc(i,3) = 0;
end
%intereaction of particles
for i = 1 : n - 1
for j = i + 1 : n
if (pos(i,1) - pos(j,1))^2 + (pos(i,2) - pos(j,2))^2 + (pos(i,3) - pos(j,3))^2 < 4*r^2
judge = 1;
distance = sqrt((pos(i,1) - pos(j,1))^2 + (pos(i,2) - pos(j,2))^2 + (pos(i,3) - pos(j,3))^2);%distance between the centers of disk i and disk j
overlap = 2 * r - distance;%overlap
f = overlap * k1; %normal force
acc(i,1) = acc(i,1) + ((pos(i,1) - pos(j,1))/distance) * f/m;
acc(i,2) = acc(i,2) + ((pos(i,2) - pos(j,2))/distance) * f/m;
acc(i,3) = acc(i,3) + ((pos(i,3) - pos(j,3))/distance) * f/m;
acc(j,1) = acc(j,1) + ((pos(j,1) - pos(i,1))/distance) * f/m;
acc(j,2) = acc(j,2) + ((pos(j,2) - pos(i,2))/distance) * f/m;
acc(j,3) = acc(j,3) + ((pos(j,3) - pos(i,3))/distance) * f/m;
end
end
end
for i = 1 : n %refresh velocity
velocity(i,1) = (velocity(i,1) + acc(i,1) * step) * (1 - k2);
velocity(i,2) = (velocity(i,2) + acc(i,2) * step) * (1 - k2);
velocity(i,3) = (velocity(i,3) + acc(i,3) * step) * (1 - k2);
end
for i = 1 : n %refresh position
pos(i,1) = pos(i,1) + velocity(i,1) * step;
pos(i,2) = pos(i,2) + velocity(i,2) * step;
pos(i,3) = pos(i,3) + velocity(i,3) * step;
end
% flag_v = 0;
for i = 1 : n
%backward wall
if pos(i,1) < r
velocity(i,1) = velocity(i,1) + ((abs(pos(i,1) - r)) * k1) * step;
end
%forward wall
if pos(i,1) > l - r
velocity(i,1) = velocity(i,1) - ((abs(pos(i,1) - (l - r))) * k1) * step;
end
%left wall
if pos(i,2) < r
velocity(i,2) = velocity(i,2) + ((abs(pos(i,2) - r)) * k1) * step;
end
%right wall
if pos(i,2) > l - r
velocity(i,2) = velocity(i,2) - ((abs(pos(i,2) - (l - r))) * k1) * step;
end
%under wall
if pos(i,3) < r
velocity(i,3) = velocity(i,3) + ((abs(pos(i,3) - r)) * k1) * step;
end
%upper wall
if pos(i,3) > l - r
velocity(i,3) = velocity(i,3) - ((abs(pos(i,3) - (l - r))) * k1) * step;
end
% if velocity(i,1)^2 + velocity(i,2)^2 > th_v^2
% flag_v = 1;
% end
end
% figure(1)
% draw the container
clf
t = [0 0 0;
l 0 0;
l l 0;
0 l 0;
0 0 l;
l 0 l;
l l l;
0 l l;];
t1 = [1 2 3 4 1 5 6 2 6 7 3 4 8 5 8 7];
i1 = 1;
x = [];
y = [];
z = [];
for i=t1
x(i1) = t(i,1);
y(i1) = t(i,2);
z(i1) = t(i,3);
i1 = i1+1;
end
if drawcount==10
drawcount=0;
end
if drawcount==0;
plot3(x,y,z,'r');
hold on
axis equal
% draw the particles
for i = 1 : n
h = drawsphere(r , pos(i,1) , pos(i,2) , pos(i,3));
h.EdgeColor = [i/n,1-i/n , 0];
end
pause(0.01)
end
drawcount=drawcount+1;
%exit the loop
% if flag_v == 0
% break;
% end
if judge == 0
break
end
%Kinetic energy of the system
% Ecc = [Ecc;
% sum(1/2.*(velocity(:,1).^2+velocity(:,2).^2))];
% T = T+step;
% Tcc = [Tcc;
% T];
% % figure(2)
% clf
% plot(Tcc,Ecc); drawnow
%
% %gif
% G = getframe(gcf);
% I = frame2im(G);
% [I,map] = rgb2ind(I,256);
% if flag == 0
% imwrite(I,map,'test.gif','gif', 'Loopcount',inf,'DelayTime',0.1);%鐢讳竴甯?
% flag = 1;
% else
% imwrite(I,map,'test.gif','gif', 'WriteMode','append','DelayTime',0.1);%鐢讳竴甯?
% end
end
plot3(x,y,z,'r');
hold on
axis equal
% draw the particles
for i = 1 : n
h = drawsphere(r , pos(i,1) , pos(i,2) , pos(i,3));
h.EdgeColor = [i/n,1-i/n , 0];
end
pause(0.01)
toc
disp('缁撴潫');
