Info
此问题已关闭。 请重新打开它进行编辑或回答。
Discrete element method array problem
4 次查看(过去 30 天)
显示 更早的评论
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.
0 个评论
回答(0 个)
此问题已关闭。
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!