Info
Diese Frage ist geschlossen. Öffnen Sie sie erneut, um sie zu bearbeiten oder zu beantworten.
Discrete element method array problem
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Antworten (0)
Diese Frage ist geschlossen.
Siehe auch
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!