2016数学建模国赛A题程序（原创）作者cclplus

杨利霞

2016数学建模国赛A题程序（原创）作者cclplus

' N+ S, o& W' H- w0 ~2 V, s2 ~, Pclear all;
close all;
; s2 ?+ T  G6 R% ^" M% vclc
format long
syms h S Fw Ff Ff1 a b c d l L F depth n pl m x1 y1 y t distance n a1 b1;
/ D4 y6 u; G0 ~! MF=[];
theta=[];
$ o6 R5 N  k+ D6 K6 cv=24; %风速
, B- J( ~  {, `9 _# S: g' @3 ~& Gl=105*10^(-3); %锚链每节链环的长度
L=22.05; %锚链的总长度
num=0; %通过更改不在海床上的链节数得到一个最优解
num1=round(L/l);
num2=0;
( d8 h* v2 E6 T, mlin=0/180*pi; %第一个链节与水平方向的夹角
lin1=90/180*pi;
; C2 A6 W5 ^% c# u; u5 o/ Z3 Mlin2=0;
m2=1200; %重物球质量
  r9 `* i$ @$ f. z7 O: O# }# }# s, Zpg=7.7*10^3; %重物球的密度（单位：kg/m^3)
depth=20; %水深
" h' H, n, P. @" g2 x0 qpl=7; %锚链单位长度的质量
vh=0; %海水流速
g=9.8; %可通过改变此语句来修改重力加速度，单位为m/s^2
p=1.025*10^3; %海水密度
M=1000; %浮标质量
m=10; %钢管质量
- Z: S  F" [- qm1=100; %设备和钢桶总质量
y=0;
d=1;
j1=0;
j2=0;
9 n" F  ?- P: N% ]) R/ @while(abs(y-d)>0.005)%在这里选择所需要的精度，
, K# i* A! a  w3 X% W* Rif (y>d)&&(num<round(L/l))
num1=num;
num=round((num1+num2)/2);
8 D. Q" J4 v8 R0 v- H" z; B: Yelseif (y<d)&&(num<round(L/l));
num2=num;
& W0 h. J. j- L2 [1 _. x$ t4 {$ E+ Bnum=round((num1+num2)/2);
elseif (y<d)&&(num==round(L/l))
lin2=lin;
lin=(lin1+lin2)/2;
0 a+ L: e+ S4 k, S0 d! e& qelseif(y>d)&&(num==round(L/l))
0 S# s0 |" r2 E+ \/ x1 `lin1=lin;
5 i: F, ?$ G& G! qlin=(lin1+lin2)/2;
# s( @& G4 L9 m6 U! V# o& `) |/ ^end
%钢桶受到的浮力
Ff1=p*g*pi*(0.3/2)^2;
; j) t% V( y! R  ]2 L4 _%钢管收到的浮力
1 G) i' d8 P( \Ff2=p*g*pi*(0.05/2)^2;
%重物球所受浮力
& l4 G; I6 ?8 G* HFfg=p*g*m2/pg;
) v8 H& |( ?' [. ?/ S# X: u%重物球所受海水水流力
Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
8 J" }( B" c! MS=2*(2-h);
$ g" t4 x, t- O* J. Y%风对浮标产生的力
Fw=0.625*S*v^2;
+ E- }1 u" u+ M& C+ f5 a%浮标在水中的体积
V=pi*(2/2)^2*h;
% Y$ c5 l' M( ^# f$ q: R0 F%浮标所受到的浮力
, s( a4 i* ^( L( u- W/ IFf=p*g*V;
%浮标受到海水的近似水流力
Fb=374*2*h*vh^2;
%钢桶受到海水的近似水流力
1 b8 L0 u+ L" N9 f2 c& uFs1=374*0.3*vh^2;
%钢管受到海水的水流力的近似值
Fs=374*0.05*vh^2;
%浮标浸没水中的高度
if num==round(L/l)
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+pl*L*g+(Fhg+4*Fs+Fs1)*tan(lin))/(p*g*pi-(1.25*v^2+374*vh^2)*tan(lin));
9 T  d8 Q' }# m  \  \& n; helse
7 K  @( o, }' ]% R7 l- H$ G# jh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
9 N3 r7 z3 Z1 ^2 g3 Z) V. h3 Y) P: uend
$ r/ s+ [1 x4 y7 `6 k; r9 r# [a=Fw+Fb;
) ~" S6 ^  C# I: x% H! ib=-M*g+Ff+(Fw+Fb)*tan(lin);
if j1==0
' W! N0 Q5 d0 e, _0 V. o. ~a=eval(a);
0 \1 z" M( t% |0 [, `/ S$ d! a! Rb=eval(b);
else
6 a! u4 x1 E/ x5 }. Send
F(1)=sqrt(a^2+b^2);
7 ~5 t0 n0 J- i7 @& n. y/ stheta(1)=atan(b/a);
n=0;
for i=1:4
- j3 U4 t) y3 n" }* J0 U%钢管受到海水的水流力
/ \& p- ~/ F  U4 ?: _. r* mFh(i)=374*0.05*sin(theta(i));
: s- t0 C. o9 t0 ~/ A& l; }' ?n=n+Fh(i);
a=Fw+Fb+n;
7 E( H; ]4 W0 E- ?  \if j1==0
3 d0 \, Q3 q  O' f: M6 z$ F/ B5 Fa=eval(a);
else
end
b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
" l) w( i0 M- AF(i+1)=sqrt(a^2+b^2);
/ E- u* Q' ^  Z, S9 U% dtheta(i+1)=atan(b/a);
end
% L" S2 n  A% C2 L( l  g! U1 m6 K+ fc=0;
* X5 W4 E+ X5 h& V: ifor i=1:5
c=c+sin(theta(i));
end
d=depth-c-h;
y1=lin;
distance=0;
' e' Q; d) e7 v4 s* z/ `. {if num==round(L/l)
y=l*sin(y1);
! f$ Q5 f5 r; kx1=Fw/sqrt(1-(sin(y1))^2);
for i=1:num-1
. Q- P4 u. Q3 @m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
! x; r% c+ N- }1 l& `m=m*l;
1 ]" s4 o" C: ^: t; q, \y=y+m;
: O( T" V5 u' K& ^4 en=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
; d! f" U& G' an=eval(n);
else
end
distance=distance+n;
9 k+ P) `) K+ X# |. E0 G* R/ r* wif j1==0
9 K' \" }. M9 d2 q) {y=eval(y);
else
end
end
else
2 |- N  e# [9 ay=y1*l;
distance=(round(L/l)-num)*l;
for i=1:num
0 r1 ~' Z; m, K- L; U  X" _x1=Fw/sqrt(1-(sin(y1))^2);
* W9 u# `0 d  P" hm=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
. p$ ^  w% H0 J( p0 by=y+m;
& y, @% E' ~/ n% Tn=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
8 j, w1 H2 ?0 S1 Yy=eval(y);
) D  _# r4 I0 p- I  U6 cn=eval(n);
6 x& l) T" L- B6 C" U  ielse
end
7 N$ t( `# A; \* v& V4 U) Sdistance=distance+n;
end
* D1 [* C$ s+ M" C1 V( k8 {1 w) l: {end
m=0;
3 m0 ]3 L7 m+ S' Mj1=1;
- K2 N7 R& f; Sj2=j2+1;
end
%钢桶受到的浮力
& ]9 i; i3 v4 D4 D# t0 DFf1=p*g*pi*(0.3/2)^2;
%钢管收到的浮力
8 j: H% k2 N) Y, |& rFf2=p*g*pi*(0.05/2)^2;
%重物球所受浮力
! l" V1 z3 J' u# L8 e/ R. eFfg=p*g*m2/pg;
%重物球所受海水水流力
Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
- j' _+ K8 d; I' j5 Z  AS=2*(2-h);
; Q* x1 {% i& T; m. p8 D- b%风对浮标产生的力
$ i" N; W5 U/ E( f4 V( zFw=0.625*S*v^2;
' T% t- u( C5 y# O* L%浮标在水中的体积
, Q4 L6 u1 {0 e4 D  I$ w' qV=pi*(2/2)^2*h;
%浮标所受到的浮力
Ff=p*g*V;
%浮标受到海水的近似水流力
" e) [# F6 b5 j% f" aFb=374*2*h*vh^2;
* @0 d0 ?! {. o! l%钢桶受到海水的近似水流力
Fs1=374*0.3*vh^2;
# w" n( H: V$ I/ L%钢管受到海水的水流力的近似值
Fs=374*0.05*vh^2;
%浮标浸没水中的高度
if num==round(L/l)
' [+ Y- u2 f. R; ]& uh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+pl*L*g+(Fhg+4*Fs+Fs1)*tan(lin))/(p*g*pi-(1.25*v^2+374*vh^2)*tan(lin));
else
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
4 q6 G7 {1 W( Q$ T3 V( m0 Na=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
F(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
+ C9 `7 q3 @# @. i9 M0 B& An=0;
0 F' X9 Z7 }; g7 C) I& mfor i=1:4
%钢管受到海水的水流力
1 s9 |* L& I$ [' T4 l) W* tFh(i)=374*0.05*sin(theta(i));
! ^/ W; L) u' E% kn=n+Fh(i);
% O$ k) t/ k6 G5 D0 ta=Fw+Fb+n;
' q/ u9 _& V+ B5 A1 x( _b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
0 C% Z7 |0 }* lF(i+1)=sqrt(a^2+b^2);
# t  X; \; @7 l  H9 e, G  utheta(i+1)=atan(b/a);
end
disp('输出钢管和钢桶的倾斜角度（角度制）')
+ V1 u$ S; x# O$ Ath=90-theta*180/pi
6 G" D: c$ o/ o  E6 gm=85*pi/180;
if theta(5)>m
% R! }( I% W9 x: cdisp('钢桶的倾斜角足够小，测量准确')
else
disp('钢桶的倾斜角过大')
end
c=0;
2 g+ `. u" P$ r: A; [8 D- e; x4 zfor i=1:5
c=c+sin(theta(i));
! i+ x  u1 T9 N; [2 i  Kend
! j( h0 [, s: l& F4 I5 Fd=depth-c-h;
y1=lin;
distance=0;
5 A- Z8 g& _/ j/ q, e# K/ uif num==round(L/l)
y=l*sin(y1);
/ y. M, ]: [$ s8 R- G9 R/ `x1=Fw/sqrt(1-(sin(y1))^2);
  ^) L4 H' C2 |* C6 {. W) E+ qfor i=1:num-1
