2016年全国大学生数学建模大赛国家一等奖优秀论文,系泊系统的设计 下载本文

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