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

杨利霞

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

# M4 f- D  W' x" gclear all;
% f2 L' W0 F8 x0 H. _. I9 N. gclose all;
. x* E1 S4 I: ]: g, s  n) p( pclc
/ E, D( Y& U% `% Tformat long
; ~0 T$ J) f% e" K" G; R4 vsyms h S Fw Ff Ff1 a b c d l L F depth n pl m x1 y1 y t distance n a1 b1;
. e0 o) S8 r( [' tF=[];
theta=[];
v=24; %风速
  P( c4 M7 K: y! x0 Rl=105*10^(-3); %锚链每节链环的长度
L=22.05; %锚链的总长度
num=0; %通过更改不在海床上的链节数得到一个最优解
num1=round(L/l);
: k6 k! E& `4 s9 |0 Ynum2=0;
: a% m, @2 v! ^' A0 Mlin=0/180*pi; %第一个链节与水平方向的夹角
, j4 _* r5 j  C7 n. Ilin1=90/180*pi;
lin2=0;
% a! v" J+ ^2 Z& B6 d* O2 Pm2=1200; %重物球质量
8 i$ r. F8 O1 V8 S+ Q( k9 O. Bpg=7.7*10^3; %重物球的密度（单位：kg/m^3)
depth=20; %水深
pl=7; %锚链单位长度的质量
+ F: J2 W( Z; W2 jvh=0; %海水流速
g=9.8; %可通过改变此语句来修改重力加速度，单位为m/s^2
p=1.025*10^3; %海水密度
M=1000; %浮标质量
; S* s1 {; ^, y! g1 `, Vm=10; %钢管质量
m1=100; %设备和钢桶总质量
y=0;
4 L! t% c+ l, x& ]( E" f: }  ?. l1 u: Z1 Xd=1;
2 K+ {& H( E1 P$ Dj1=0;
: B, D- Z9 [' |7 Q5 ^  p9 C6 Nj2=0;
while(abs(y-d)>0.005)%在这里选择所需要的精度，
if (y>d)&&(num<round(L/l))
9 Q" j  T5 p( \# h% H6 d8 vnum1=num;
" C2 [0 _0 t2 I- P& W5 ]$ Onum=round((num1+num2)/2);
elseif (y<d)&&(num<round(L/l));
num2=num;
num=round((num1+num2)/2);
elseif (y<d)&&(num==round(L/l))
lin2=lin;
lin=(lin1+lin2)/2;
elseif(y>d)&&(num==round(L/l))
lin1=lin;
' Y/ f# P" Z2 K4 S7 s' tlin=(lin1+lin2)/2;
end
. r. z5 v# d; l: ]) o# Z%钢桶受到的浮力
1 h2 L7 q# l! `# ]/ m! WFf1=p*g*pi*(0.3/2)^2;
%钢管收到的浮力
Ff2=p*g*pi*(0.05/2)^2;
. }" T7 p7 }9 C0 B! }% M3 J& a%重物球所受浮力
' K7 w$ n& C) G# J8 RFfg=p*g*m2/pg;
%重物球所受海水水流力
Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
S=2*(2-h);
; A' r5 N3 e4 ~, v& h%风对浮标产生的力
9 \9 \( `1 |, F% @1 PFw=0.625*S*v^2;
%浮标在水中的体积
. ^+ i7 y5 r- B! P. R) b1 mV=pi*(2/2)^2*h;
%浮标所受到的浮力
) J0 C6 q- m) d$ r1 dFf=p*g*V;
%浮标受到海水的近似水流力
Fb=374*2*h*vh^2;
%钢桶受到海水的近似水流力
Fs1=374*0.3*vh^2;
%钢管受到海水的水流力的近似值
% @/ z9 N/ l  vFs=374*0.05*vh^2;
1 S$ u3 T; M; q) Y! E" c%浮标浸没水中的高度
if num==round(L/l)
% b  E8 P. A+ u# r3 E- jh=(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));
, U. _3 j3 T5 S( j( Xelse
2 Q, _& i9 C( g* o; z6 Uh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
9 ^( _% F5 Y8 p$ f7 k, ga=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
if j1==0
a=eval(a);
b=eval(b);
else
end
F(1)=sqrt(a^2+b^2);
0 z- ^: n' p+ g; P" V$ T2 ^9 atheta(1)=atan(b/a);
n=0;
( b2 G3 b2 i+ p$ u8 zfor i=1:4
%钢管受到海水的水流力
Fh(i)=374*0.05*sin(theta(i));
+ k$ V: f& ^0 d& G# In=n+Fh(i);
a=Fw+Fb+n;
if j1==0
a=eval(a);
else
0 f& n% U1 G. q: q* cend
/ I/ B4 ]& y3 g) Q; Jb=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
1 x+ m! X7 ~# U! ~6 H% D& mF(i+1)=sqrt(a^2+b^2);
# ^+ d! F( o$ k0 I/ {theta(i+1)=atan(b/a);
end
- ~8 O  `. e* [5 C' Q1 Y6 Ec=0;
for i=1:5
c=c+sin(theta(i));
end
$ L9 K4 m5 E  W8 Jd=depth-c-h;
) a( w' }+ z! T5 ^9 O# W! l" l9 ty1=lin;
distance=0;
& u) C( b  ?4 a( Vif num==round(L/l)
y=l*sin(y1);
4 u1 [, ^: |, N) {8 M1 Rx1=Fw/sqrt(1-(sin(y1))^2);
for i=1:num-1
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
m=m*l;
y=y+m;
9 A( B! F; X8 p3 @n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
n=eval(n);
6 W) G# k+ N( ~: @else
end
, o3 [' W' f( gdistance=distance+n;
3 y' d" D  |' `  B- i) @& `) cif j1==0
/ \9 Q0 Y* A0 Z' p! \3 dy=eval(y);
else
end
3 V" \- M# [# A' k  Y' L" gend
else
y=y1*l;
( u; |- K( o5 |3 I- L9 udistance=(round(L/l)-num)*l;
for i=1:num
x1=Fw/sqrt(1-(sin(y1))^2);
/ f- ?. p1 W* W! Om=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
  M  i( P/ d. |* b% j- Iy=y+m;
' ]  h6 _1 ~3 ^( [* Q* Z  ~n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
y=eval(y);
" I* S: {; [  `* j- W: `n=eval(n);
else
end
3 L# |' J0 I- b# tdistance=distance+n;
( C, M' A! U- z% B0 kend
7 }& B' |; S; s% O6 E+ F5 Mend
( b" Y6 O7 D1 W4 E. s# v  Rm=0;
j1=1;
! p, v: t. {5 [j2=j2+1;
end
%钢桶受到的浮力
Ff1=p*g*pi*(0.3/2)^2;
%钢管收到的浮力
0 ?0 Y9 q! D9 f: IFf2=p*g*pi*(0.05/2)^2;
9 q5 I' n; t) b1 }' o, V%重物球所受浮力
, R# V5 T/ `7 R" lFfg=p*g*m2/pg;
%重物球所受海水水流力
# t4 ]/ g0 A3 J" [Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
; R; b7 {' e! o; Y# C* \3 y%风对浮标受力面的投影面积
S=2*(2-h);
0 q1 X+ X( c7 i1 C1 f* ~3 L4 a/ S' f5 E  d%风对浮标产生的力
Fw=0.625*S*v^2;
7 D) |* V) y; H6 N%浮标在水中的体积
V=pi*(2/2)^2*h;
# t) o% `1 y) Y* f%浮标所受到的浮力
6 B5 r* l0 W+ lFf=p*g*V;
%浮标受到海水的近似水流力
1 H4 R/ x% D, h6 r  `Fb=374*2*h*vh^2;
%钢桶受到海水的近似水流力
5 f5 N$ o: E+ G1 sFs1=374*0.3*vh^2;
%钢管受到海水的水流力的近似值
# V+ m5 v2 a1 y) j5 z  k' iFs=374*0.05*vh^2;
%浮标浸没水中的高度
if num==round(L/l)
9 d3 e- P! I+ I3 N9 b5 {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));
& }6 P7 @! \, D- k) _/ Uelse
) I" w- b7 W! Z  e: n+ Xh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
' Z4 z' F7 d7 Z. {5 m6 Rend
# Q" ]6 Y) n$ b9 N) ]+ da=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
F(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
n=0;
for i=1:4
) ~* S4 ~) ?6 A- C# D5 S5 L- m%钢管受到海水的水流力
0 V2 a* P; }9 u7 F+ ]Fh(i)=374*0.05*sin(theta(i));
n=n+Fh(i);
a=Fw+Fb+n;
$ x, |7 Q) K# V+ r# `: {  wb=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
F(i+1)=sqrt(a^2+b^2);
theta(i+1)=atan(b/a);
end
disp('输出钢管和钢桶的倾斜角度（角度制）')
/ ?0 Y8 U: m  Z! Xth=90-theta*180/pi
m=85*pi/180;
- ^5 C' ?$ L" Iif theta(5)>m
disp('钢桶的倾斜角足够小，测量准确')
' h2 T- J' k# x" v3 ^% ?1 Velse
3 Q' A  E# o5 J8 N* Bdisp('钢桶的倾斜角过大')
end
c=0;
7 W' [" j" C' }, p; k; afor i=1:5
4 p) s% k3 y  I+ {! Sc=c+sin(theta(i));
end
7 d' p1 A1 K7 ~+ }3 k- {0 s- Pd=depth-c-h;
$ s" R  w+ h8 y+ x% Wy1=lin;
! }; v  A$ ]  n7 L3 @distance=0;
, z( _/ u! ^- k* Y7 y! nif num==round(L/l)
y=l*sin(y1);
x1=Fw/sqrt(1-(sin(y1))^2);
for i=1:num-1
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
4 C. Q5 |& w, R- D' x3 dm=m*l;
) z# d; Z1 `9 ay=y+m;
) o- p; R+ ]" l8 j5 x, `' D5 an=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
distance=distance+n;
plot(distance,y,'o')
hold on
end
2 d- ~) j8 N2 x# ]  e9 L9 v/ g- ielse
: K$ v; L5 W9 b) }7 j0 my=y1*l;
for i=1:round(L/l)-num
distance=i*l;
y=0;
plot(distance,y,'o')
hold on
$ Y5 Y$ h) q4 s* o" ^1 z# qgrid on
end
1 L, R  f$ |, J* x  ^0 Y* e, q) ~for i=1:num
x1=Fw/sqrt(1-(sin(y1))^2);
, m7 I% W: V6 D2 K* ~# H! O# _m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
4 Y+ K0 b" |! V  ^; N! By=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
8 n! j2 D4 [& z0 d4 J2 }8 z$ zif j1==0
y=eval(y);
n=eval(n);
else
end
distance=distance+n;
) A7 n; a  _* m  {4 o/ B- {/ W( Fplot(distance,y,'o')
2 W5 y9 @9 R$ C& xhold on
end
8 [' g3 q/ s7 i1 `5 Qend
m=0;
for i=1:5
m=m+cos(theta(i));
& @: a7 \! [- F: i& X/ p& mend
%浮标的运动半径
- x6 s! u; s5 Sdisp('输出浮标的运动半径')
/ f- ]& i3 e: h5 v2 _( ians=distance+m
# M( Q5 r9 I' R# d# o' S
/ G7 T6 s9 u4 N  S5 V) s1 F  |
/ }1 K+ H$ g5 D# y# ~
' L; P" B: g8 G$ e' Y2 ]+ z. Z* H


) v7 R' ^3 @1 r8 w  w
+ D( }  ]& L7 v

571334077

2333333333333333333333