clear;
clc
tic;
global r %距日距离
global len %整段绳长
global m %绳质量
global M1 %主星质量
global M2 %副星质量
global d %系绳直径
global E
global sgm %应力
global loss %阻尼损耗
global V %系绳电势
global n0 %电子密度
global v %太阳风速度
global mm %太阳风离子质量
global vac %真空介电常数
global qsw %库仑力强度
global rad %自旋频率
global M
global L0
global time %系统时间
global n %离散数目
global K %每段刚度
global C %每段阻尼
global Fsw %库仑力
global Fw %离心力
global N %拉力
global p
global y0
%%
%参数赋值
n=6;
r=1.49597870e11;
len=1e3;
m=0.01;
M1=120;
M2 =0.1;
d=5e-5;
E=7e10;
sgm=1e7;
loss=1e-4;
V =2e4; %电势大小
n0=7.3e6; %粒子密度
v =4e5; %粒子速度
mm=1.67e-27;
vac=8.854e-12;
rad=0.005;
time=0;
L0=[125,250,250,250,125];
M=[M1,m/4,m/4,m/4,m/4,M2];
K=zeros(1,n-1);
for i=1:n-1
K(i)=E*pi*d^2/(4*L0(i));
end
C=zeros(1,n-1);
for i=1:n-1
C(i)=(E*pi*d^2*m*loss/(4*len*L0(i)^2))^(1/2);
end
%K1=1.099557428756428;K2=0.549778714378214;K3=0.549778714378214;K4=0.549778714378214;K5=1.099557428756428;
%C1=2.965882571858067e-06;C2=1.482941285929034e-06;C3=1.482941285929034e-06;C4=1.482941285929034e-06;C5=2.965882571858067e-06;
qsw =0.18*V*sqrt(vac*n0*mm*v^2);
%%
%y0为初始位置和速度 x,y,z,x*,y*,z*.
y0=zeros(1,6*n);
y0(9)=125;
y0(15)=375;
y0(21)=625;
y0(27)=875;
y0(33)=1000;
%%
%库伦力
Fsw=zeros(3,n-2);
for j=1:n-2
Fsw(1,j)=-qsw*len/4;
end
%离心力
Fw=zeros(3,n);
Fw(3,1)=-(m/2+M2)*(rad^2)*len;
Fw(3,n)=M2*(rad^2)*len;
p=[0, 125, 375, 625, 875, 1000];
for h=2:n-1
Fw(3,h)=M(2)*(rad^2)*p(h);
end
% %绳拉力
% N=zeros(3,n-1);
% for i=1:n-1
%
% N(1,i) = (K(i)*(sqrt((y(i*n+1)-y((i-1)*n+1))^.2+(y(i*n+2)-y((i-1)*n+2))^.2+(y(i*n+3)-y((i-1)*n+3))^.2)-L0(i))+C(i)*(sqrt((y(i*n+4)-y((i-1)*n+4))^.2+(y(i*n+5)-y((i-1)*n+5))^.2+(y(i*n+6)-y((i-1)*n+6))^.2)))...
% *(y(i*n+1)-y((i-1)*n+1))/(sqrt((y(i*n+1)-y((i-1)*n+1))^.2+(y(i*n+2)-y((i-1)*n+2))^.2+(y(i*n+3)-y((i-1)*n+3))^.2));
% N(2,i) = (K(i)*(sqrt((y(i*n+1)-y((i-1)*n+1))^.2+(y(i*n+2)-y((i-1)*n+2))^.2+(y(i*n+3)-y((i-1)*n+3))^.2)-L0(i))+C(i)*(sqrt((y(i*n+4)-y((i-1)*n+4))^.2+(y(i*n+5)-y((i-1)*n+5))^.2+(y(i*n+6)-y((i-1)*n+6))^.2)))...
% *(y(i*n+2)-y((i-1)*n+2))/(sqrt((y(i*n+1)-y((i-1)*n+1))^.2+(y(i*n+2)-y((i-1)*n+2))^.2+(y(i*n+3)-y((i-1)*n+3))^.2));
% N(3,i) = (K(i)*(sqrt((y(i*n+1)-y((i-1)*n+1))^.2+(y(i*n+2)-y((i-1)*n+2))^.2+(y(i*n+3)-y((i-1)*n+3))^.2)-L0(i))+C(i)*(sqrt((y(i*n+4)-y((i-1)*n+4))^.2+(y(i*n+5)-y((i-1)*n+5))^.2+(y(i*n+6)-y((i-1)*n+6))^.2)))...
% *(y(i*n+3)-y((i-1)*n+3))/(sqrt((y(i*n+1)-y((i-1)*n+1))^.2+(y(i*n+2)-y((i-1)*n+2))^.2+(y(i*n+3)-y((i-1)*n+3))^.2));
% end
w=zeros(1,6*n);
for i=1:6*n
w(i)=1e-8;
end
options=odeset('RelTol',1e-4,'AbsTol',w);
[t,y]=ode45('one_km',0:1:3600,y0,options);
%%
% load matlab.mat
% save gg;
toc;
a1=[y(1,1) y(1,7) y(1,13) y(1,19) y(1,25) y(1,31)];
a2=[y(1,3) y(1,9) y(1,15) y(1,21) y(1,27) y(1,33)];
b1=[y(611,1) y(611,7) y(611,13) y(611,19) y(611,25) y(611,31)];
b2=[y(611,3) y(611,9) y(611,15) y(611,21) y(611,27) y(611,33)];
c1=[y(1250,1) y(1250,7) y(1250,13) y(1250,19) y(1250,25) y(1250,31)];
c2=[y(1250,3) y(1250,9) y(1250,15) y(1250,21) y(1250,27) y(1250,33)];
g1=[y(1900,1) y(1900,7) y(1900,13) y(1900,19) y(1900,25) y(1900,31)];
g2=[y(1900,3) y(1900,9) y(1900,15) y(1900,21) y(1900,27) y(1900,33)];
e1=[y(2480,1) y(2480,7) y(2480,13) y(2480,19) y(2480,25) y(2480,31)];
e2=[y(2480,3) y(2480,9) y(2480,15) y(2480,21) y(2480,27) y(2480,33)];
f1=[y(3120,1) y(3120,7) y(3120,13) y(3120,19) y(3120,25) y(3120,31)];
f2=[y(3120,3) y(3120,9) y(3120,15) y(3120,21) y(3120,27) y(3120,33)];
h1=[y(3600,1) y(3600,7) y(3600,13) y(3600,19) y(3600,25) y(3600,31)];
h2=[y(3600,3) y(3600,9) y(3600,15) y(3600,21) y(3600,27) y(3600,33)];
figure(1);
plot(t,y(:,1),'-b','LineWidth',2);
axis([0 3600 -30 0]);
grid on
title('近端x方向位移');xlabel('时间(s)');ylabel('R0x(m)');
set(gcf,'color',[1 1 1])
figure(2);
plot(t,y(:,3),'-b','LineWidth',2);
axis([0 3600 -0.005 0.02]);
grid on
title('近端z方向位移');xlabel('时间(s)');ylabel('R0z(m)');
set(gcf,'color',[1 1 1])
figure(3);
plot(t,y(:,31),'-b','LineWidth',2);
axis([0 3600 -200 0]);
grid on
title('远端x方向位移');xlabel('时间(s)');ylabel('R5x(m)');
set(gcf,'color',[1 1 1])
figure(4);
plot(t,y(:,33),'-b','LineWidth',2);
axis([0 3600 975 1005]);
grid on
title('远端z方向位移');xlabel('时间(s)');ylabel('R5z(m)');
set(gcf,'color',[1 1 1])
figure(5);
plot(a1,a2,'-k*',b1,b2,':k*',c1,c2,'-r*',g1,g2,':r*',e1,e2,'-b*',f1,f2,'--b*',h1,h2,'-g*','LineWidth',2,'markersize',4);
legend(['0s',sprintf('\n')],['611s',sprintf('\n')],['1250s',sprintf('\n')] ,['1900s',sprintf('\n')],['2480s',sprintf('\n')],['3120s',sprintf('\n')] ,['3600s',sprintf('\n')],'location','best');
axis([-200 0 0 1000]);
grid on
title('各点x-z平面振幅');xlabel('X0(m)');ylabel('Z0(m)');
set(gcf,'color',[1 1 1])
figure(6);
plot(t,radd(:,1),':b',t,radd(:,2),'--r',t,radd(:,3),'-k','LineWidth',2,'markersize',4);
s=legend(['0.005rad/s',sprintf('\n')],['0.006rad/s',sprintf('\n')],['0.007rad/s',sprintf('\n')],'location','best');
axis([0 3600 -30 0]);
grid on
title('不同自旋频率下的近端x方向位移');xlabel('时间(s)');ylabel('R0x(m)');
set(gcf,'color',[1 1 1])
figure(7);
plot(t,radd(:,4),':b',t,radd(:,5),'--r',t,radd(:,6),'-k','LineWidth',2,'markersize',4);
s=legend(['0.005rad/s',sprintf('\n')],['0.006rad/s',sprintf('\n')],['0.007rad/s',sprintf('\n')],'location','best');
axis([0 3600 -0.01 0.02]);
grid on
title('不同自旋频率下的近端z方向位移');xlabel('时间(s)');ylabel('R0z(m)');
set(gcf,'color',[1 1 1])
figure(8);
plot(t,radd(:,7),':b',t,radd(:,8),'--r',t,radd(:,9),'-k','LineWidth',2,'markersize',4);
s=legend(['0.005rad/s',sprintf('\n')],['0.006rad/s',sprintf('\n')],['0.007rad/s',sprintf('\n')],'location','best');
axis([0 3600 -250 0]);
grid on
title('不同自旋频率下的远端x方向位移');xlabel('时间(s)');ylabel('R5x(m)');
set(gcf,'color',[1 1 1])
figure(9);
plot(t,radd(:,10),':b',t,radd(:,11),'--r',t,radd(:,12),'-k','LineWidth',2,'markersize',4);
s=legend(['0.005rad/s',sprintf('\n')],['0.006rad/s',sprintf('\n')],['0.007rad/s',sprintf('\n')],'location','best');
axis([0 3600 975 1010]);
grid on
title('不同自旋频率下的远端z方向位移');xlabel('时间(s)');ylabel('R5z(m)');
set(gcf,'color',[1 1 1])
评论1