Problem with solving discrete element method using leap frog method
显示 更早的评论
Dear Sir or Madam,
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
采纳的回答
Alan Stevens
2021-10-21
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Testing Frameworks 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
