Need to stop particle interaction during simulation. Any ideas and hints would be useful

Hi I am trying to sumulate a process where the particle are not supposed to interact with each other. As you can see in the image, the particles overlap each other on the surface. Is there any way I could stop the particles from interacting with each other , ie, stop them from overlapping and form a continuous rows of circles, where the next level of circles would sit on top of the previous layer and not on the surface? Any help advice and suggestions would be useful. if you need more information, let me know.

[Walter Roberson]

How are you representing the particles? Is your simulation in MATLAB itself or in Simulink or SimState or ?

[Hari]

Yes its done completely in Matlab itself.

[Hari]

Is there any way that i can include the entire thing into a matrix lattice of nXn, all initialized to 0 and then when a ball sits on that position, that particular grid is converted to 1 and then include a code where the status of the grid is checked before the next ball sits on it?

[Marc]

There is many ways to approach this. I recommend looking into discrete element modeling. eDEM is a commercially available solution that is very $$$$$.

Can you show the code that leads to this plot?

Otherwise it is very hard to make suggestions.

[Hari]

function y = biliard2_1(x)

% User inputs

n = input('enter the number of particles ') ;                    % no of particles in the system
r=input('enter the radius of the particles  ')           ;             % radius of each particle
N=input('enter the number of iterations desired ');          % entering the number of iterations desired
A=input('enter the y-coordinate of the wave ');

% defining the square wave

t = 0:0.05:100 ;                                   
sq_wave=0.5*A;
sq_wave =sq_wave+ 0.5*A*square(2*pi*0.01*t);       % 'A' is the amplitude of the square wave (0 to A) equivalent to the surface of the silicon, % '50' is the frequency of the square waves equivalent to the spacing between the pores
w=(sq_wave)';                                   % converting the column matrix of the wave to a row vector
time=(t)';                                          % converting the column matrix of the time to a row vector,which is basically the x-axis
wave=[time,w];                                  % combining the x-y coordinates to a single vector which has the [x,y] coordinates of the wave.

dt = 0.5     ;             % speed of the particles
ma=1;
mb=1;
x0=[0 0 0 0];

% Calculations used in the code

l=   100;                 

p = (0.5*rand(n,2)+0.5)*l;           % rand(m,n) gives a matrix of a mXn size,
% multiplied with the mean free path length,
%determines the initial position of the
% particles in the 1st column being the
% x-axis and 2nd column being the y-axis
v = rand(n,2)  ;          % determines the speed of the individual particles, multiply with 'x' to get the actual velocity
set_flag=zeros(n,1);

