Discrete element method array problem

Question

Thin Rupar Win 2021-10-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1567613-discrete-element-method-array-problem

关闭： Sabin 2023-1-5

n_part=4;
kn=5;
kt=2/7*kn;
m=0.3;
g=9.81;
rad(1:n_part)=0.5;
y=zeros();
position_x=mode([1:n_part],5)+1;
velocity_x=rand(size(position_x))-0.5;
position_y=sort(position_x);
velocity_y=rand(size(position_y))-0.5;
w_x=rand(size(position_x));
w_y=rand(size(position_y));
y(1:6:6*n_part-5)=position_x;
y(2:6:6*n_part-4)=velocity_x;
y(3:6:6*n_part-3)=w_x;
y(4:6:6*n_part-2)=position_y;
y(5:6:6*n_part-1)=velocity_y;
y(6:6:6*n_part)=w_y;
y1=zeros();
y2=zeros();
timestep=100;
dt=0.001;
acceleration_x=zeros(1,n_part);
acceleration_y=zeros(1,n_part);
f=zeros();
v_half_x=rand(size(position_x));
v_half_y=rand(size(position_y));
% simulation is based on velocity verlet (or) leap frog method.
% inside is acceleratoin term only. not combine with other terms to find
% angular velocity.
for i=1:timestep
    % position x
    v_half_x(i+1)=velocity_x(i)+0.5*dt*acceleration_x(i);
    position_x(i+1)=position_x(i)+v_half_x(i)*dt;
    
    % position_y
    v_half_y(i+1)=velocity_y(i)+0.5*dt*acceleration_y(i);
    position_y(i+1)=position_y(i)+v_half_y(i)*dt;
    
    for i_part=1:n_part
        r_1=[y(6*i_part-5);y(6*i_part-2)];
        v_1=[y(6*i_part-4);y(6*i_part-1)];
        %v_half1=[v_half_x;v_half_y];
        w_1=[y(6*i_part-3);y(6*i_part)];
        rad1=rad(i_part);
        for j_part=i_part+1:n_part
            r_2=[y(6*j_part-5);y(6*j_part-2)];
            v_2=[y(6*j_part-4);y(6*j_part-1)];
            %v_half2=[v_half_x;v_half_y];
            w_2=[y(6*j_part-3);y(6*j_part)];   % later update with equation,
            rad2=rad(j_part);
            if(norm(r_1-r_2)< (rad(i_part)+rad(j_part)))
                % this is substituted equation for testing. the orginal
                % with only force term.
                unit_vector=(r_1-r_2)./norm(r_1-r_2);
                v_i_j=v_1-v_2+((rad1.*w_1)+(rad2.*w_2)).*unit_vector;
                v_i_j_t=-unit_vector.*(unit_vector.*v_i_j);
                t_i_j=v_i_j_t./norm(v_i_j_t);
                X_dot=(v_1-v_2)-((rad1.*w_1)+(rad2.*w_2)).*t_i_j;
                normal_vector=(v_1-v_2).*unit_vector;
                tangential_vector=((v_1-v_2).*t_i_j)-((rad1.*w_1)+(rad2.*w_2));
                F_n=kn.*normal_vector;                         % normal force
                F_t=kt.*tangential_vector;                     % Tangent force
                Contact_force=F_n+F_t;                         % contact force
                F_T=Contact_force+(m*g);
                acceleration_x(i+1)=F_T(1,:)./m;
                acceleration_y(i+1)=F_T(2,:)./m;
            end
        end
    end
    
    
    
    
    velocity_y(i+1)=v_half_y(i)+0.5*dt*acceleration_y(i+1);
    velocity_x(i+1)=v_half_x(i)+0.5*dt*acceleration_x(i+1);
    
    
    
end

Hi,

I am new to solving the problem with the matlab. As this is part of the assignment for the school,I would know how to get update to y() value with leap frog algorithm. Thank you very much for your help.