0 r1 D9 [# X9 |  Hm=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
) H/ J. O* @9 k0 ?$ Zm=m*l;
3 U$ Y; b$ g, Ky=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
: s% v( z; G# a3 p& D' O! H5 Sdistance=distance+n;
plot(distance,y,'o')
hold on
end
else
. O7 E) L( h# K' V" s8 K& ly=y1*l;
* q! q7 t8 m' G9 |$ ?  ^for i=1:round(L/l)-num
' X6 g- ~. V  T2 M6 Xdistance=i*l;
( W. g' g, a0 q2 I" ]: ay=0;
plot(distance,y,'o')
hold on
grid on
end
for i=1:num
( k9 Y# f3 b5 Y! C: @x1=Fw/sqrt(1-(sin(y1))^2);
m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
& ^' y& J: N8 Z6 L' u: @" Zy=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
8 k3 P: d5 t7 s0 A4 @+ Y/ b- ~" pif j1==0
y=eval(y);
0 Z  c* x- _* D  Sn=eval(n);
else
end
" \# w( a+ i* odistance=distance+n;
5 Z! l- S, e9 C) a( e6 i: C, Cplot(distance,y,'o')
hold on
end
end
m=0;
for i=1:5
m=m+cos(theta(i));
# H, B6 X' G3 Y4 J% S8 J2 Zend
%浮标的运动半径
disp('输出浮标的运动半径')
' n. P* c( @% M% R3 Dans=distance+m
( k7 S: Q6 F: K



1 U! G1 k0 f& ]8 I, z8 ]8 p/ a
4 h2 W$ b! H+ D
, J2 ?; n  y# q+ d" D
1 |; ]8 M4 g1 I7 c  e

571334077

2333333333333333333333
: F1 e" u5 t$ g  l9 n( f/ S