clear; clc; close all;
%% 参数设置
tMax = 300000; % 模拟时间步数
Ks=[5000 5000];
Kb=[500,500];
Nnodes=[30 40];
Nfila=length(Ks);
dt=0.01;
U=cell(1,Nfila); V=cell(1,Nfila);
Fx=cell(1,Nfila); Fy=cell(1,Nfila);
Fsx=cell(1,Nfila); Fsy=cell(1,Nfila);
Fbx=cell(1,Nfila); Fby=cell(1,Nfila);
DisL=cell(1,Nfila); DisR=cell(1,Nfila);
dis2=@(x1,y1,x2,y2) sqrt((x1-x2).^2 + (y1-y2).^2);
for i=1:Nfila
Rcx{i}=zeros(Nnodes(i),1)+(i-1.5)*40;
Rcy{i}=[0:Nnodes(i)-1]';
U{i}=Rcx{i}*0; V{i}=Rcy{i}*0;
Fx{i}=Rcx{i}*0; Fy{i}=Rcy{i}*0;
Fsx{i}=Rcx{i}*0; Fsy{i}=Rcy{i}*0;
Fbx{i}=Rcx{i}*0; Fby{i}=Rcy{i}*0;
DisL{i}=Rcx{i}*0; DisR{i}=Rcy{i}*0;
end
for j=1:tMax
%% 计算拉力,弯曲力
for i=1:Nfila
Fsx{i}=Rcx{i}*0; Fsy{i}=Rcy{i}*0;
Fbx{i}=Rcx{i}*0; Fby{i}=Rcy{i}*0;
Fx{i}=Rcx{i}*0; Fy{i}=Rcy{i}*0;
end
for i=1:Nfila
%拉力
for n=2:Nnodes(i)-1
disL{i}(n)=dis2(Rcx{i}(n-1),Rcy{i}(n-1),Rcx{i}(n),Rcy{i}(n));
disR{i}(n)=dis2(Rcx{i}(n+1),Rcy{i}(n+1),Rcx{i}(n),Rcy{i}(n));
Fsx{i}(n)= -Ks(i)*(disL{i}(n)-1)*(Rcx{i}(n)-Rcx{i}(n-1))/disL{i}(n)+...
-Ks(i)*(disR{i}(n)-1)*(Rcx{i}(n)-Rcx{i}(n+1))/disR{i}(n);
Fsy{i}(n)= -Ks(i)*(disL{i}(n)-1)*(Rcy{i}(n)-Rcy{i}(n-1))/disL{i}(n)+...
-Ks(i)*(disR{i}(n)-1)*(Rcy{i}(n)-Rcy{i}(n+1))/disR{i}(n);
end
n=Nnodes(i);
Fsx{i}(n)=Ks(i)*(disR{i}(n-1)-1)*(Rcx{i}(n-1)-Rcx{i}(n))/disR{i}(n-1);
Fsy{i}(n)=Ks(i)*(disR{i}(n-1)-1)*(Rcy{i}(n-1)-Rcy{i}(n))/disR{i}(n-1);
n=1;
Fsx{i}(n)=Ks(i)*(disL{i}(n+1)-1)*(Rcx{i}(n+1)-Rcx{i}(n))/disL{i}(n+1);
Fsy{i}(n)=Ks(i)*(disL{i}(n+1)-1)*(Rcy{i}(n+1)-Rcy{i}(n))/disL{i}(n+1);
%弯曲力
for n=1:Nnodes(i)
if n>=3 && n<=Nnodes(i)-2
Fbx{i}(n)=-Kb(i)*(Rcx{i}(n-2)-4*Rcx{i}(n-1)+6*Rcx{i}(n)-4*Rcx{i}(n+1)+Rcx{i}(n+2));
Fby{i}(n)=-Kb(i)*(Rcy{i}(n-2)-4*Rcy{i}(n-1)+6*Rcy{i}(n)-4*Rcy{i}(n+1)+Rcy{i}(n+2));
end
if n==1
Fbx{i}(n)=-Kb(i)*(Rcx{i}(n)+Rcx{i}(n+2)-2*Rcx{i}(n+1));
Fby{i}(n)=-Kb(i)*(Rcy{i}(n)+Rcy{i}(n+2)-2*Rcy{i}(n+1));
end
if n==Nnodes(i)
Fbx{i}(n)=-Kb(i)*(Rcx{i}(n)+Rcx{i}(n-2)-2*Rcx{i}(n-1));
Fby{i}(n)=-Kb(i)*(Rcy{i}(n)+Rcy{i}(n-2)-2*Rcy{i}(n-1));
end
if n==2
Fbx{i}(n)=-Kb(i)*(-2*Rcx{i}(n-1)+5*Rcx{i}(n)-4*Rcx{i}(n+1)+Rcx{i}(n+2));
Fby{i}(n)=-Kb(i)*(-2*Rcy{i}(n-1)+5*Rcy{i}(n)-4*Rcy{i}(n+1)+Rcy{i}(n+2));
end
if n==Nnodes(i)-1
Fbx{i}(n)=-Kb(i)*(-2*Rcx{i}(n+1)+5*Rcx{i}(n)-4*Rcx{i}(n-1)+Rcx{i}(n-2));
Fby{i}(n)=-Kb(i)*(-2*Rcy{i}(n+1)+5*Rcy{i}(n)-4*Rcy{i}(n-1)+Rcy{i}(n-2));
end
end
%% 位置更新
for n=3:Nnodes(i)
% Fx{i}(n)=Fsx{i}(n)+Fbx{i}(n)-0.01*U{i}(n)+0.01*sin(2*pi*j/13000);
Fx{i}(n)=Fsx{i}(n)+Fbx{i}(n)-0.02*U{i}(n);%+0.005*sin(2*pi*j/(30000));
Fy{i}(n)=Fsy{i}(n)+Fby{i}(n)-0.02*V{i}(n);
if n==Nnodes(i)
Fx{i}(Nnodes(i))=Fx{i}(Nnodes(i))+0.4*((sin(2*pi*j/10000)>0)-0.5);
end
U{i}(n)=U{i}(n)+Fx{i}(n)*dt;
V{i}(n)=V{i}(n)+Fy{i}(n)*dt;
Rcx{i}(n)=Rcx{i}(n)+U{i}(n)*dt;
Rcy{i}(n)=Rcy{i}(n)+V{i}(n)*dt;
end
end
%
%
if mod(j,200)==0
clf
clc
j
axis([-100 100,0 100])
set(gcf,'Color',[0.8 0.8 0.8])
hold on
for i=1:Nfila
plot(Rcx{i}(:),Rcy{i}(:),'r.')
% quiver(Rcx{i}(:),Rcy{i}(:),Fx{i}(:),Fy{i}(:))
end
hold off
drawnow
end
end