8．附录
1、第一问和第二问程序 160911
%系泊系统的第一问和第二问 %搜索方法 clear clc
close all
g=9.8;%重力加速度 v=36;%风速
rou=1025;%海水密度
ma=1000;da=2;ha=2;%浮标 lb=1;mb=10;db=0.05;%钢管 mc=100;lc=1;dc=0.3;%钢桶 ball=1200;%重物球
ld=0.105;md=7;Ld=22.05;nd=Ld/ld;md=md*ld*g;%链环 me=600;%锚
Ffb=rou*g*pi*db^2/4*lb;%钢管的浮力 mb=mb*g;%钢管的重力 t=1;
for ball=120000:2:3000
for h1=0.759%0.9:0.01:1.1%0.3:0.05:1%0.31:0.01:1.5%0.68:0.0005:0.72 for h2=0.781%h1:0.01:1.5%0.75:0.0005:0.79 L=h1; R=0;
%------浮标-------------------------
S1=da*(ha-(h1+h2)/2)*da/sqrt((h2-h1)^2+da^2); S2=pi*da^2/8*(h2-h1)/da; S=S1+S2;%风载荷
Fw=0.625*S*v^2;%风力
V=pi*da^2/8*(h1+h2);%吃水体积 Ffa=rou*g*V;%浮标的浮力 Ta1=Ffa-ma*g;%竖直方向 Ta2=Fw;%水平方向
Ta=sqrt(Ta1^2+Ta2^2);%浮标的总拉力
thetaa=atand(Ta1/Ta2); %浮标拉力的角度 alphaa=atand((h2-h1)/da);%浮标的倾斜角 %------浮标-------------------------
29
%------钢管------------------------- for k=1:4 Tb2(k)=Ta2;%水平方向 if k==1
Tb1(k)=Ta1+Ffb-mb;%竖直方向 else
Tb1(k)=Tb1(k-1)+Ffb-mb;%竖直方向 end
thetab(k)=atand(Tb1(k)/Tb2(k));%拉力角度 Tb(k)=sqrt(Tb1(k)^2+Tb2(k)^2);%拉力大小
if k==1%钢管的倾斜角
alphab(k)=atand((Ta*sind(thetaa)+Tb(k)*sind(thetab(k)))... /(Tb(k)*cosd(thetab(k))+Ta*cosd(thetaa))); else
alphab(k)=atand((Tb(k-1)*sind(thetab(k-1))+Tb(k)*sind(thetab(k)))... /(Tb(k)*cosd(thetab(k))+Tb(k-1)*cosd(thetab(k-1)))); end
if k==1
xb(k)=ld*cosd(alphaa);%钢管k的x长度 else
xb(k)=ld*cosd(alphab(k-1));%钢管k的x长度 end
R=R+xb(k);%目前游动总长度
L=L+sind(alphab(k))*lb; k=k+1; end
%------钢管-------------------------
%------钢桶------------------------- Ffc=rou*g*pi*dc^2/4*lc;%钢桶的浮力 Tc2=Ta2;%水平方向
if (Tb1(4)+Ffc-(mc+ball)*g)<0 break;
end
Tc1=Tb1(4)+Ffc-(mc+ball)*g;%竖直方向 Tc=sqrt(Tc1^2+Tc2^2);%钢桶的总拉力
thetac=atand(Tc1/Tc2); %钢桶拉力的角度
alphac=atand((Tc*sind(thetac)+Tb(4)*sind(thetab(4))+ball*g)... /(Tb(4)*cosd(thetab(4))+Tc*cosd(thetac)));%钢桶的倾斜角
30
% if alphac<85 % continue; % end
R=R+lc*cosd(alphac);%目前游动总长度 L=L+lc*sind(alphac);
%------钢桶-------------------------
%------链环------------------------- % for k=1:nd % if k==1
% Td(k)=Tc-md*sind(thetac);%拉力
% thetad(k)=(-md*cosd(thetac))/Tc+thetac;%拉力角度 % else
% Td(k)=Td(k-1)-md*sind(thetad(k-1));%拉力
% thetad(k)=-md*cosd(thetad(k-1))/Td(k-1)+thetad(k-1);%拉力角度
% end
% options = odeset('RelTol',1e-4,'AbsTol',[1e-8 1e-8]); % if k==1
% [A,B] =ode45(@rigid,[0,0.105],[Tc thetac/180*pi],options); % else
% [A,B] = ode45(@rigid,[nd*(k-1) nd*k],[Td(k-1),thetad(k-1)/180*pi],options); % end
% Td(k)=B(end,1);
% thetad(k)=B(end,2)*180/pi;
for k=1:nd
Td2(k)=Tc2;%水平方向拉力 if k==1
Td1(k)=Tc1-md;%竖直方向拉力 else
Td1(k)=Td1(k-1)-md;%拉力 end
thetad(k)=atand(Td1(k)/Td2(k));%拉力角度 Td(k)=sqrt(Td1(k)^2+Td2(k)^2);%拉力大小
if k==1
xd(k)=0-ld*cosd(thetad(k));%锚链k的x坐标 yd(k)=0-ld*sind(thetad(k));%锚链k的y坐标 else
31
xd(k)=xd(k-1)-ld*cosd(thetad(k));%锚链k的x坐标 yd(k)=yd(k-1)-ld*sind(thetad(k));%锚链k的y坐标 end
R=R+ld*cosd(thetad(k));%目前游动总长度 L=L+ld*sind(thetad(k));%
if abs(L-18)<0.1 && thetad(k)>0 && thetad(k)<30%第一二问 disp('d'); numk(t)=k; BALL(t)=ball;
THETA(t)=thetad(k); H1(t)=h1;H2(t)=h2; t=t+1; break; else
continue; end
% if thetad(k)<0
% L=L-ld*sind(thetad(k)); % if abs(L-18)<0.1 % disp('d'); % H1=h1;H2=h2; % break; % else
% continue;
% end % end
k=k+1; end
%------链环-------------------------
end end end
% plot(Td); % figure
% plot(thetad)
[thetasort,s]=sort(THETA);
H1best=H1(s(1));H2best=H2(s(1));Ballbest=BALL(s(1)); numk=numk(s(1))
32