Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Particle Simulation in 3D, help please!

  1. Dec 3, 2009 #1
    I am trying to simulate a 3D particle simulation, where there are a 1000 particles of random x,y,z postn varying area between 2 n 4, and mass... There are different forces acting on each particles, Force due to heat, wind, and gravity, we calc acceleration from that, after which change in distance and hence new x,y,z positions are obtained. I have written most of it (at least I thnk I did),

    I need help updating these new x,y,z postion at time range (1:10), I am sort of stuck there.

    I also need to lock down this 3D space to a 100x100m box, it seems to go beserk,...

    and finally, these particles should get stuck to the top and bottom if they reach there... and if it goes beyond left or right, it should appear in the opposite side.

    Here's the code, let me knw what you guys think.

    NB : particle data file here : http://bit.ly/6qJaNb [Broken]

    Code (Text):



    %n = no.of particles  REPLACE 1000 with n

    data=load('particledata.dat');

    scatter3(data(:,1),data(:,2),data(:,3),data(:,5),'filled'), view(-60,60)

    m=data(:,4);

    g=9.81;

    x=data(:,1);
    y=data(:,2);
    z=data(:,3);

    F_w=zeros(3,1,1000);

    F_g=zeros(3,1,1000);

    F_h=zeros(3,1,1000);


    q=0.08;

    A=data(:,5);


               
               
    for t=2:1:10
        for w=1:1:1000;
         
                F_g_z(w)=-m(w)*g*z(w);   % force due to grav

                F_g(1,1,w)=0;
                F_g(2,1,w)=0;
                F_g(3,1,w)=F_g_z(w);

                rx=50-y(w);
                ry=x(w)-50;

                F_w_x(w)=1.6*(exp(-(sqrt(rx^2+ry^2)-25/20)^2))*rx;    % force wind
                F_w_y(w)=1.6*(exp(-(sqrt(rx^2+ry^2)-25/20)^2))*ry;
                F_w(1,1,w)=F_w_x(w);
                F_w(2,1,w)=F_w_y(w);
                F_w(3,1,w)=0;

                T(w)=1/16*(70*x(w)+(50-y(w))^2);

                F_h1(w)=q*A(w)*T(w)*z(w);             %force heat

                F_h(1,1,w)=0;
                F_h(2,1,w)=0;
                F_h(3,1,w)=F_h1(w);

                F=F_g+F_w+F_h;
           
                a=F/m(w);                     % acceleration
               
               
                u(1,1,w,1)=a(1,1,w);
                u(2,1,w,1)=a(2,1,w);
                u(3,1,w,1)=a(3,1,w);
               
                u(1,1,w,t)=u(1,1,w,t-1)+a(1,1,t)*t;    % to store u in x direction, of w-th particle with time t
                u(2,1,w,t)=u(2,1,w,t-1)+a(2,1,t)*t;
                u(3,1,w,t)=u(3,1,w,t-1)+a(3,1,t)*t;
               
    %             end
       
               
    %             fprintf('%u\n',u);

                s(1,1,w,t)=u(1,1,w,t)*t+(((1/2)*a(1,1,w))*t^2);   % s=ut+1/2at^2;
                s(2,1,w,t)=u(2,1,w,t)*t+(((1/2)*a(2,1,w))*t^2);
                s(3,1,w,t)=u(3,1,w,t)*t+(((1/2)*a(3,1,w))*t^2);
               
               
    %             ds=s_fin-s_ini;
                ds(1,1,w,t)=s(1,1,w,t)-s(1,1,w,t-1);
                ds(2,1,w,t)=s(2,1,w,t)-s(2,1,w,t-1);
                ds(3,1,w,t)=s(3,1,w,t)-s(3,1,w,t-1);

                x_change(w,t)=ds(1,1,w,t);
                y_change(w,t)=ds(2,1,w,t);
                z_change(w,t)=ds(3,1,w,t);
               
                x_final(w,1)=0;
                y_final(w,1)=0;
                z_final(w,1)=0;
               
                x_final(w,t)=x_final(w,t-1)+x_change(w,t);

                y_final(w,t)=y_final(w,t-1)+y_change(w,t);

                z_final(w,t)=z_final(w,t-1)+z_change(w,t);
               
                x(w)=x_final(w,t);
                y(w)=y_final(w,t);
                z(w)=z_final(w,t);
               
    %             fprintf('%z\n',z(w));
           
    % %            
    %             scatter3(x_final(w,t),y_final(w,t),z_final(w,t),'filled'), view(-60,60)
    %             pause(0.1);
               
               
               
               
    %             NEED TO UPDATE particledata / 3D to change accordingly
    %             HOW TO?
    %             ballp(i)=(x(i),y(i),z(i))
    %             dt = 0.5
    %             t=0.0
    %
    %
    %             while True:
    %
    %               t = t + dt
    %               ballpos(i) = ballpos(i) + (ballp(i)/m(i))*dt
    %
    %               if not (side > ball.x > -side):   (0,Y,Z) OR (100,Y,Z)
    %                 ball.p.x = -ball.p.x
    %               if not (side > ball.y > -side):   (X,0,Z) OR (X,100,Z)
    %                 ball.p.y = -ball.p.y
    %               if not (side > ball.z > -side):
    %                 ball.p.z = -ball.p.z   HERE BALL IS STUCK TO TOP N BOTTOM  for any (X,Y,0) or (X,Y,100)

     

        end
               
    end

    %
    % for w=1:1000
    %     for t=2:1:10
    %
    %     scatter3(x_final(w,t),y_final(w,t),z_final(w,t),'filled'), view(-60,60)
    %     pause(0.1);
    %     end
    % end
    %  
    %    

       
           



     
     
    Last edited by a moderator: May 4, 2017
  2. jcsd
  3. Dec 3, 2009 #2
    I dont understant what you means :( sry. When I glance at your code, I see that you can find x/y/z_final. You just need to draw a figure again to obtain new pos.

    do you mean you want to see fig in 2d? If that, view(2) can help you.

    It calls periodic condition. You can do this simply. Lets assume that z-pos ranges from 0 to 100 and there is another matrix to indicate whether particles have its z-pos out of this range (0 - no, 1 - yes). If (1), skip calculating this particle.

    For left-right side, you can put a condition after calculating and modifying data for this special move. It should be more complicated.

    p/s: I can only say this because I dont have time to read your code and completely understand it. Ill try it again when I am free :D
     
  4. Dec 3, 2009 #3


    Let me elaborate... I have x,y,z position of a 1000 particles. Now, over time these particles will be moving, due to the different forces applied on them, I want to simulate how these particles move around this 100x100x100m box. The t(1:10) was just time range going form 1 to 10, in steps of 1.

    the 100mx100m was a mistake, it is supposed to be 100x100x100m

    I have updated the code, to the one previously used, to make things a little clear.

    Thanx!
     
  5. Dec 4, 2009 #4
    It is a simple code that i tried to write. In this code, I used your coordinate data and created new position by using random integers. If z-pos is out of [0 100] range, this particle will not move next time.

    Code (Text):

    clc
    clear all

    %
    data = load('particledata.dat');
    x = data(:, 1);
    y = data(:, 2);
    z = data(:, 3);
    a_z = zeros(length(z), 1);

    scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

    for t = 1 : 100
        for j = 1 : length(x)
            if a_z(j) == 1
                continue
            end
           
            d = randi(10, 1, 3);
            s = randi(2, 1, 3);
           
            x(j) = x(j) + d(1) * (-1)^s(1);
            y(j) = y(j) + d(2) * (-1)^s(2);
            z(j) = z(j) + d(3) * (-1)^s(3);
           
            if x(j) <= 0 || x(j) >= 100
                C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
                x(j) = 100 * C;
            end
           
            if y(j) <= 0 || y(j) >= 100
                C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
                y(j) = 100 * C;
            end
           
            %if p is at bottom or top, it cant move anymore.
            if z(j) < 0 || z(j) > 100
                C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
                z(j) = 100 * C;
                a_z(j) = 1;
            end
        end
       
        scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
        pause(1)
    end
     
    Hope this helps.
     
  6. Dec 4, 2009 #5

    Thnx for the code, Its seems pretty straight forward with random x,y,z positions, its just when trying to implement forces and mapping them wrt to time, it seemed a bit too confusing... I will give it a ago, and come bak if I have any problems, Thanx a lot for ur help!
     
  7. Dec 4, 2009 #6

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    logix88,

    Your code has a *big* problem in the way you are propagating the particles from one time step to the next. If you want realistic physics, those particles need to have position and velocity. The acceleration changes the velocity, and the velocity in turn changes the position.

    Make sure you do things in that order. If you update the position first and velocity second you will have what is called the Euler integration technique. Simply switch the order (velocity, then position) and you will have what is called symplectic Euler integration technique. The difference between the two is immense. Basic Euler integration is truly abysmal. Symplectic Euler, while not that great, represents a vast improvement in realism.
     
  8. Dec 4, 2009 #7
    I thought thats what I did, first I calculate Forces due to heat, wind, gravity for each particle, and I get acceleration from those values. Initial vel is 0, I then use the acceleration to obtain u=u0+at and then plug that into s=ut+1/2*at^2. Then obtain the change in distance for each x,y,z between previous time step and current time step and finally update the particles position. Do you mean, that I should update final position before distance? What other way can I do that?


    Code (Text):

    clc
    clear all

    %
    data = load('particledata.dat');
    x = data(:, 1);
    y = data(:, 2);
    z = data(:, 3);
    m = data(:, 4);
    a_z = zeros(length(z), 1);
    g = 9.81;
    q = 0.08;

    A=data(:,5);

    F_w=zeros(3,1,1000);

    F_g=zeros(3,1,1000);

    F_h=zeros(3,1,1000);

    u=zeros(3,1,1000,100);

    scatter3(x, y, z, data(:,5),'filled'), view(-60,60)

    for t = 2 : 100
        for j = 1 : length(x)
    %         if a_z(j) == 1
    %             continue
    %         end
           
           
            F_g_z(j)=-m(j)*g*z(j);

            F_g(1,1,j)=0;
            F_g(2,1,j)=0;
            F_g(3,1,j)=F_g_z(j);

            rx=50-y(j);
            ry=x(j)-50;

            F_w_x(j)=1.6*(exp(-(((sqrt(rx^2+ry^2))-25)/20))^2)*rx;
            F_w_y(j)=1.6*(exp(-(((sqrt(rx^2+ry^2))-25)/20))^2)*ry;
            F_w(1,1,j)=F_w_x(j);
            F_w(2,1,j)=F_w_y(j);
            F_w(3,1,j)=0;

            T(j)=1/16*(70*x(j)+(50-y(j))^2);

            F_h1(j)=q*A(j)*T(j)*z(j);

            F_h(1,1,j)=0;
            F_h(2,1,j)=0;
            F_h(3,1,j)=F_h1(j);

            F=F_g+F_w+F_h;
        %    
        %     a=a_g+a_w+a_h;
            a=F/m(j);

            u(1,1,j,1)=a(1,1,j);
            u(2,1,j,1)=a(2,1,j);
            u(3,1,j,1)=a(3,1,j);

            u(1,1,j,t)=u(1,1,j,t-1)+a(1,1,j)*t;
            u(2,1,j,t)=u(2,1,j,t-1)+a(2,1,j)*t;
            u(3,1,j,t)=u(3,1,j,t-1)+a(3,1,j)*t;

    %             end


    %             fprintf('%u\n',u);

            s(1,1,j,t)=u(1,1,j,t)*t+(((1/2)*a(1,1,j))*t^2);
            s(2,1,j,t)=u(2,1,j,t)*t+(((1/2)*a(2,1,j))*t^2);
            s(3,1,j,t)=u(3,1,j,t)*t+(((1/2)*a(3,1,j))*t^2);


    %             ds=s_fin-s_ini;
            ds(1,1,j,t)=s(1,1,j,t)-s(1,1,j,t-1);
            ds(2,1,j,t)=s(2,1,j,t)-s(2,1,j,t-1);
            ds(3,1,j,t)=s(3,1,j,t)-s(3,1,j,t-1);

            x_change(j,t)=ds(1,1,j,t);
            y_change(j,t)=ds(2,1,j,t);
            z_change(j,t)=ds(3,1,j,t);

            x_final(j,1)=0;
            y_final(j,1)=0;
            z_final(j,1)=0;

            x_final(j,t)=x_final(j,t-1)+x_change(j,t);

            y_final(j,t)=y_final(j,t-1)+y_change(j,t);

            z_final(j,t)=z_final(j,t-1)+z_change(j,t);

            x(j)=x_final(j,t);
            y(j)=y_final(j,t);
            z(j)=z_final(j,t);
           
                   
            if x(j) <= 0 || x(j) >= 100
                C = (abs(x(j)) - x(j)) / (2 * abs(x(j)));
                x(j) = 100 * C;
            end
           
            if y(j) <= 0 || y(j) >= 100
                C = (abs(y(j)) - y(j)) / (2 * abs(y(j)));
                y(j) = 100 * C;
            end
           
            %if p is at bottom or top, it cant move anymore.
            if z(j) < 0 || z(j) > 100
                C = (abs(z(j)) - z(j)) / (2 * abs(z(j)));
                z(j) = 100 * C;
                a_z(j) = 1;
            end
        end
       
        scatter3(x, y, z, data(:,5),'filled'), view(-60,60)
        pause(0.1)
    end

     
     
  9. Dec 4, 2009 #8

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    This assumes constant acceleration. You don't have constant acceleration.

    What you should be doing is something like this:

    [tex]\aligned
    \vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\
    \vec x(t+\Delta t) &= \vec x(t) + \vec u(t+\Delta t) \Delta t
    \endaligned[/tex]

    You going from t=2 to 100 (and why 2?) in steps of 1. So in your case, [itex]\Delta t = 1[/itex].


    Do not do this:

    [tex]\aligned
    \vec u(t+\Delta t) &= \vec u(t) + \vec a(t)\Delta t \\
    \vec x(t+\Delta t) &= \vec x(t) + \vec u(t)\Delta t
    \endaligned[/tex]

    That second one is the basic Euler method. The first one is symplectic Euler.
     
  10. Dec 4, 2009 #9

    u(1,1,j,1)=a(1,1,j);
    u(2,1,j,1)=a(2,1,j);
    u(3,1,j,1)=a(3,1,j);

    u(1,1,j,t+1)=u(1,1,j,t)+a(1,1,j)*t;
    u(2,1,j,t+1)=u(2,1,j,t)+a(2,1,j)*t;
    u(3,1,j,t+1)=u(3,1,j,t)+a(3,1,j)*t;


    x_final(j,1)=x(j);
    y_final(j,1)=y(j);
    z_final(j,1)=z(j);

    x_final(j,t+1)=x_final(j,t)+u(1,1,j,t+1);
    y_final(j,t+1)=y_final(j,t)+u(2,1,j,t+1);
    z_final(j,t+1)=z_final(j,t)+u(3,1,j,t+1);

    x(j)=x_final(j,t+1);
    y(j)=y_final(j,t+1);
    z(j)=z_final(j,t+1);

    Is that what it should be like?
     
    Last edited: Dec 4, 2009
  11. Dec 4, 2009 #10

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    No. Why are you doing this?
    Think of it this way: Suppose you had time going from 102 to 200? None of your forces depend on time. When time is zero is just a label. Moving that label doesn't change the physics.
     
  12. Dec 4, 2009 #11
    Hey, I am sry I dont understand what you said, I am not a physicist, I knw I have some prb understanding the simple physics involved in this.

    TBH, I am more confused now than before...
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Particle Simulation in 3D, help please!
  1. 3d graph help (Replies: 3)

  2. 3D plot help (Replies: 1)

Loading...