Discrete element method array problem

7 ビュー (過去 30 日間)
Thin Rupar Win
Thin Rupar Win 2021 年 10 月 20 日
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 件)

Community Treasure Hunt

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

Start Hunting!

Translated by