for tt=1:N              % total number of iterations for which the simulation has to run for each run there is a new rand
    % number generated for the position of the
    % particles, change the x-axis from 0-50 to
    % 0-100

      plot(p(:,1),p(:,2),'r o',t,sq_wave,'b-','markersize',30);
      grid on;
      axis([0  l  0  l]) ;          % forming the plot with each particle being in the position according to 'p'
      ;
      drawnow;

      for i = 1 : n               % iterating over all the 'n' particles

          for j = i + 1 : n           % iterating over all the 'i+1' particle which is interacting with the other

              d = sqrt((p(i, 1) - p(j, 1)) ^ 2 + (p(i, 2) - p(j, 2)) ^ 2);  % calculating the diameter of 2 particles and then calculating the distance of separation between them

              if d < 2 * r                                                  % if the distance between the 2 particles is less than twice the radius, ie, they are colliding
                  va = sqrt(v(i,1)^2 + v(i,2)^2)  ;                        % calculating the x-velocity of each individual particle
                  vb = sqrt(v(j,1)^2 + v(j,2)^2) ;                          % calculating the y-velocity of each individual interacting particle
                  gama = atan((p(i,2)-p(j,2))/(p(i,1)-p(j,1))) ;            % calculating the inverse tangent in radians [pi/2,-pi/2] for each element. calculate the angle between the 2 interacting particles
                  alpha = atan2(v(i,1),v(i,2)) - gama    ;                  % P = atan2(Y,X) returns an array P of same size as X,Y containing the element-by-element, four-quadrant inverse tangent of Y and X
                  beta  = atan2(v(j,1),v(j,2)) - gama   ;                   % calculating the angle between the y-component of the velocity and
                  [x,fval] = fsolve(@gh,x0)   ;                             % Solves the function given in 'gh' starting from x0 and returns the value of the function at the solution x.

                  v(i,1)=x(1)*cos(x(3)+gama);                                 % uses the 'x' value obtained after solving the equation in the below function
                  v(i,2)=x(1)*sin(x(3)+gama);
                  v(j,1)=x(2)*cos(x(4)+gama);
                  v(j,2)=x(2)*sin(x(4)+gama);

                  p(i,1)=p(i,1)+v(i,1)*dt;                                            % x-position of 'i' particle = initial x-position of 'i' particle + distance covered in time 'dt' at speed v(i,1)
                  p(i,2)=p(i,2)+v(i,2)*dt ;                                           % y-position of 'i' particle = initial y-position of 'i' particle + distance covered in time 'dt' at speed v(i,2)

                  p(j,1)=p(j,1)+v(j,1)*dt ;                                           % x-position of 'j' particle ( other interacting particle) = initial x-position of 'j' particle + distance covered in time 'dt' at speed v(j,1)
                  p(j,2)=p(j,2)+v(j,2)*dt ;                                           % y-position of 'j' particle ( other interacting particle) = initial y-position of 'j' particle + distance covered in time 'dt' at speed v(j,1)
              end

          end

          if p(i, 1) < 0         % defines the dimesions within which the particles are moving randomly, beyond which the direction changes, then the direction of the particle is reversed to the opposite x-direction
              v(i, 1) = -v(i, 1);
          end
          if p(i, 1) > l
              v(i, 1) = -v(i, 1);
          end
          if p(i, 2) < 0          % if the y-coordinate of the each 'i' particle is (), then the direction of the particle is reversed to the opposite y-direction
              v(i, 2) = -v(i, 2);
          end
          if p(i, 2) > l
              v(i, 2) = -v(i, 2);
          end

          %%main code where the particles interact with the surface and this is the section where everyhting should be there

          if p(i,1)<55                        % checking the x-position of the particles and updating the position flag each time it crosses this mark

              set_flag(i)=true;

          else
              set_flag(i)=false;
          end

          if (p(i,2)<(A+5)&& set_flag(i)==true)           % the region on top of  the wave
              rand_nbr=rand(1);
              sticking_coeff=0.6;                         % possibility of adding an equation to determine the sticking coefficient, rs= 0.25*sticking_coeff*mean molecular velocity*conc of intermediates
              if rand_nbr>sticking_coeff                  % if the random number generated is greater than the sticking coeffiecient, then the particle will rebound into the gas phase
                  v(i,2)=-v(i,2);

              else
                  v(i,2)=0;                               % particle sticks to the surface and undergoes surface chemistry
                  v(i,1)=0;
                  rectangle('Position',[p(i,1),20,10,10],'facecolor','r');

              end
          elseif (p(i,2)<(A+5)&& set_flag(i)==false)          % the region of the trench, boucning off the bottom and then bouncing off
              v(i,2)=v(i,2);
              v(i,1)=v(i,1);
              %                         rand_nbr1=rand(1);
              %                         sticking_coeff1=0.4;                             % sticking coefficient at the bottom of the trench
              %                         if rand_nbr1>sticking_coeff1
              %                             v(i,2)=-v(i,2);
              %                         end

          end

      end

      p=p+v*dt  ;                                                         % final position = initial position + velocity*time  required for the movement of the particles

end

      function y = gh(x)

          e=1;
          y = [ma*va*cos(alpha) + mb*vb*cos(beta) - ma*x(1)*cos(x(3)) - mb*x(2)*cos(x(4));
              va*sin(alpha) - x(1)*sin(x(3));
              vb*sin(beta)  - x(2)*sin(x(4));
              x(1)*cos(x(3)) - x(2)*cos(x(4)) - e*(va*cos(alpha) - vb*cos(beta))];    % matrix of y of size (4X1) which returns the value of x in the above function
          % x0=[0 0 0 0] which then takes the value of [x] given by [y(1) y(2) y(3) y(4)]

      end
  end

Answer 1

[Hari]

Yes, its completely done in Matlab itself.