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

杨利霞

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

1 P4 g$ m/ K4 r* g

' r) i+ P# ]" z# R( ~clear all;
close all;
clc
9 P$ v/ H" C0 [' J: gformat long
5 W: ], B. P5 I8 H3 V, }1 x: Ksyms h S Fw Ff Ff1 a b c d l L F depth n pl m x1 y1 y t distance n a1 b1;
F=[];
% J5 C+ F* z; I: Ctheta=[];
2 l8 M( t# Q4 l4 Y' |( G8 D: [v=24; %风速
. m- c  x+ O. M6 a  @4 g/ ol=105*10^(-3); %锚链每节链环的长度
L=22.05; %锚链的总长度
num=0; %通过更改不在海床上的链节数得到一个最优解
; c% V, M. e6 ~" znum1=round(L/l);
) T3 j' P- ^/ I. h# T. C& h* \num2=0;
lin=0/180*pi; %第一个链节与水平方向的夹角
lin1=90/180*pi;
lin2=0;
m2=1200; %重物球质量
pg=7.7*10^3; %重物球的密度（单位：kg/m^3)
+ j9 @+ ?# P$ idepth=20; %水深
) i* c; R* k- G  w/ F1 E0 Qpl=7; %锚链单位长度的质量
vh=0; %海水流速
g=9.8; %可通过改变此语句来修改重力加速度，单位为m/s^2
# j! i+ j3 I2 bp=1.025*10^3; %海水密度
5 f  v8 v  ^/ A0 D# VM=1000; %浮标质量
m=10; %钢管质量
m1=100; %设备和钢桶总质量
4 d% a: K' \+ Q& f0 Hy=0;
7 h% i' K1 k# f* }5 [0 wd=1;
7 U0 _0 E3 ~8 N0 B% u! jj1=0;
j2=0;
9 I3 l9 b5 [. j5 i0 ?& gwhile(abs(y-d)>0.005)%在这里选择所需要的精度，
7 B$ t& L- S9 ], P2 G# Y; g; Y8 dif (y>d)&&(num<round(L/l))
$ D! h& |9 _0 h3 i6 T0 b) `# B2 Tnum1=num;
9 k. n1 H- v; u# ~( y  ^num=round((num1+num2)/2);
; h6 X- j2 X* E0 E  g: h- Helseif (y<d)&&(num<round(L/l));
num2=num;
4 ~/ _; ^% m0 z$ c9 R' K! dnum=round((num1+num2)/2);
elseif (y<d)&&(num==round(L/l))
lin2=lin;
5 l; _5 A+ G6 x2 L8 u' }5 I6 wlin=(lin1+lin2)/2;
6 i+ w" K( B( p+ Q2 velseif(y>d)&&(num==round(L/l))
5 e0 G# X9 t6 m2 V$ }0 }8 Ilin1=lin;
lin=(lin1+lin2)/2;
4 z  P. f- t7 h1 H$ C' Wend
, H8 N/ W. Y, w& u- X! k%钢桶受到的浮力
; T* W( x3 n; R. E& BFf1=p*g*pi*(0.3/2)^2;
( i( }4 U$ Y2 K7 J%钢管收到的浮力
% S$ q' J( Y, D' {5 NFf2=p*g*pi*(0.05/2)^2;
%重物球所受浮力
Ffg=p*g*m2/pg;
%重物球所受海水水流力
  l1 Y9 O- l9 u  i! X. _; _Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
S=2*(2-h);
: e' V  S! R# d" X%风对浮标产生的力
Fw=0.625*S*v^2;
%浮标在水中的体积
7 `4 |7 P4 V$ ^+ @" y- rV=pi*(2/2)^2*h;
' Y( j; P  o* b2 _+ j%浮标所受到的浮力
Ff=p*g*V;
%浮标受到海水的近似水流力
/ ]  @& ?/ L8 W: J. FFb=374*2*h*vh^2;
%钢桶受到海水的近似水流力
Fs1=374*0.3*vh^2;
  b- c) E, X. ]1 H. d%钢管受到海水的水流力的近似值
