Info

Cette question est clôturée. Rouvrir pour modifier ou répondre.

Discrete element method array problem

1 vue (au cours des 30 derniers jours)
Thin Rupar Win
Thin Rupar Win le 20 Oct 2021
Clôturé : Sabin le 5 Jan 2023
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.

Réponses (0)

Cette question est clôturée.

Community Treasure Hunt

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

Start Hunting!

Translated by