Code covered by the BSD License

# Plasticine simulation

### Maxim Vedenyov (view profile)

Plasticine circle by its own mass compressed on table. 2d

test3.m
tm=10;
dt=0.01;
ta=0:dt:tm;
Lta=length(ta);

d0=0.01; % distance of stability
ss=0.03; % square size
m=1e-3;
k=1; % spring coefficient
fcr=2*sqrt(k/m); % critical friction coeffcient
f=2*fcr; % friction coeficent

r=zeros(2,0);
g=[0; -10]; % gravitation

R=0.1;
R2=R^2;
ph=-0.1; % plate level
% make circle:
for yt=-R:d0:R
for xt=-R:d0:R
if (xt^2+yt^2<=R2)
r=[r [xt;yt]];
end
end
end
hr=plot(r(1,:),r(2,:),'k.');
hold on;
hpl=plot([-0.15 0.15],ph*[1 1],'r--');
axis equal;

v=zeros(size(r)); % velocity
N=size(r,2);
Na=1:N;

bi=false(1,N);
d=NaN(1,N);
Nx=NaN(N,N);
Ny=NaN(N,N);
for tc=2:Lta
t=ta(tc);
D=NaN(N,N); % distance matrix
% NaN when not need to calulate force
for n1=1:N-1
rt=r(:,n1);
ind2=n1+1:N;
bi=(abs(rt(1)-r(1,:))<=ss)&(abs(rt(2)-r(2,:))<=ss)&(Na>n1); % rest points in square
d(bi)=sqrt((rt(1)-r(1,bi)).^2+(rt(2)-r(2,bi)).^2);
D(n1,bi)=d(bi);
D(bi,n1)=D(n1,bi)';
d((d==0)&bi)=NaN; % points in the same place make no force between each other
Nx(n1,bi)=(r(1,bi)-rt(1))./d(bi);
Nx(bi,n1)=-Nx(n1,bi)';
Ny(n1,bi)=(r(2,bi)-rt(2))./d(bi);
Ny(bi,n1)=-Ny(n1,bi)';
end

%v(1,i1)=v(1,i1)+dt*sum((Nx(i1,:)*(D(i1,:)-d0))/m);
%Nx(i1,i2)*(v(1,i1)-v(1,i2))
%Ny(i1,i2)*(v(2,i1)-v(2,i2))
% pr=Nx(i1,i2)*(v(1,i1)-v(1,i2))+Ny(i1,i2)*(v(2,i1)-v(2,i2))
% Nx(i1,i2)*pr
% Ny(i1,i2)*pr
% v(1,i1)=v(1,i1)+(dt*f/m)*sum(Nx(i1,:).*(Nx(i1,:).*(v(1,i1)-v(1,:))+Ny(i1,i2)*(v(2,i1)-v(2,:))));

dvx=bsxfun(@minus,v(1,:),v(1,:)');

%A=(D-d0)/m; % acelerations from edges
%v=v+dt*;
end