前些天把数模做完感觉有点空虚…就没有及时的总结, 今天中秋节了, 放了三天假, 刚写了一下笔记, 感觉有点精神了, 就赶紧把国赛的总结写了吧.

正文

我们选的是A题, 总体的过程是基本一二天就到10点左右就回去了, 第二天8点钟到, 最后一天通了宵, 感觉还不算很累, 精神也还算饱满, 估计是因为这个题确实有点…不知道怎么说..说难也不难, 说简单也不简单, 我作为编程的话工作还是比较轻松的, 原理到不怎么难.

A题 系泊系统的设计

近浅海观测网的传输节点由浮标系统、系泊系统和水声通讯系统组成（如图1所示）。某型传输节点的浮标系统可简化为底面直径2m、高2m的圆柱体，浮标的质量为1000kg。系泊系统由钢管、钢桶、重物球、电焊锚链和特制的抗拖移锚组成。锚的质量为600kg，锚链选用无档普通链环，近浅海观测网的常用型号及其参数在附表中列出。钢管共4节，每节长度1m，直径为50mm，每节钢管的质量为10kg。要求锚链末端与锚的链接处的切线方向与海床的夹角不超过16度，否则锚会被拖行，致使节点移位丢失。水声通讯系统安装在一个长1m、外径30cm的密封圆柱形钢桶内，设备和钢桶总质量为100kg。钢桶上接第4节钢管，下接电焊锚链。钢桶竖直时，水声通讯设备的工作效果最佳。若钢桶倾斜，则影响设备的工作效果。钢桶的倾斜角度（钢桶与竖直线的夹角）超过5度时，设备的工作效果较差。为了控制钢桶的倾斜角度，钢桶与电焊锚链链接处可悬挂重物球。

系泊系统的设计问题就是确定锚链的型号、长度和重物球的质量，使得浮标的吃水深度和游动区域及钢桶的倾斜角度尽可能小。

问题1 某型传输节点选用II型电焊锚链22.05m，选用的重物球的质量为1200kg。现将该型传输节点布放在水深18m、海床平坦、海水密度为1.025×103kg/m3的海域。若海水静止，分别计算海面风速为12m/s和24m/s时钢桶和各节钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。
问题2 在问题1的假设下，计算海面风速为36m/s时钢桶和各节钢管的倾斜角度、锚链形状和浮标的游动区域。请调节重物球的质量，使得钢桶的倾斜角度不超过5度，锚链在锚点与海床的夹角不超过16度。
问题3 由于潮汐等因素的影响，布放海域的实测水深介于16m~20m之间。布放点的海水速度最大可达到1.5m/s、风速最大可达到36m/s。请给出考虑风力、水流力和水深情况下的系泊系统设计，分析不同情况下钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。

问题3 由于潮汐等因素的影响，布放海域的实测水深介于16m~20m之间。布放点的海水速度最大可达到1.5m/s、风速最大可达到36m/s。请给出考虑风力、水流力和水深情况下的系泊系统设计，分析不同情况下钢桶、钢管的倾斜角度、锚链形状、浮标的吃水深度和游动区域。

思路

我是编程的对建模过程就不多说了, 我们用的是集中质量法, 然后一路迭代下去, 先由浮标受力分析一直传递下去.

计算程序太多就不放上来了, 因为包括仿真求值之类的其实都是基于两个程序的基础的, 分别是2D和3D的情况下迭代, 其他模型的检验和计算都是基于这两个算法的, 后面只是把他们封装了一下

对于2D的情况, 不考虑水速,

迭代算法如下

clc,clear
%浮标参数
r=1;  %浮标底半径
h=2;   %浮标高度
mfloat=1000; %质量

g=9.8;   %重力加速度
vair=36; %风速
vsea=0; %流速
row=1025; %海水密度

theta=zeros(1,215);  %初始化各杆与水平方向夹角

m=zeros(1,215); %初始化各段质量
m([1:4])=ones(1,4)*10; % 钢管质量
mball=1200; %重物质量
m(5)=100+mball;%重物和仪器总重
m([6:215])=7*0.105;%每段缆索质量

V=zeros(1,215);
V([1:4])=pi*0.025^2*1;
V(5)=pi*0.15^2*1;
V([6:215])=0;

%d=((1000*g+sum(m(:))*g)/row/g-sum(V(:)))/pi; %吃水深度
d=0.773;
%d = 1.4

f=0.625*(h-d)*2*r*vair^2%海风载荷

T=zeros(1,215);%每段的拉力

l([1:4])=1;
l(5)=1;
l(6:215)=0.105; %各段长度

%计算各段角度
theta(1)=atan((row*pi*r^2*d*g-mfloat*g)/f)
T(1)=f/cos(theta(1))

