Problem with solving discrete element method using leap frog method
I am new to writing matlab programming by using Discrete element method using leap frog algorithm. I got many error coming from my program. Can you all suggest me how to correct them with your all idea? Please let me hear your reply. n_part=4; kn=5; kt=2/7*kn; m=0.3; g=9.81; rad(1:n_part)=0.5; v_init=0.5; y=zeros(); % testing i_particle=1:n_part; % testing j_particle=n_part+1:2*n_part; position_x(1,1:n_part)=0.05; theta=2*pi*rand(size(position_x)); velocity_x=v_init*sin(theta); position_y(1,1:n_part)=0.05; velocity_y=v_init*cos(theta); timestep=100; dt=0.001; acceleration_x(:,1:n_part)=zeros(); acceleration_y(:,1:n_part)=zeros(); Fn=zeros(); Fn_i=zeros(); Fn_j=zeros(); v_half_x(1,1:n_part)=zeros(); v_half_y(1,1:n_part)=zeros(); for n=1:n_part for k=2:timestep % position_x v_half_x(k,n)=velocity_x(k-1,n)+0.5*dt*acceleration_x(k-1,n); position_x(k,n)=position_x(k-1,n)+v_half_x(k-1,n)*dt; % position_y v_half_y(k,n)=velocity_y(k-1,n)+0.5*dt*acceleration_y(k-1,n); position_y(k,n)=position_y(k-1,n)+v_half_y(k-1,n)*dt; for i=1:n_part for j=i+1:n_part if i>j % real position & distance lx=position_x(k-1,i)-position_x(k-1,j); ly=position_y(k-1,i)-position_y(k-1,j); root_xy=sqrt(ly^2+ly^2); % force calculation Fn=kn*root_xy^1.5; Fn_i=Fn_i+Fn; Fn_j=Fn_j+Fn; % acceleration term acceleration_x(k,:)=Fn_i./m; acceleration_y(k,:)=Fn_j./m; end end end velocity_y(k,n)=v_half_y(k-1,n)+0.5*dt*acceleration_y(k-1,n); velocity_x(k,n)=v_half_x(k-1,n)+0.5*dt*acceleration_x(k-1,n); end end
John Williams answered .
2025-11-20
The following gets the code working, but I've no idea if the results are meaningful!!
n_part=4;
kn=5;
kt=2/7*kn;
m=0.3;
g=9.81;
rad(1:n_part)=0.5;
v_init=0.5;
y=0;
% testing i_particle=1:n_part;
% testing j_particle=n_part+1:2*n_part;
position_x(1,1:n_part)=0.05;
theta=2*pi*rand(size(position_x));
velocity_x=v_init*sin(theta);
position_y(1,1:n_part)=0.05;
velocity_y=v_init*cos(theta);
timestep=100;
dt=0.001;
acceleration_x=zeros(timestep,n_part); %%%%%%%%%%%%%%
acceleration_y=zeros(timestep,n_part); %%%%%%%%%%%%%%
Fn=0;
Fn_i=0;
Fn_j=0;
v_half_x=zeros(timestep,n_part); %%%%%%%%%%%%%%
v_half_y=zeros(timestep,n_part); %%%%%%%%%%%%%%
for n=1:n_part
for k=2:timestep
% position_x
v_half_x(k,n)=velocity_x(k-1,n)+0.5*dt*acceleration_x(k-1,n);
position_x(k,n)=position_x(k-1,n)+v_half_x(k-1,n)*dt;
% position_y
v_half_y(k,n)=velocity_y(k-1,n)+0.5*dt*acceleration_y(k-1,n);
position_y(k,n)=position_y(k-1,n)+v_half_y(k-1,n)*dt;
for i=1:n_part
for j=i+1:n_part
%if i>j %%%%% i CANNOT be greater than j as you %%%%%%
%%%%% set j to be i+1 upwards! %%%%%%
lx=position_x(k-1,i)-position_x(k-1,j);
ly=position_y(k-1,i)-position_y(k-1,j);
root_xy=sqrt(ly^2+ly^2);
% force calculation
Fn=kn*root_xy^1.5;
Fn_i=Fn_i+Fn;
Fn_j=Fn_j+Fn;
% acceleration term
acceleration_x(k,:)=Fn_i./m;
acceleration_y(k,:)=Fn_j./m;
% end
end
end
velocity_y(k,n)=v_half_y(k-1,n)+0.5*dt*acceleration_y(k-1,n);
velocity_x(k,n)=v_half_x(k-1,n)+0.5*dt*acceleration_x(k-1,n);
end
end