Problem with solving discrete element method using leap frog method

15 ビュー (過去 30 日間)
Thin Rupar Win
Thin Rupar Win 2021 年 10 月 21 日
コメント済み: Thin Rupar Win 2021 年 10 月 21 日
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
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
  1 件のコメント
Thin Rupar Win
Thin Rupar Win 2021 年 10 月 21 日
Thank you very much for your idea. I acknowledge your help. Have a nice day.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeMatrix Computations についてさらに検索

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by