for i=1:214
    theta(i+1)=atan((T(i)*sin(theta(i))+row*V(i)*g-m(i)*g)/(T(i)*cos(theta(i))));
    T(i+1)=(T(i)*cos(theta(i)))/cos(theta(i+1));
end

x=0;
y=0;
for i=1:215
    x=x+l(i)*cos(theta(i));
    y=y+l(i)*sin(theta(i));
end

xb=zeros(1,216);xb(1)=0;
yb=zeros(1,216);yb(1)=d-2;

for i=2:216
    if yb(i)<-18
        xb(i)=xb(i-1)-l(i-1);
        yb(i)=-18;
    else
        xb(i)=xb(i-1)-l(i-1)*cos(theta(i-1));
        yb(i)=yb(i-1)-l(i-1)*sin(theta(i-1));
    end
end

figure
plot(xb,yb)


%theta(1)*180/pi
%theta(5)*180/pi
%theta(6)*180/pi
%theta(215)*180/pi
%theta*180/pi

3d的情况下

clc,clear
%浮标参数
r=1;  %浮标底半径
h=2;   %浮标高度
mfloat=1000; %质量

g=10;   %重力加速度
vair=36;thetaair=0; %风速风向
vsea=1.5; thetasea=pi/2; %流速流向
row=1025; %海水密度
countripe=210;%缆索段数
count=countripe+5;%总的杆段数
rowripe=7;%缆索线密度
lenripe=0.105;%缆索每段长度
theta=zeros(1,count);  %初始化各杆投影角度
phi=zeros(1,count); %初始化各杆与竖直方向夹角

m=zeros(1,count); %初始化各段质量
m([1:4])=ones(1,4)*10; % 钢管质量
mball=1200; %重物质量
m(5)=100+mball;%重物和仪器总重
m([6:count])=lenripe*rowripe;%每段缆索质量

V=zeros(1,count);
V([1:4])=pi*0.025^2*1;
V(5)=pi*0.15^2*1;
V([6:count])=0;

d=((1000*g+sum(m(:))*g)/row/g-sum(V(:)))/pi; %吃水深度
%d=0.773;

f=0.625*(h-d)*2*r*vair^2;%海风载荷

T=zeros(1,count);%每段的拉力

l([1:4])=1;
l(5)=1;
l([6:count])=lenripe; %各段长度

%计算侧面积
S=zeros(1,count);
S([1:4])=0.05*1;
S(5)=0.3*1;
S([6:count])=0;

%计算各段角度
a=row*pi*r^2*d*g-mfloat*g;
b=0.625*2*r*(h-d)*vair^2*cos(thetaair)+374*2*r*d*vsea^2*cos(thetasea);
c=0.625*2*r*(h-d)*vair^2*sin(thetaair)+374*2*r*d*vsea^2*sin(thetasea);
T(1)=sqrt(a^2+b^2+c^2);
phi(1)=acos(a/T(1));
d=b/T(1)/sin(phi(1));e=c/T(1)/sin(phi(1));
if d>=0&e>=0
    theta(1)=real(acos(d));
elseif d>=0&e<0
    theta(1)=real(acos(d))-pi;
elseif d<0&e>=0
    theta(1)=real(acos(d));
else
    theta(1)=-real(acos(d));
end

for i=1:(count-1)
    a=T(i)*cos(phi(i))+row*V(i)*g-m(i)*g;
    b=T(i)*sin(phi(i))*cos(theta(i))+374*S(i)*cos(phi(i))*vsea^2;
    c=T(i)*sin(phi(i))*sin(theta(i))+374*S(i)*sin(phi(i))*vsea^2;
    T(i+1)=sqrt(a^2+b^2+c^2);
    phi(i+1)=acos(a/T(i+1));
    d=b/T(i+1)/sin(phi(i+1));e=c/T(i+1)/sin(phi(i+1));
    if d>=0&e>=0
        theta(i+1)=real(acos(d));
    elseif d>=0&e<0
        theta(i+1)=real(acos(d))-180;
    elseif d<0&e>=0
        theta(i+1)=real(acos(d));
    else
        theta(i+1)=-real(acos(d));
    end
end

x=0;
y=0;
z=d;
for i=1:count
    x=x+l(i)*sin(phi(i))*cos(theta(i));
    y=y+l(i)*sin(phi(i))*sin(theta(i));
    z=z+l(i)*cos(phi(i));
end

xb=0;
yb=0;
zb=-d;

for i=2:(count+1)
    zb(i)=zb(i-1)-l(i-1)*cos(phi(i-1));    
    if zb(i)<-18
        xb(i)=xb(i-1)-l(i-1);
        yb(i)=yb(i-1);
        zb(i)=-18;
    else
        xb(i)=xb(i-1)-l(i-1)*sin(phi(i-1))*cos(theta(i-1));
        yb(i)=yb(i-1)-l(i-1)*sin(phi(i-1))*sin(theta(i-1));
    end
end

figure
plot3(xb,yb,zb)