Fs=374*0.05*vh^2;
- _* l0 @) K  e0 s& B%浮标浸没水中的高度
) s1 m& P% x3 }% B: O) Vif 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));
else
+ D1 w4 v* S9 x4 ~. mh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
% B/ b- p& v+ Mend
a=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
if j1==0
1 P2 [+ f7 [) Q" V7 W/ [$ ^# Ha=eval(a);
b=eval(b);
else
2 t/ }- [8 R0 k% |" M: }) R* Xend
3 H4 a; C' c7 V% x7 |$ O6 X' d* X! zF(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
  q& f. `% S: Hn=0;
% |0 f$ j# X; n" M$ N9 z; Efor i=1:4
%钢管受到海水的水流力
% c" r$ r( {  d" v& P" l, g: hFh(i)=374*0.05*sin(theta(i));
! _& \$ d; \  Zn=n+Fh(i);
a=Fw+Fb+n;
* I0 w6 V- A$ X& h: {if j1==0
a=eval(a);
else
end
: _$ H. ]% ]' z$ Eb=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
c=0;
9 y- n. r- M: c# U/ L3 M  \for i=1:5
c=c+sin(theta(i));
7 M2 ?. x( k; Wend
8 C2 g% t) B- g5 I1 r  K) Pd=depth-c-h;
y1=lin;
distance=0;
if num==round(L/l)
y=l*sin(y1);
x1=Fw/sqrt(1-(sin(y1))^2);
  O9 }3 c0 A% o# r: wfor i=1:num-1
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
% O$ \, ~% @; H( O: {9 v1 S; A$ om=m*l;
y=y+m;
! S# D! y: M3 A1 U: Mn=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
n=eval(n);
else
end
' a; q6 p  r9 z5 _5 A4 Tdistance=distance+n;
& [2 @( ]6 P5 A% P$ ?4 W/ tif j1==0
( R' E& [* d# Q0 C" q+ uy=eval(y);
else
end
end
else
2 g1 j* H  R' r/ j* My=y1*l;
; q$ R, a; D( W8 P- r' D% W- `distance=(round(L/l)-num)*l;
for i=1:num
; B9 z# ~; H( v: b  s4 f1 i3 ox1=Fw/sqrt(1-(sin(y1))^2);
m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
/ Y5 P1 C6 Z" {1 G/ @y=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
  z4 R# D0 p$ H% c9 z3 }if j1==0
y=eval(y);
n=eval(n);
else
end
distance=distance+n;
end
end
6 _9 v4 z# l  c+ N1 \m=0;
j1=1;
& k6 e. E( [; ~2 P- v, U5 t7 Oj2=j2+1;
5 R4 b; n/ t; o) l9 B8 pend
% h6 W  N( [: m& {/ e- G; m%钢桶受到的浮力
Ff1=p*g*pi*(0.3/2)^2;
; y0 s$ A% _9 G2 X) N%钢管收到的浮力
Ff2=p*g*pi*(0.05/2)^2;
%重物球所受浮力
Ffg=p*g*m2/pg;
! t( |# j" A. S: Y% m%重物球所受海水水流力
. ]! D6 u- q* I5 C0 uFhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%风对浮标受力面的投影面积
0 e, Z6 D4 x* H* kS=2*(2-h);
; j9 A! ^5 E1 Z5 a( n) D4 J4 w, u%风对浮标产生的力
7 I# L2 k: ~5 d# ZFw=0.625*S*v^2;
; |5 q0 h: v  j* |# T7 D%浮标在水中的体积
V=pi*(2/2)^2*h;
5 N. W7 Y, Q; F+ p: d9 u%浮标所受到的浮力
Ff=p*g*V;
& K$ Z2 D; N4 H2 `9 v; V7 n8 d%浮标受到海水的近似水流力
Fb=374*2*h*vh^2;
! Y2 N- O% m4 Q6 C%钢桶受到海水的近似水流力
2 s* ]8 u$ a4 w' NFs1=374*0.3*vh^2;
# o; B# A! F7 K% E- u%钢管受到海水的水流力的近似值
; C+ ~8 w: o0 v+ DFs=374*0.05*vh^2;
%浮标浸没水中的高度
if num==round(L/l)
2 b6 M* _& L) k4 W/ ^' Qh=(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
4 K! }) W/ B/ {8 A6 dh=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
a=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
F(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
1 I& J* z: @  f! q# N4 On=0;
for i=1:4
2 v/ w! W- Q+ r%钢管受到海水的水流力
Fh(i)=374*0.05*sin(theta(i));
9 e3 y2 V, B( ?# in=n+Fh(i);
+ i. B- v( I: z7 S0 Q8 r1 Qa=Fw+Fb+n;
b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
/ w2 k( `& @% ?F(i+1)=sqrt(a^2+b^2);
2 G) ?* Y2 A- V+ Ytheta(i+1)=atan(b/a);
end
disp('输出钢管和钢桶的倾斜角度（角度制）')
/ [+ [& f9 }2 V& C: Hth=90-theta*180/pi
% A! v  }% I7 mm=85*pi/180;
if theta(5)>m
- n; k! V4 a) v$ h7 Edisp('钢桶的倾斜角足够小，测量准确')
8 c( p: O  b' v+ p$ M5 telse
* G0 X# l& [# A8 mdisp('钢桶的倾斜角过大')
' c8 d3 U  K" _/ D( J0 fend
c=0;
for i=1:5
c=c+sin(theta(i));
& F. E0 C! Q  s, ?: X1 c% C0 Hend
d=depth-c-h;
y1=lin;
' F# @0 s# F8 Wdistance=0;
if num==round(L/l)
- B4 Q* f: w  v7 ?& g: m& ey=l*sin(y1);
x1=Fw/sqrt(1-(sin(y1))^2);
; Q, N5 X5 A4 Bfor i=1:num-1
$ H4 @# \# r' O8 K. I9 Pm=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
; E5 I$ o  }# rm=m*l;
y=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
6 p* J+ a8 b- z% f4 c6 m; n) jdistance=distance+n;
( t6 M3 t- ~  Z" gplot(distance,y,'o')
hold on
& ~8 e3 g3 o1 l: t' z6 F$ A; ?2 _end
else
& p" C0 i) v# Y4 ny=y1*l;
. c6 S# D3 Y1 c  E% p+ F! O  kfor i=1:round(L/l)-num
distance=i*l;
& M  j) C! a2 Dy=0;
) \5 q5 P0 r' _, o# ]plot(distance,y,'o')
hold on
3 h) {. w7 K: c* L. V) i  dgrid on
8 K$ u, f4 w( [' d2 D9 E7 Cend
for i=1:num
x1=Fw/sqrt(1-(sin(y1))^2);
; O4 y& A9 C, ~  @. ^m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
3 [( g) L) Q: z) |5 j0 q$ ay=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
+ W4 s: A1 Q. w- ~if j1==0
1 N: V' Z$ u; M# {y=eval(y);
' q4 u5 `' Z$ i+ m0 n0 en=eval(n);
else
end
distance=distance+n;
" T* P: {$ N# c8 C2 y+ A; X  _plot(distance,y,'o')
hold on
1 i8 D' M/ Y$ B) t. S# {end
end
m=0;
5 ~- e' r  `- ^( i' l$ c% ?for i=1:5
( t, b+ ?: c. E9 {; P; Um=m+cos(theta(i));
  X, W" z. z; Q/ |7 B& ?" Iend
% o; R* U1 ^6 v6 ~' L4 d' ?%浮标的运动半径
( x8 g  y9 o& I, mdisp('输出浮标的运动半径')
ans=distance+m

9 [! k: O5 m! F* N, e0 d0 f
7 F/ Y) ^! S$ x% I1 `
5 Y, q) s* Q* |; w2 T7 j& h
; ^7 G9 C; L* J0 Q. A" M( m: F! o

571334077

2333333333333333333333